我正试图对我的用户定义分布进行蒙特卡罗模拟研究,我已经从一篇文章中修改了我的R
代码。以下是我的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) 在log(1 - theta)部分:产生NaNs
(2) 在log(1 - theta * (1 + x^rho)^(-lambda) * exp(-beta * … )部分:产生NaNs
还有:
(3) 控制台没有显示均值、RMSE和偏差的任何值。
请帮我解决这些问题,谢谢。
x
小于零时,log(x)
将会产生NaN
。例如,如果theta > 1
,则1 - theta < 0
,会抛出一个错误。因此,您需要考虑参数值的限制。 - duckmayr