Hello, Everyone!
I used stan for simulation research and encountered a problem. The error is as follows:
define BOOST_NO_CXX11_RVALUE_REFERENCES
^
:0:0: note: this is the location of the previous definition
Error in FUN(X[[i]], …) : Stan does not support NA (in x) in data
My program is as follows：
timestart <- Sys.time()
options(digits=22)
library(rstan)
pnum<-100
inum<-20
var<-c(4,16,64)
iter=2
inmc.vc<-array(0,dim=c(iter,length(var)))
inmc.se<-array(0,dim=c(iter,length(var)))
vcbias<-array(0,dim=c(iter,length(var)))
vcrmse<-array(0,dim=c(iter,length(var)))
inmc.ci50<-array(0,dim=c(iter,length(var)))
inmc.ci950<-array(0,dim=c(iter,length(var)))
inmc.cover<-array(0,dim=c(iter,length(var)))
for (k in 1:iter){
pz<-rnorm(pnum,0,(var[1])^0.5)
iz<-rnorm(inum,0,(var[2])^0.5)
piz<-array(rnorm((pnum*inum),0,(var[3])^0.5),dim=c(pnum,inum))
x<-array(0,dim=c(pnum,inum))
for(p in 1:pnum){
for(i in 1:inum){
x[p,i]<-0+pz[p]+iz[i]+piz[p,i]
}
}
file_path<-“C:/Users/Administrator/Desktop/ms/data/data_”
write.csv(x,file=paste(file_path,k,".csv",sep = “”),row.names = F )
library(simFrame)
a<-read.csv(file=paste(file_path,k,".csv",sep = “”),header = TRUE )
n<-1
for(n in 1:1000){
set.seed(n)
nc<-NAControl(NArate=0.05) ##MCAR missing rate
x<-setNA(a,nc) ###insert into missing data from x
file_path<-“C:/Users/Administrator/Desktop/ms/datamiss/datamiss_”
write.csv(x,file=paste(file_path,k,".csv",sep = “”),row.names = F )
}
gtdata<-list(np=pnum,ni=inum,vc=var,x=x)
mcmc.model<-‘data{
int<lower=1> np;
int<lower=1> ni;
real x[np,ni];
real vc[3];
}
parameters{
real<lower=0> sigmap;
real<lower=0> sigmai;
real<lower=0> sigmae;
real u;
vector[np] P;
vector[ni] I;
}
transformed parameters{
real sigmap2;
real sigmai2;
real sigmae2;
sigmap2=sigmap^2;
sigmai2=sigmai^2;
sigmae2=sigmae^2;
}
model{
//define variable
real mu;
//likelihood
for (p in 1:np){
P[p]~normal(0,sigmap);
}
for (i in 1:ni){
I[i]~normal(0,sigmai);
}
for (p in 1:np){
for (i in 1:ni){
mu=u+P[p]+I[i];
x[p,i]~normal(mu,sigmae);
}
}
//prior
u~normal(0,sqrt(1000));
target+=inv_gamma_lpdf(sigmap2|2,4);
target+=inv_gamma_lpdf(sigmai2|2,16);
target+=inv_gamma_lpdf(sigmae2|2,64);
}’
require(rstan)
require(reshape)
mcmc.fit<-stan(model_code=mcmc.model,data=gtdata,chains=2,iter=1000,warmup=500,algorithm=“NUTS”, pars=c(“sigmap2”,“sigmai2”,“sigmae2”))
result<-monitor(mcmc.fit,probs=c(0.05,0.95),print=F)
inmc.vc[k,]<-result[1:3,1]
inmc.se[k,]<-result[1:3,3]
inmc.ci50[k,]<-result[1:3,4]
inmc.ci950[k,]<-result[1:3,5]
inmc.cover[k,]<-ifelse(inmc.ci50[k,]<var&var<inmc.ci950[k,],1,0)
vcrmse[k,1]<-(inmc.vc[k,1]-var[1])^2
vcrmse[k,2]<-(inmc.vc[k,2]-var[2])^2
vcrmse[k,3]<-(inmc.vc[k,3]-var[3])^2
}
timeend <- Sys.time()
runningtime <- timeend-timestart
sink(“C:/Users/Administrator/Desktop/ms/output.txt”,split=T)
print(“vc: p,i,pi”)
colMeans(inmc.vc)
print(“se: p,i,pi”)
colMeans(inmc.se)
print(“vc bias: p,i,pi”)
vcbias<-(colMeans(inmc.vc)-var)
print(vcbias)
print(“vc rmse: p,i,pi”)
print(sqrt((colSums(vcrmse))/iter))
print(“Credible Interval: 5%~95%,cover”)
colMeans(inmc.ci50)
colMeans(inmc.ci950)
colMeans(inmc.cover)
print(runningtime)
sink()
I want everyone to help solve it.
How should I modify the statement and let the missing data be estimated as an unknown parameter of the model?
My research is to estimate the results of missing data.