Вопросы о R2Winbugs (ожидаемый оператор сбора c) - PullRequest
Новая
       91

Вопросы о R2Winbugs (ожидаемый оператор сбора c)

0 голосов
/ 13 апреля 2020

Мой код синтаксически правильный. Тем не менее, я сталкиваюсь с неправильным сообщением в заголовке, подробности в следующем: 

## MAR data analysis
library(MASS) #Load the MASS package
library(R2WinBUGS) #Load the R2WinBUGS package
library(boa) #Load the boa package
load('F:/Final_project/Y_osn.Rdata')
load('F:/Final_project/index_1_Y_osn.Rdata')
load('F:/Final_project/index_2_Y_osn.Rdata')
Y1<- data.matrix(Y_osn[index_1_Y_osn,]);
Y2<-data.matrix(Y_osn[index_2_Y_osn,]);
N1<-length(Y1[,1]); 
N2<-length(Y2[,1]); 
PP<-length(Y_osn[1,]);
TT<-2# spilt the data into 2 parts: one for setting initial value,
#phi<-matrix(data=c(1.0,0.3,0.3,1.0),ncol=3) #The covariance matrix of xi
Ro<-matrix(data=c(1.0,0.3,0.3,0.3,1,0.3,0.3,0.3,1.0), ncol=3)
#yo<-matrix(data=NA,nrow=N,ncol=P); 
NY=PP # dimension of Y
Neta=1 # dimension of eta
Nxi=3 # dimension of xi
Ngam=3 # dimension of gamma


#Matrices save the Bayesian Estimates and Standard Errors
Eu<-matrix(data=NA,nrow=TT,ncol=NY)
SEu<-matrix(data=NA,nrow=TT,ncol=NY)
Elam<-matrix(data=NA,nrow=TT,ncol=NY-Neta-Nxi)
SElam<-matrix(data=NA,nrow=TT,ncol=NY-Neta-Nxi)
Egam<-matrix(data=NA,nrow=TT,ncol=Ngam)
SEgam<-matrix(data=NA,nrow=TT,ncol=Ngam)
Ephx<-matrix(data=NA,nrow=TT,ncol=6)
SEphx<-matrix(data=NA,nrow=TT,ncol=6)
Esgd<-numeric(TT); 
SEsgd<-numeric(TT);
Esgm<-matrix(data=NA,nrow=TT,ncol=NY-12)
SEsgm=matrix(NA,nrow=TT,ncol=NY-12)

#Arrays save the HPD intervals
uby=array(NA, c(TT,9,2))
#ubeta=array(NA, c(TT,2))
lam=array(NA, c(TT,6,2))
gam=array(NA, c(TT,3,2))
sgd=array(NA, c(TT,2))
phx=array(NA, c(TT,3,2))
#co=matrix(data=NA,nrow=N,ncol=1)




## calculate the preassigned interval
index_m1<-c(2,3,11,22);
index_m2<-c(2,10,21,32);
num_index<-c(2,4,2,3);
count_f1<-matrix(NA,nrow=index_m2[1]-index_m1[1]+1,ncol=num_index[1]);#trusted 2 category eta[2]
count_f2<-matrix(NA,nrow=index_m2[2]-index_m1[2]+1,ncol=num_index[2]);#eta[3,4]& xxi1
count_f3<-matrix(NA,nrow=index_m2[3]-index_m1[3]+1,ncol=num_index[3]);#xxi2
count_f4<-matrix(NA,nrow=index_m2[4]-index_m1[4]+1,ncol=num_index[4]);#xxi3

## calculate the normal quantile
quantile_f1<-matrix(NA,nrow=index_m2[1]-index_m1[1]+1,ncol=num_index[1]);#trusted 2 category eta[2]
quantile_f2<-matrix(NA,nrow=index_m2[2]-index_m1[2]+1,ncol=num_index[2]);#eta[3,4]& xxi1
quantile_f3<-matrix(NA,nrow=index_m2[3]-index_m1[3]+1,ncol=num_index[3]);#xxi2
quantile_f4<-matrix(NA,nrow=index_m2[4]-index_m1[4]+1,ncol=num_index[4]);#xxi3

