2
我試圖實現在ř(如在this whitepaper描述)其涉及解決與兩個未知參數的式子的β-幾何概率模型兩個參數。在該示例中,他們使用Excel來執行此操作,從alpha = beta = 1
開始將值限制爲alpha > 0.0001 < beta
。求解最大似然與下約束
我幾乎在R中實現了這一點,但我似乎無法爲我解決任何求解器問題。請幫忙。
# probability mass function
P = function (alpha, beta, t) {
out = numeric(length=length(t))
for (i in seq_along(t)) {
if (t[i]==1) {
out[i] = alpha/(alpha + beta)
} else {
out[i] = ((beta + t[i] - 2)/(alpha + beta + t[i] - 1)) * P(alpha, beta, t[i] - 1)
}
}
out
}
# survival probability function
S = function(alpha, beta, t) {
if (t==1) {
1 - P(alpha, beta, t=t)
} else {
S(alpha, beta, t - 1) - P(alpha, beta, t=t)
}
}
# log likelihood function
LL = function(alpha, beta=1, t, n) {
sum(n * log(P(1,1,t))) + (sum(n[1:length(n)-1]) * log(S(alpha, beta, t=t[length(t)])))
}
# make some test data
n = c(239L, 2650L, 1063L, 1014L, 473L, 1304L)
t = 1:6
# log likelihood
LL(alpha=1, beta=1, n=n, t=t)
# use numerical optimization to find values of alpha and beta
optim(c(1,1), fn=LL, n=n, t=t)
require(stats4)
mle(LL, start=list(alpha=1, beta=1), t=t, n=n)