llb_321 发表于 2021-3-4 20:09
换句话说,你的函数找不到最优参数
谢谢你的回复,确实是我简化的函数设定的问题,我把简化函数改了,但还是会有问题,可以麻烦你再看一下吗?我的目的是估计“para”的值,使U_Cur/(U_Cur+U_Tgt)最大,其中U_Cur和U_Tgt是exp(A)/(exp(A)+exp(B))的形式。这是我重新设定的函数:
- library("maxLik")
- integrand <− function(x,para)
- {
- # print(para)
-
- # dv, dv_T, dn, dn_T are the actual input data, here only gives one set of value for easy debugging
- dv <- -1.4527280
- dv_T <- 0.2038909
- dn <- 17.5174383
- dn_T <- 26.6993828
-
- beta10 = para[1];
- beta11 = para[2];
- beta12 = para[3];
- alpha1 = para[4];
- beta20 = para[5];
- beta21 = para[6];
- beta22 = para[7];
- alpha2 = para[8];
-
- U_Cur = exp(beta10 + beta11*dv + beta12*dn + alpha1*x);
- U_Tgt = exp(beta20 + beta21*dv_T + beta22*dn_T + alpha2*x);
-
- # print(U_Cur)
- # print(U_Tgt)
-
- P <- (U_Cur/(U_Cur+U_Tgt))
-
- # print(P)
-
- return(P)
- }
- integrand2 <− function(para1) {integrate( integrand, lower = -2 , upper = 2, para=para1)$value }
- MLE<-maxLik(logLik=integrand2,start=c(0.5, 0.2, 0.1, 0.3, -0.5, 0.2, 0.1, 0.6))
复制代码我通过调试,发现报错的原因是U_Cur或U_Tgt有时可以是inf或0。我把U_Cur/(U_Cur+U_Tgt)加个log再变换,或者变成把U_Cur=exp(A)/(exp(A)+exp(B))变成1-(exp(B)/(exp(A)+exp(B)))等方法都试了。还是会报错。请问还有别的什么办法吗?