## calculate the interval 
interval_index_f1<-matrix(NA,nrow=index_m2[1]-index_m1[1]+1,ncol=num_index[1]+1);#trusted 2 category eta[2]
interval_index_f2<-matrix(NA,nrow=index_m2[2]-index_m1[2]+1,ncol=num_index[2]+1);#eta[3,4]& xxi1
interval_index_f3<-matrix(NA,nrow=index_m2[3]-index_m1[3]+1,ncol=num_index[3]+1);#xxi2
interval_index_f4<-matrix(NA,nrow=index_m2[4]-index_m1[4]+1,ncol=num_index[4]+1);#xxi3

## frequency statistic

for (j in index_m1[1]:index_m2[1]){
  count_f1[(j-index_m1[1]+1),]<-as.numeric(table(Y_osn[,j]));
  quantile_f1[(j-index_m1[1]+1),]<-cumsum(count_f1[(j-index_m1[1]+1),])/sum(count_f1[(j-index_m1[1]+1),]);
  interval_index_f1[(j-index_m1[1]+1),1]<--3000;
  interval_index_f1[(j-index_m1[1]+1),num_index[1]+1]<-3000;
  interval_index_f1[(j-index_m1[1]+1),num_index[1]]<-0
}

for (j in index_m1[2]:index_m2[2])
{
  count_f2[(j-index_m1[2]+1),]<-as.numeric(table(Y_osn[,j]));
  quantile_f2[(j-index_m1[2]+1),]<-cumsum(count_f2[(j-index_m1[2]+1),])/sum(count_f2[(j-index_m1[2]+1),]);
  interval_index_f2[(j-index_m1[2]+1),1]<--3000;
  interval_index_f2[(j-index_m1[2]+1),num_index[2]+1]<-3000;
  for (k in 2: num_index[2]){
    interval_index_f2[(j-index_m1[2]+1),k]<-qnorm(quantile_f2[(j-index_m1[2]+1),k-1]);
  }
}


for (j in index_m1[3]:index_m2[3]){
  count_f3[(j-index_m1[3]+1),]<-as.numeric(table(Y_osn[,j]));
  quantile_f3[(j-index_m1[3]+1),]<-cumsum(count_f3[(j-index_m1[3]+1),])/sum(count_f3[(j-index_m1[3]+1),]);
  interval_index_f3[(j-index_m1[3]+1),1]<--3000;
  interval_index_f3[(j-index_m1[3]+1),num_index[3]+1]<-3000;
  interval_index_f3[(j-index_m1[3]+1),num_index[3]]<-0
}

for (j in index_m1[4]:index_m2[4]){
  count_f4[(j-index_m1[4]+1),]<-as.numeric(table(Y_osn[,j]));
  quantile_f4[(j-index_m1[4]+1),]<-cumsum(count_f4[(j-index_m1[4]+1),])/sum(count_f4[(j-index_m1[4]+1),]);
  interval_index_f4[(j-index_m1[4]+1),1]<--3000;
  interval_index_f4[(j-index_m1[4]+1),num_index[4]+1]<-3000;
  for (k in 2: num_index[4]){
    interval_index_f4[(j-index_m1[4]+1),k]<-qnorm(quantile_f4[(j-index_m1[4]+1),k-1])
  }
}

DIC=numeric(TT) #DIC values
#Parameters to be estimated
#parameters<-c("uby","ubeta","lam","gam","sgd","phx")
parameters<-c("uby","lam","gam","sgm","sgd","phx")

#Initial values for the MCMC in WinBUGS

init1=list(uby=rep(0,NY),lam=rep(0,NY-Neta-Nxi),gam=rep(0,Ngam),psi=rep(1,NY-12),psd=1,phi=matrix(c(1.0,0.3,0.3,0.3,1,0.3,0.3,0.3,1.0),nrow=3),xi=matrix(rep(1,N1*3),ncol=3))

init2=list(uby=rep(1,NY),lam=rep(1,NY-Neta-Nxi),gam=rep(1,Ngam),psi=rep(0.5,NY-12),psd=2,phi=matrix(c(1.0,0.1,0.1,0.1,1,0.1,0.1,0.1,1.0),nrow=3),xi=matrix(rep(0,N1*3),ncol=3))

