
library(mvtnorm)
source('/Users/zzhang4/zzy/Classes/Applied\ Bayesian\ Methods/Spring2010/rtobugs.r')
grm.linear<-function(N, T, R){	
## Constants and parameters
mL<-60
vL<-200
mS<-3
vS<-3
vLS<-0
vE<-100

mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))

y<-array(NA, dim=c(N, T))
LS<-array(NA, dim=c(N, 2))
p<-array(NA, dim=c(N, T))

for (i in 1:N){
  LS[i,]<-rmvnorm(1, mu, sigma)
  for (j in 1:T){
    y[i, j] <- LS[i,1] + LS[i,2]*(j-1) + rnorm(1, 0, sqrt(vE))
    p[i, j] <- pnorm(-2 + LS[i,1]*.02 + LS[i,2]*.02)
  }
}

## generate missing data here
m<-array(runif(N*T), dim=c(N,T))
M<- (m