On 03/05/2014, 11:49 AM, thanoon younis wrote:> dear all members
> i have error in the code below "Error in Y[i, 9] = 0.9 * XI[i, 2] +
eps[9]
> : subscript out of bounds" is there anyone helps me please.
You created Y with 8 columns, then you refer to column 9.
Duncan Murdoch
>
> many thanks in advance
> thanoon
>
>
> llibrary(mvtnorm) #Load mvtnorm package
> library(R2WinBUGS) #Load R2WinBUGS package
>
> N=500 #Sample size
> BZ=numeric(N) #Fixed covariate in structural equation
> XI=matrix(NA, nrow=N, ncol=2) #Explanatory latent variables
> Eta=numeric(N) #Outcome latent variables
> Y=matrix(NA, nrow=N, ncol=8) #Observed variables
>
> #The covariance matrix of xi
> phi=matrix(c(1, 0.3, 0.3, 1), nrow=2)
>
> #Estimates and standard error estimates
> Eu=matrix(NA, nrow=100, ncol=10); SEu=matrix(NA, nrow=100, ncol=10)
> Elam=matrix(NA, nrow=100, ncol=7); SElam=matrix(NA, nrow=100, ncol=7)
> Eb=numeric(100); SEb=numeric(100)
> Egam=matrix(NA, nrow=100, ncol=5); SEgam=matrix(NA, nrow=100, ncol=5)
> Esgm=matrix(NA, nrow=100, ncol=10); SEsgm=matrix(NA, nrow=100, ncol=10)
> Esgd=numeric(100); SEsgd=numeric(100)
> Ephx=matrix(NA, nrow=100, ncol=3); SEphx=matrix(NA, nrow=100, ncol=3)
>
> R=matrix(c(1.0, 0.3, 0.3, 1.0), nrow=2)
>
> parameters=c("u", "lam", "b",
"gam", "sgm", "sgd", "phx")
>
> init1=list(u=rep(0,10), lam=rep(0,7), b=0, gam=rep(0,5), psi=rep(1,10),
> psd=1, phi=matrix(c(1, 0, 0, 1), nrow=2))
>
> init2=list(u=rep(1,10), lam=rep(1,7), b=1, gam=rep(1,5), psi=rep(2,10),
> psd=2, phi=matrix(c(2, 0, 0, 2), nrow=2))
>
> inits=list(init1, init2)
>
> eps=numeric(10)
>
> for (t in 1:100) {
> #Generate Data
> for (i in 1:N) {
> BZ[i]=rt(1, 5)
>
> XI[i,]=rmvnorm(1, c(0,0), phi)
>
> delta=rnorm(1, 0, sqrt(0.36))
> Eta[i]=0.5*BZ[i]+0.4*XI[i,1]+0.4*XI[i,2]+0.3*XI[i,1]*XI[i,2]
> +0.2*XI[i,1]*XI[i,1]+0.5*XI[i,2]*XI[i,2]+delta
>
> eps[1:3]=rnorm(3, 0, sqrt(0.3))
> eps[4:7]=rnorm(4, 0, sqrt(0.5))
> eps[8:10]=rnorm(3, 0, sqrt(0.4))
> Y[i,1]=Eta[i]+eps[1]
> Y[i,2]=0.9*Eta[i]+eps[2]
> Y[i,3]=0.7*Eta[i]+eps[3]
> Y[i,4]=XI[i,1]+eps[4]
> Y[i,5]=0.9*XI[i,1]+eps[5]
> Y[i,6]=0.7*XI[i,1]+eps[6]
> Y[i,7]=0.5*XI[i,1]+eps[7]
> Y[i,8]=XI[i,2]+eps[8]
> Y[i,9]=0.9*XI[i,2]+eps[9]
> Y[i,10]=0.7*XI[i,2]+eps[10]
> }
>
> #Run WinBUGS
> data=list(N=500, zero=c(0,0), z=BZ, R=R, y=Y)
>
> model<-bugs(data,inits,parameters,
> model.file="C:/Simulation/model.txt",
> n.chains=2,n.iter=10000,n.burnin=4000,n.thin=1,
> bugs.directory="c:/Program Files/WinBUGS14/",
> working.directory="C:/Simulation/")
>
> #Save Estimates
> Eu[t,]=model$mean$u; SEu[t,]=model$sd$u
> Elam[t,]=model$mean$lam; SElam[t,]=model$sd$lam
> Eb[t]=model$mean$b; SEb[t]=model$sd$b
> Egam[t,]=model$mean$gam; SEgam[t,]=model$sd$gam
> Esgm[t,]=model$mean$sgm; SEsgm[t,]=model$sd$sgm
> Esgd[t]=model$mean$sgd; SEsgd[t]=model$sd$sgd
> Ephx[t,1]=model$mean$phx[1,1]; SEphx[t,1]=model$sd$phx[1,1]
> Ephx[t,2]=model$mean$phx[1,2]; SEphx[t,2]=model$sd$phx[1,2]
> Ephx[t,3]=model$mean$phx[2,2]; SEphx[t,3]=model$sd$phx[2,2]
> }
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>