2014-04-01 74 views
2

我試圖實現在ř(如在this whitepaper描述)其涉及解決與兩個未知參數的式子的β-幾何概率模型兩個參數。在該示例中,他們使用Excel來執行此操作,從alpha = beta = 1開始將值限制爲alpha > 0.0001 < beta求解最大似然與下約束

我幾乎在R中實現了這一點,但我似乎無法爲我解決任何求解器問題。請幫忙。

RFiddle here

# 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) 

回答

2

默認情況下,optim將最大限度地減少,但是你想最大化LL。此外,要使用像L-BFGS-B,它使用綁定信息的方法:

optim(c(1, 1), function(x) -LL(x[1], x[2], n=n, t=t), method="L-BFGS-B", 
     lower=c(0.0001, 0.0001)) 
# $par 
# [1] 0.0001 9679.3562 
# 
# $value 
# [1] 17075.64 
# 
# $counts 
# function gradient 
#  87  87 
# 
# $convergence 
# [1] 0 
# 
# $message 
# [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" 

我們可以驗證我們已經改善了數似然:

LL(1, 1, n=n, t=t) 
# [1] -27659.45 

LL(0.0001, 9679.3562, n=n, t=t) 
# [1] -17075.64