User Tools

Site Tools


lab:projects:02robust_growth_curve_modeling_using_student-t_distribution:bugs_codes_for_data_simulation

Data generation

Model{
  # Model specification for linear growth curve model
  for (i in 1:300){
    LS[i,1:2]~dmnorm(muLS[i,1:2], Inv_cov[1:2,1:2])
    muLS[i,1]<-6
    muLS[i,2]<-1
    for (t in 1:5){
      y[i, t] ~ dnorm(muY[i,t], Inv_Sig_e2[i])
      muY[i,t]<-LS[i,1]+LS[i,2]*(t-1)/4
    }
    w[i]~dgamma(a, a)
    Inv_Sig_e2[i]<-w[i]*Pre_e2
  }

  a<-1.5
  Pre_e2<- 2
}

# The (naive) starting values for model parameters.
list(Inv_cov= structure(.Data = c(.5,0,0,2),.Dim=c(2,2)))

Model{
  # Model specification for linear growth curve model
  for (i in 1:500){
    LS[i,1:2]~dmnorm(muLS[i,1:2], Inv_cov[1:2,1:2])
    muLS[i,1]<-6
    muLS[i,2]<-1
    for (t in 1:5){
      y[i, t] ~ dnorm(muY[i,t], Inv_Sig_e2)
      muY[i,t]<-LS[i,1]+LS[i,2]*(t-1)/4
    }
  }

  Inv_Sig_e2<-2/3

}

# The (naive) starting values for model parameters.
list(Inv_cov= structure(.Data = c(.5,0,0,2),.Dim=c(2,2)))

## model with outliers

Model{
  # Model specification for linear growth curve model
  for (i in 1:N){
    z[i]~dcat(lambda[1:2])
    LS[i,1:2]~dmnorm(muLS[i,1:2], Inv_cov[1:2,1:2])
    muLS[i,1]<-bL[1]
    muLS[i,2]<-bS[1]
    for (t in 1:5){
      y[i, t] ~ dnorm(muY[i,t], Inv_Sig_e2[z[i]])
      muY[i,t]<-LS[i,1]+LS[i,2]*t
    }
  }
  Inv_Sig_e2[1]<- 1/Sig_e2[1]
  Inv_Sig_e2[2]<- 1/Sig_e2[2]
  Inv_cov[1:2,1:2]<-inverse(Cov[1:2,1:2])
}

# The (naive) starting values for model parameters.
list(lambda=c(.95, .05), bL=c(60), bS = c(3), Sig_e2 = c(100, 150), Cov= structure(.Data = c(200,0,0,3),.Dim=c(2,2)), N=500)


## model with outliers

R codes

## Simulated t data
y<-read.table('data-generated-student-coda.txt')

ymat<-array(y[, 2], dim=c(1000,5, 500))

t(ymat[2,,])->temp

plot(1:5, temp[1,1:5], ylim=c(30, 150), type='l')
for (i in 2:200){
	lines(1:5, temp[i,1:5])
	}

## write data into txt files

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('./Tdata/TdataN500I',i,'.txt', sep='')
 write.table(temp, filename, row.names=F, col.names=F)
}


## for normal data

y<-read.table('data-generated-normal-coda.txt')

ymat<-array(y[, 2], dim=c(1000,5, 500))

t(ymat[2,,])->temp

plot(1:5, temp[1,1:5], ylim=c(30, 150), type='l')
for (i in 2:200){
	lines(1:5, temp[i,1:5])
	}

## write data into txt files

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('./Ndata/NdataN500I',i,'.txt', sep='')
 write.table(temp, filename, row.names=F, col.names=F)
}	



## for normal data with outliers

y<-read.table('data-generated-outliers-coda.txt')

ymat<-array(y[, 2], dim=c(1000,5, 500))

t(ymat[3,,])->temp

boxplot(temp[,1],temp[,2],temp[,3],temp[,4],temp[,5])

plot(1:5, temp[1,1:5], ylim=c(30, 150), type='l')
for (i in 2:200){
	lines(1:5, temp[i,1:5])
	}

## write data into txt files

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('./Odata/OdataN500I',i,'.txt', sep='')
 write.table(temp, filename, row.names=F, col.names=F)
}

Page Tools