我正在尝试对我的用户自定义发行版做一个蒙特卡罗模拟研究,我修改了我的R
代码,我在一篇文章中看到的另一个。
library(rootSolve)
library(Matrix)
library(bbmle)
alpha=4
beta=0.3
rho=2
lambda=0.5
theta=0.2
samp=1000
par1=c(alpha,beta,rho,lambda,theta)
####Define BXIIWG quantile
BXIIWG_quantile=function(alpha,beta,rho,lambda,theta,u){
f=function(x){
beta*x^alpha+lambda*log(1+x^rho)+log(1-u)
}
x=uniroot(f,c(0,100),tol=0.0001)$root
return(x)
}
####Define BXIIWG log-likelihood
BXIIWG_neglogl=function(alpha,beta,rho,lambda,theta){
-sum(log(alpha*beta*x^(alpha-1)*(1+x^rho)+rho*lambda*x^(rho-1)))+(lambda+1)*sum(log(1+x^rho))-n*log(1-theta)+beta*sum(x^alpha)+2*sum(log(1-theta*(1+x^rho)^(-lambda)*exp(-beta*x^alpha)))
}
####Define simulation process of BXIIWG
BXIIWG_simulation=function(size=c(25,50,100,200,400,800),samp,par1){
Mean=vector()
RMSE=vector()
Bias=vector()
for (iter_size in 1:length(size)){
coef1=matrix(NA,samp,5)
colnames(coef1)=c('alpha','beta','rho','lambda','theta')
for (nsamp in 1:samp){
tryCatch(
{
x1_BXIIWG=NULL
q=runif(size[iter_size],0,1)
x1=sapply(q,BXIIWG_quantile,
alpha=par1[1],beta=par1[2],rho=par1[3],lambda=par1[4],theta=par1[5])
###BXIIWG for x1
x1_BXIIWG<-mle2(BXIIWG_neglogl,
start=list(alpha=par1[1],beta=par1[2],rho=par1[3],lambda=par1[4],theta=par1[5]),
method="L-BFGS-B",data=list(x=x1),
lower=c(alpha=0,beta=0,rho=0,lambda=0,theta=0),
upper=c(alpha=Inf,beta=Inf,rho=Inf,lambda=Inf,theta=1),use.ginv=TRUE)
coef[nsamp,]=coef(x1_BXIIWG)
},error=function(e){}
)
}
Mean[length(size)*(0:4)+iter_size]=apply(coef1,2,mean,na.rm=TRUE)
RMSE[length(size)*(0:4)+iter_size]=apply((coef1-matrix(rep(par1,nsamp),
ncol=5,byrow=T))^2,2,function(x){sqrt(mean(x,na.rm=TRUE))})
}
Bias-as.vector(sapply(1:5,function(x){Bias[(length(size)*(x-1)+1):
(length(size)*x)]=Mean[(length(size)*(x-1)+1):(length(size)*x)]-par1[x]}))
samplesize=as.vector(t(mapply(rep,size,5)))
return(cbind(samplesize,Mean,RMSE,Bias))
}
BXIIWGsim1<-BXIIWG_simulation(25,1000,par1)
Mean1<-BXIIWGsim1$Mean
RMSE1<-BXIIWGsim1$RMSE
Bias1<-BXIIWGsim1$Bias
当我运行这些代码时,我获得了以下错误;
(1)In log(1 - theta) : NaNs produced
(2)In log(1 - theta * (1 + x^rho)^(-lambda) * exp(-beta * ... : NaNs produced
以及
(3)No values of Mean,RMSE and Bias are displayed in console
拜托,我需要你帮我解决这些问题。谢谢。
2条答案
按热度按时间c3frrgcw1#
log(x)
的Hi由于负值而产生NaN。假设:
x = c(1,2,3,-4,-5)
log(x)
:生产的NaN
1.需要重新检查数据集。
2.如果有几个否定的条目,那么就试着删除观察结果[再次视情况而定]。
3.使用
abs()
:abs(x)
:答复1,2,3,4,5
.a64a0gku2#
还有一个更简单的解决方案,尝试使用
log1p()
(它计算 log(1+x))。