如何解决R 函数产生 NAN，诊断统计::优化
我运行此代码是为了计算零膨胀模型中 Beta 二项式分布的 MLE。 但是，对于特定数据，我遇到了一些错误。请你帮帮我好吗？ 这是 R 代码：
type = "zi"; lowerbound = 0.01; upperbound = 10000;
n=27.59554;alpha1=19.22183;alpha2=41.90441;
x=c(9,12,13,11,7,6,2,9,10,11)
N = length(x)
t = x[x > 0]
m = length(t)
neg.log.lik <- function(y)
{
n1 = y[1]
a1 = y[2]
b1 = y[3]
logA = lgamma(a1 + n1 + b1) + lgamma(b1)
logB = lgamma(a1 + b1) + lgamma(n1 + b1)
ans = m * log(1 - exp(logB - logA)) + m * logA - m *
lgamma(n1 + 1) - sum(lgamma(t + a1)) - sum(lgamma(n1 -
t + b1)) - m * lgamma(a1 + b1) + sum(lgamma(t + 1)) +
sum(lgamma(n1 - t + 1)) + m * lgamma(a1)
return(ans)
}
gp <- function(y)
{
#n1=27.59554;a1=19.22183;b1=41.90441;
n1 = y[1]
a1 = y[2]
b1 = y[3]
logA = lgamma(a1 + n1 + b1) + lgamma(b1)
logB = lgamma(a1 + b1) + lgamma(n1 + b1)
dn = -m * exp(logB - logA) * (digamma(n1 + b1) - digamma(a1 +
n1 + b1))/(1 - exp(logB - logA)) - m * digamma(n1 +
1) - sum(digamma(n1 + b1 - t)) + sum(digamma(n1 -
t + 1)) + m * digamma(a1 + n1 + b1)
da = -m * exp(logB - logA) * (digamma(a1 + b1) - digamma(a1 +
n1 + b1))/(1 - exp(logB - logA)) - sum(digamma(t +
a1)) - m * digamma(a1 + b1) + m * digamma(a1 + n1 +
b1) + m * digamma(a1)
db = -m * exp(logB - logA) * (digamma(a1 + b1) + digamma(n1 +
b1) - digamma(a1 + n1 + b1) - digamma(b1))/(1 - exp(logB -
logA)) + m * digamma(b1) - sum(digamma(n1 - t + b1))-
m * digamma(a1 + b1) + m * digamma(a1 + n1 + b1)
return(c(dn,da,db))
}
estimate = stats::optim(par = c(n,alpha1,alpha2),fn = neg.log.lik,gr = gp,method = "L-BFGS-B",lower = c(max(x) - lowerbound,lowerbound,lowerbound),upper = c(upperbound,upperbound,upperbound))
错误是： stats::optim(par = c(n,中的错误： 由 optim 提供的非有限值 另外： 警告信息： 在 digamma(n1 + b1 - t) 中：产生 NaN
你有什么想法吗？欣赏任何建议
解决方法
问题在于，由于 digamma(0)，dn 在某些时候会变成 NaN。
一个选择是像我一样考虑这种可能性。但是您应该探索在该事件中该怎么做。
您在这里遇到了问题 digamma(n1 + b1 - t) 在某些时候它会产生 digamma(0)，所以 NaN 所以失败。
type = "zi"; lowerbound = 0.01; upperbound = 10000;
n=27.59554;alpha1=19.22183;alpha2=41.90441;
x=c(9,12,13,11,7,6,2,9,10,11)
N = length(x)
t = x[x > 0]
m = length(t)
neg.log.lik <- function(y)
{
n1 = y[1]
a1 = y[2]
b1 = y[3]
logA = lgamma(a1 + n1 + b1) + lgamma(b1)
logB = lgamma(a1 + b1) + lgamma(n1 + b1)
ans = m * log(1 - exp(logB - logA)) + m * logA - m *
lgamma(n1 + 1) - sum(lgamma(t + a1)) - sum(lgamma(n1 -
t + b1)) - m * lgamma(a1 + b1) + sum(lgamma(t + 1)) +
sum(lgamma(n1 - t + 1)) + m * lgamma(a1)
return(ans)
}
gp <- function(y)
{
#n1=27.59554;a1=19.22183;b1=41.90441;
n1 = y[1]
a1 = y[2]
b1 = y[3]
logA = lgamma(a1 + n1 + b1) + lgamma(b1)
logB = lgamma(a1 + b1) + lgamma(n1 + b1)
dn = -m * exp(logB - logA) * (digamma(n1 + b1) - digamma(a1 +
n1 + b1))/(1 - exp(logB - logA)) - m * digamma(n1 +
1) - sum(digamma(n1 + b1 - t)) + sum(digamma(n1 -
t + 1)) + m * digamma(a1 + n1 + b1)
if(is.na(dn)){
dn=-99999999
}
da = -m * exp(logB - logA) * (digamma(a1 + b1) - digamma(a1 +
n1 + b1))/(1 - exp(logB - logA)) - sum(digamma(t +
a1)) - m * digamma(a1 + b1) + m * digamma(a1 + n1 +
b1) + m * digamma(a1)
db = -m * exp(logB - logA) * (digamma(a1 + b1) + digamma(n1 +
b1) - digamma(a1 + n1 + b1) - digamma(b1))/(1 - exp(logB -
logA)) + m * digamma(b1) - sum(digamma(n1 - t + b1))-
m * digamma(a1 + b1) + m * digamma(a1 + n1 + b1)
print("dn")
print(dn)
print("da")
print(da)
print("db")
print(db)
return(c(dn,da,db))
}
estimate = stats::optim(par = c(n,alpha1,alpha2),fn = neg.log.lik,gr = gp,method = "L-BFGS-B",lower = c(max(x) - lowerbound,lowerbound,lowerbound),upper = c(upperbound,upperbound,upperbound))
另一种选择是添加少量：digamma(n1 + b1 - t)+0.001
type = "zi"; lowerbound = 0.01; upperbound = 10000;
n=27.59554;alpha1=19.22183;alpha2=41.90441;
x=c(9,11)
N = length(x)
t = x[x > 0]
m = length(t)
neg.log.lik <- function(y)
{
n1 = y[1]
a1 = y[2]
b1 = y[3]
logA = lgamma(a1 + n1 + b1) + lgamma(b1)
logB = lgamma(a1 + b1) + lgamma(n1 + b1)
ans = m * log(1 - exp(logB - logA)) + m * logA - m *
lgamma(n1 + 1) - sum(lgamma(t + a1)) - sum(lgamma(n1 -
t + b1)) - m * lgamma(a1 + b1) + sum(lgamma(t + 1)) +
sum(lgamma(n1 - t + 1)) + m * lgamma(a1)
return(ans)
}
gp <- function(y)
{
#n1=27.59554;a1=19.22183;b1=41.90441;
n1 = y[1]
a1 = y[2]
b1 = y[3]
logA = lgamma(a1 + n1 + b1) + lgamma(b1)
logB = lgamma(a1 + b1) + lgamma(n1 + b1)
dn = -m * exp(logB - logA) * (digamma(n1 + b1) - digamma(a1 +
n1 + b1))/(1 - exp(logB - logA)) - m * digamma(n1 +
1) - sum(digamma(n1 + b1 - t)+0.001) + sum(digamma(n1 -
t + 1)) + m * digamma(a1 + n1 + b1)
da = -m * exp(logB - logA) * (digamma(a1 + b1) - digamma(a1 +
n1 + b1))/(1 - exp(logB - logA)) - sum(digamma(t +
a1)) - m * digamma(a1 + b1) + m * digamma(a1 + n1 +
b1) + m * digamma(a1)
db = -m * exp(logB - logA) * (digamma(a1 + b1) + digamma(n1 +
b1) - digamma(a1 + n1 + b1) - digamma(b1))/(1 - exp(logB -
logA)) + m * digamma(b1) - sum(digamma(n1 - t + b1))-
m * digamma(a1 + b1) + m * digamma(a1 + n1 + b1)
print("dn")
print(dn)
print("da")
print(da)
print("db")
print(db)
return(c(dn,upperbound))
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