inits<-list(init1, init2)
# Run WINBUGS

  #Input data set for WinBUGS
  data<-list(N=N1,zero2=c(0,0,0),P=PP,R=Ro,z=Y1,thdone=interval_index_f1,thdtwo=interval_index_f2,thdthree=interval_index_f3,thdfour=interval_index_f4)
  #Call WinBUGS
    model<-bugs(data,inits,parameters,model.file="F:/Final_project/ignore/model1.txt",
              n.chains=2,n.iter=200,n.burnin=10,n.thin=1,
              bugs.directory="E:/math/math software/WinBUGS14/",
              working.directory="F:/Final_project/ignore/",debug=TRUE)
  #Save Bayesian Estimates
  Eu[1,]<-model$mean$uby; 
  Elam[1,]<-model$mean$lam;
  Egam[1,]<-model$mean$gam
  Ephx[1,1]<-model$mean$phx[1,1]; 
  Ephx[1,2]<-model$mean$phx[1,2];
  Ephx[1,3]<-model$mean$phx[1,3];
  Ephx[1,4]<-model$mean$phx[2,2]; 
  Ephx[1,5]<-model$mean$phx[2,3];
  Ephx[1,6]<-model$mean$phx[3,3];  
  Esgd[1]<-model$mean$sgd
  Esgm[1]<-model$mean$sgm
  #Save Standard Errors
  SEu[1,]<-model$sd$uby; 
  SElam[1,]<-model$sd$lam;
  SEgam[1,]<-model$sd$gam;
  SEphx[1,1]<-model$sd$phx[1,1]; 
  SEphx[1,2]<-model$sd$phx[1,2];
  SEphx[1,3]<-model$sd$phx[1,3];
  SEphx[1,4]<-model$sd$phx[2,2]; 
  SEphx[1,5]<-model$sd$phx[2,3];
  SEphx[1,6]<-model$sd$phx[3,3]; 
  SEsgd[1,]<-model$sd$sgd
  SEsgm[1,]<-model$sd$sgm
  #Save DIC value
  DIC[1]=model$DIC #DIC

И моя модель: 

model {
for(i in 1:N) {
#Measurement equation model
y[i,1]~dnorm(mu[i,1],psi[1])
y[i,2]~dnorm(mu[i,2],1)I(thdone[z[i,2]],thdone[z[i,2]+1])

for(j in 3:10){
y[i,j]~dnorm(mu[i,j],psi[j-1])I(thdtwo[j-2, z[i,j]],thdtwo[j-2, z[i,j]+1])
}

for (j in 11:21){
y[i,j]~dnorm(mu[i,j],1)I(thdthree[j-10, z[i,j]],thdthree[j-10, z[i,j]+1])
}

for (j in 22:32){
y[i,j]~dnorm(mu[i,j],psi[j-12])I(thdfour[j-21, z[i,j]+1],thdfour[j-21, z[i,j]+2])
}

mu[i,1]<-uby[1]+eta[i]
for (l in 2:4){
mu[i,l]<-uby[l]+lam[l-1]*eta[i]
}

mu[i,5]<-uby[5]+xi[i,1]
for (l in 6:10){
mu[i,l]<-uby[l]+lam[l-2]*xi[i,1]
}

mu[i,11]<-uby[11]+xi[i,2]
for (l in 12:21){
mu[i,l]<-uby[l]+lam[l-3]*xi[i,2]
}

mu[i,22]<-uby[22]+xi[i,3]
for (l in 23:32){
mu[i,l]<-uby[l]+lam[l-4]*xi[i,3]
}

#Structural equation model
xi[i,1:3]~dmnorm(zero2[1:3],phi[1:3,1:3])
eta[i]~dnorm(etamu[i],psd)
etamu[i]<-gam[1]*xi[i,1]+gam[2]*xi[i,2]+gam[3]*xi[i,3]
} #End of i

for(i in 1:3){ zero2[i]<-0 }
#Priors inputs for loadings and coefficients
#setting psi first

for (j in 1:20){
psi[j]~dgamma(9,4)
sgm[j]<-1/psi[j];
}
for (i in 1:28){
lam[i]~dnorm(1,4.0);
}
for (i in 1:P){ 
uby[i]~dnorm(0,4.0) 
}
var.gam<-4.0*psd
gam[1]~dnorm(0,var.gam); 
gam[2]~dnorm(0,var.gam)
gam[3]~dnorm(0,var.gam)
#Priors inputs for precisions
psd~dgamma(9,3); 
sgd<-1/psd;

phi[1:3,1:3]~dwish(R[1:3,1:3], 10)
phx[1:3,1:3]<-inverse(phi[1:3,1:3])
} #end

Итак, что не так с сообщением "Ожидаемая коллекция оператор c ", я не смог выяснить, какая часть не инициализирована. Пожалуйста, помогите мне найти его.

...