# Read data
y <- scan("schiz.asc")
y <- matrix(y,nrow=17,byrow=T)
y <- log(y)
n <- nrow(y); # number of individuals
n1 <- 11 # number of nonschizophrenics
n2 <- 6  # number of schizophrenics
m <- ncol(y) # number of obs. per individual
# Observed indicator variable
S.j <- NULL
for (i in 1:nrow(y)) if (i<12) S.j[i]=0 else S.j[i] = 1 
# Crude estimates
alpha.0 <- rowMeans(y)
sigma2.y.0 <- mean(rowVars(y[1:11,]))
mu.0 <- mean(alpha.0*(1-S.j))
beta.0 <- mu.0-mean(alpha.0*S.j)
sigma2.alpha.0 <- var(alpha.0)
lambda.0 <- 1/3
tau.0 <- 1
# Starting points
mu <- mu.0/rchisq(100,1)
sigma2.y <- sigma2.y.0*rchisq(100,4)
beta <- beta.0/rchisq(100,1)
lambda <- lambda.0/rchisq(100,1)
tau <- tau.0/rchisq(100,1)

alpha <- matrix(NA,100,nrow(y))
for (i in 1:100) alpha[i,] <- alpha.0/rchisq(17,1)

# ECM 
log.post <- function() {
         first <- matrix(NA,n,1)
	 second <- matrix(NA,n,m)
         for (i in 1:n){ 
             first[i] <- dnorm(alpha.s[i],mu.s+beta.s*S.j[i],sqrt(s2.a))

           for (j in 1:m) {
	   second[i,j] <- (1-l.s*S.j[i])*dnorm(y[i,j],alpha.s[i],sqrt(s2.a)) +
	     l.s*S.j[i]*dnorm(y[i,j],alpha.s[i]+tau.s,sqrt(s2.y))  } }

	 third <- prod(first)*prod(second)
	 return(third) }
	   

z.ij <- function(i,j,l,a,t,s){ 

     pp <- 1/(1+(1-l)/l*exp(-.5/s*(2*(y[i,j]-a)*t-t^2)))
     return(pp) }


# Gibbs
n.chains <- 2
n.iterations <- 1500
# Storage
mu.s <- beta.s <- s2.a <- l.s <- tau.s <- s2.y <- matrix(NA,n.iterations+1,n.chains)
alpha.s <- array(NA,c(n.iterations+1,n,n.chains))
  # Initialize chains
    mu.s[1,] <- c(0,0)
    s2.y[1,] <- c(0.18^2,0.19^2)
    beta.s[1,] <- c(0.32,0.30)
    l.s[1,] <- c(.10,.12)
    tau.s[1,] <- c(0.82,0.84)
    s2.a[1,] <- c(0.14^2,0.16^2)
    for (i in 1:n.chains) alpha.s[1,,i] <- al+(i-1)/2

 
    # Update alpha.j

      for (k in 1:n) {
       a.var <- 1/s2.a[i-1,j] + 30/s2.y[i-1,j]
       a.mean <- ((mu.s[i-1,j]+beta.s[i-1,j]*S.j[k])/s2.a[i-1,j] +
               sum(y[k,]-xi[k,]*tau.s[i-1,j])/s2.y[i-1,j])/a.var
       alpha.s[i,k,j] <- rnorm(1,a.mean,sqrt(a.var)) 
                                                    }
    # Update lambda

      h <- sum(xi[12:n,])
      l.s[i,j] <- rbeta(1,h+1,180-h+1)

    # Update sigma2.y

      nonc <- sum((y-alpha.s[i,,j]%*%matrix(1,1,m)-xi*tau.s[i-1,j])^2)
      s2.y[i,j] <- nonc/rchisq(1,508)

    # Update sigma2.a

      s2.a[i,j] <- sum((alpha.s[i,,j]-mu.s[i-1,j]-beta.s[i-1,j]*S.j)^2)/rchisq(1,15)

    # Update tau

      tau.mean <- sum((y-alpha.s[i,,j]%*%matrix(1,1,m))*xi)/sum(xi)
      tau.var <- s2.y[i,j]/sum(xi)
      tau.s[i,j] <- abs(rnorm(1,tau.mean,sqrt(tau.var)))

    # Update mu

      mu.s[i,j] <- rnorm(1,sum(alpha.s[i,,j]-beta.s[i-1,j]*S.j)/17,
		      sqrt(s2.a[i,j]/17))

    # Update beta

      beta.s[i,j] <- rnorm(1,sum(alpha.s[i,12:17,j]-mu.s[i])/6,sqrt(s2.a[i,j]/6))

    } # go to next iteration

  } # go to next chain


# Calculate Potential scale reduction

names1 <- c("a1","a2","a3","a4","a5","a6","a7","a8","a9","a10","a11","a12",
	"a13","a14","a15","a16","a17")

names2 <- c("mu","beta","tau","lambda","s.y","s.a")


output <- array(NA,c(n.iterations,n.chains,length(names2)))
output[,,1] <- mu.s[2:(n.iterations+1),]
output[,,2] <- beta.s[2:(n.iterations+1),]
output[,,3] <- tau.s[2:(n.iterations+1),]
output[,,4] <- l.s[2:(n.iterations+1),]
output[,,5] <- sqrt(s2.y[2:(n.iterations+1),])
output[,,6] <- sqrt(s2.a[2:(n.iterations+1),])

for (i in 1:length(names2)){
    I <- gandr.conv(output[,,i])
    print(c(names2[i],round(I$quantiles,2),round(I$confshrink,3)))  }