2014-06-27 46 views
0

以下代码是查找广义β分布的3个参数。optimizing()函数参数估计误差

beta1 = function(a,b,t) { beta(a+(1/t),b) } 
beta2 = function(a,b,t) { beta(a+(2/t),b) } 

geb11 = function(a,b,t) { beta2(a,b,t)/beta(a,b) } 
geb12 = function(a,b,t) { (beta1(a,b,t)-beta2(a,b,t))/beta(a,b) } 
geb22 = function(a,b,t) { 1 + (beta2(a,b,t)-2*beta1(a,b,t))/beta(a,b) } 

gbetloglik = function(v) { 
    a = v[1]; b = v[2]; t = v[3] 
    loglik = n1*log(geb11(a,b,t)) + n3*log(geb12(a,b,t)) + n2*log(geb22(a,b,t)) 
    return(-loglik) 
} 

n1 = 127; n2 = 111; n3 = 262 
abt = optim(c(5,5,1),gbetloglik,lower=c(0.01,0.01,0.1),method="L-BFGS-B")$par 

然后它给出了一个错误,如下面的许多警告。

Error in optim(c(2, 8, 1), gbetloglik, lower = c(0.01, 0.01, 0.1), method = "L-BFGS-B") : 
    L-BFGS-B needs finite values of 'fn' 

警告消息是

1: In beta(a + (2/t), b) : underflow occurred in 'beta' 
2: In beta(a, b) : underflow occurred in 'beta' 
3: In beta(a + (1/t), b) : underflow occurred in 'beta' 
... 

N1,N2和N3是由一些随机数生成场景生成,和(N1,N2,N3)判定是否发生了错误。 (对于某些(N1,N2,N3)不会发生错误)

我不明白,这个错误的原因是什么,什么是“溢在‘测试版’的意思。

回答

1

作警告消息计算表明公测时有下溢

试试这个(我们已经添加了##行):

gbetloglik = function(v) { 
    a = v[1]; b = v[2]; t = v[3] 
    if (beta(a, b) == 0) return(1e6) ## 
    loglik = n1*log(geb11(a,b,t)) + n3*log(geb12(a,b,t)) + n2*log(geb22(a,b,t)) 
    return(-loglik) 
} 
+0

看来,并称命令后不经常发生的错误,但它仍然存在.. – user67275

+0

然后适当地处理beta的其他事件。 –

+0

你的建议解决了'错误'。但是仍然有'beta'警告下溢 – user67275