READ THIS FIRST:

This file contains *five files* of R code that implement the slice sampler for the simplex parameterization for the smoothed ANOVA example in section 4.1 of He, Hodges, & Carlin, "Re-considering the variance parameterization in multiple precision models", Bayesian Analysis, tentatively accepted.  

Before attempting to run this code, you need to split this one file back into the original five files.  The beginning of each file is denoted by 

XXXXXXXXXXXXXXXXXXXXX BEGIN [program name] XXXXXXXXXXXXXX.

The end of each file is denoted by

XXXXXXXXXXXXXXXXXXXXX END [program name] XXXXXXXXXXXXXX.

The main file, listed first below, calls the other four files, which contain code to prepare the dataset ("data.prepare.r"), miscellaneous R functions used in the main program ("SS.splx.fun.r"), code for computing effective sample size ("ess.r"), and for handling output ("SS.splx.output.r").  These latter four files are included in that order after the main file.  The file "data.prepare.r" includes the dataset itself, a 2^3 design with 6 replicates per cell, originally given in Hodges & Sargent (2001 Biometrika, 88:367-379).  


XXXXXXXXXXXXXXXXXXXXXXXXXXX  BEGIN Main R program XXXXXXXXXXXXXXXXXXXXXXXX

source("data.prepare.r")
source("SS.splx.fun.r")
source("ESS.r")

ar<-br<-0.01

#Slice Sampler for BH reparameterization
#"gamma" in this program is "lambda" in He, Hodges, Carlin
# starting values
beta0<-rep(1/(nbeta+1),nbeta) #sum of beta=1
gamma0<-10
# chain specifications
nchain<-1
niter<-10000
count<-2
update<-100
burn<-5000

#par(mfrow=c(3,1))
 
#ptm <- proc.time() #record system time
  
  #1 INITIALIZE
  beta<-array(NA, dim = c(niter,nbeta,nchain)) 
  #N distinct precisions, nchains, niterations.

  U<-gamma<-matrix(NA,nrow = niter, ncol = nchain)
  beta[1,,]<-beta0; gamma[1,]<-gamma0
  U[1,]<-runif(1,0, exp(loglike(gamma0,beta[1,,1]) ) )

  #Run the Slice Sampler
  for (j in 1:nchain){
    for (i in 2:niter){
     #print(c(j,i))
     if (i == burn) ptm <- proc.time() #record system time

  #1 GAMMA
     gamma[i,j]<-SS.gamma(U[i-1,j],beta[i-1,,j])

  #2 BETA
     #beta has to be blockwise updated, because sum(beta.c)=1.
     beta[i,,j]<-SS.beta(U[i-1,j],gamma[i,j]) 
     
  #3 U
     U[i,j]<-runif(1,0, exp(loglike(gamma[i,j],beta[i,,j]) ) )
    }
  }
  elapse<- proc.time() - ptm

##########
# OUTPUT #
##########

source("SS.splx.output.r")


XXXXXXXXXXXXXXXXXXXXXXXXXXX  END Main R program XXXXXXXXXXXXXXXXXXXXXXXX




XXXXXXXXXXXXXXXXXXXXXXXXXXX  BEGIN data.prepare.r XXXXXXXXXXXXXXXXXXXXXXXX

#######################################################
# Data Preparation for run soft materials data set in #
# Hodges & Sargent (2001 Biometrika, Section 6)       #
#######################################################
m1p1c1<-c(4.19,3.99,3.39,3.86,3.44,4.04)
m1p1c0<-c(3.79,3.74,3.30,3.81,3.92,3.86)

m1p0c1<-c(4.19,4.70,3.49,4.60,4.15,4.41)
m1p0c0<-c(4.32,4.40,4.50,4.35,4.91,4.22)

m0p1c1<-c(3.32,3.26,3.72,3.64,3.73,3.38)
m0p1c0<-c(3.33,3.40,3.30,3.15,3.48,3.06)

m0p0c1<-c(3.20,3.09,3.41,3.65,3.04,3.42)
m0p0c0<-c(2.88,2.95,4.82,3.36,3.88,2.79)

y<-c(m1p1c1,m1p1c0,m1p0c1,m1p0c0,m0p1c1,m0p1c0,m0p0c1,m0p0c0)

############################################################
# This data is 2^3 format with 6 replications, so c=8, n=6
# 3 2-way interactions, 1 3-way interaction
############################################################
c=8;n=6
cn<-c*n

################################
#create orthogonal design matrix
################################
one<-rep(1,8)
mould<-c(rep(1,4),rep(-1,4))
pigment<-rep(c(1,1,-1,-1),2)
cure<-rep(c(1,-1),4)

mp<-mould*pigment
mc<-mould*cure
pc<-pigment*cure
mpc<-mould*pigment*cure

H<-cbind(one,mould,pigment,cure,mp,mc,pc,mpc)
X<-kronecker(H,rep(1,6))

##############################
# FIND A1, A2, Py
##############################

#to make A1'%*%A1=cn*I_M; A2'%*%A2=cn*I_N
A1<-X[,1:4]
M<-4
A2<-X[,5:8]
N<-4
# there are only 3 free beta parameters, the last 
# beta = 1-beta1-beta2-beta3nbeta<-3 

yStar<-t(A2)%*%y
#when there is no replication, Py=0, else Py!=0.
Py<-ssq(y)-ssq(t(A1)%*%y)/cn-ssq(t(A2)%*%y)/cn 

ys2<-yStar^2
nj<-rep(1,4)

###############################
# Prior Information
###############################

aEta<-bEta<-0.01 #G(.01,.01) prior for eta
#the prior is induced by the uniform on degree of freedoms rho1,...,rho4.
f.prior<-function(gamma,beta) {
   beta.s<-1-sum(beta)
   beta.c<-c(beta,beta.s) #complete list of beta1-beta4
   A<-matrix(0,4,4)
   for (i in 1:4) {A[i,1]<-beta.c[i]/48/(1+gamma*beta.c[i]/48)^2}
   A[1,2]<-gamma/48/(1+gamma*beta[1]/48)^2
   A[2,3]<-gamma/48/(1+gamma*beta[2]/48)^2
   A[3,4]<-gamma/48/(1+gamma*beta[3]/48)^2
   A[4,2:4]<- -gamma/48/(1+gamma*beta.c[4]/48)^2
   return(abs(det(A)))
}

# This prior is induced, so it is ugly. This prior is used to compare 
# the results to Hodges & Sargent's paper.


XXXXXXXXXXXXXXXXXXXXX END data.prepare.r XXXXXXXXXXXXXX.



XXXXXXXXXXXXXXXXXXXXX BEGIN SS.splx.fun.r XXXXXXXXXXXXXX.

###########
#2. ssq  ##
###########
ssq<-function(x) {
   sum(x^2)
}

#####################
#3. inner product  ##
#####################

inner<-function(x,y) {
   t(x)%*%y
}

#################
#4. vector norm #
#################

norm<-function(x) {
   sqrt(ssq(x))
}

#####################
#5. log likelihood  #
#####################
# the argument "adj" adjusts the log likelihood up or down
# to avoid numerical problems.

loglike<-function(gamma,beta,adj=50){
   beta.s<-1-sum(beta)
   beta.c<-c(beta,beta.s) #complete list of beta
   G<- log(1+cn/(gamma*beta.c))
   b<-aEta+(cn-M)/2
   R<-bEta+.5*ssq(y)-.5*ssq(t(A1)%*%y)/cn-.5*sum(ys2/(cn+gamma*beta.c))

   prior <- dgamma(gamma, ar, rate=br, log=T)
   logl <- -0.5*sum(G)-b*log(R)
   return(prior+logl+adj)
}

#################
#6. draw gamma  #
#################
SS.gamma<-function(u,beta) {   
   repeat{
      try<-runif(1,0,2000)
      logL<-loglike(try,beta)
      if (log(u)<logL) break
   }
   return (try)
}

#################
#7. draw beta   #
#################
SS.beta<-function(u,gamma) {
   
   beta<-rep(NA,nbeta)
   repeat{   
      beta<-randBeta()      
      logL<-loglike(gamma,beta)
      if (log(u)<logL) break
   }
   return (beta)
}

#########################
#8. generate r.v. beta  #
#########################

randBeta<-function() {
   #generate beta s.t. sum(beta)=1
   exp<-rexp(nbeta+1)
   S.exp<-sum(exp)
 
   #each individual beta is uniformly distributed given S.exp
   beta<-exp[1:nbeta]/S.exp 
   return(beta) 
}


XXXXXXXXXXXXXXXXXXXXX END SS.splx.fun.r XXXXXXXXXXXXXX.




XXXXXXXXXXXXXXXXXXXXX BEGIN ess.r XXXXXXXXXXXXXX.

##########################################################################
#ESS related functions
#Reference: (1) MCMC in practice, page 50. Chapman & Hall / CRC press. 
#               Richardson & Spielhalter
#           (2) Geyer 1992. Practical MCMC, Stat. Science 7:473-483.
##########################################################################

# Gamma(m) = gamma(2*m) + gamma(2*m+1), where gamma(.) is autocorrelation
# with lag (.)
Gm = function(x, m) {
     n = length(x)
     if (m==0) {g2m = var(x)*(n-1)/n} else {
        g2m = acf(x, plot=F, lag.max=2*m, type="covariance")$acf[2*m+1]
     }
     g2m1 = acf(x, plot=F, lag.max=2*m+1, type="covariance")$acf[2*m+2]
     ans = g2m + g2m1
     return(ans)
}

#find the largest m* such that Gamma(m*) is positive. 
f.m.pos = function(x) {
  m <- 0
  repeat {
     G = Gm(x,m)
     if (G < 1e-6) {m.pos = m; break}
     m = m + 1
  }
  return(m.pos)
}

#find the largest m* such that Gamma(m*) is positive and monotonely decreasing. 
f.m.mon = function(x) {
  m <- 0
  G1 = Gm(x,m)
  repeat {
     G2 = Gm(x,m+1)

     if (G2 < 1e-6 || G2 >= G1) {m.mon = m; break;}
     m = m + 1
     G1 = G2
  }
  return(m.mon)
}

#find the largest m* such that Gamma(m*) is positive, monotonely 
# decreasing and convex. 
f.m.con = function(x) {
  m <- 0
  G1 = Gm(x,m); G2 = Gm(x,m+1)
  repeat {     
     G3 = Gm(x, m+2)
     is.convex = (G2 < 0.5*(G1+G3))
     if (G2 < 1e-6 || G2 >= G1 || !is.convex) {m.con = m; break;}
     m = m + 1
     G1 = G2
     G2 = G3
  }
  return(m.con)
}
#In general, m.convex < m.mon < m.positive. 

ESS<-function(x){
  n = length(x)
  #effective sample size for parameter x

  #m.pos = f.m.pos(x)
  #m.mon = f.m.mon(x)
  m.con = f.m.con(x)
  ess.temp = function(m) {
    if (m==0) {a = 0} else { 
       a = acf(x,plot=F,lag.max=2*m)$acf[2:(2*m+1)]
    }
    kappa <- 1 + 2*sum(a)
    ESS = n / kappa
    return(ESS)
  }
  #ESS.pos = round(ess.temp(m.pos))
  #ESS.mon = round(ess.temp(m.mon))
  #ESS.con = round(ess.temp(m.con))

  ESS = round(ess.temp(m.con))

  #ESS = list()
  #ESS$ESS.positive = ESS.pos; 
  #ESS$ESS.monotone = ESS.mon; 
  #ESS$ESS.convex = ESS.con

  return(ESS)
}


XXXXXXXXXXXXXXXXXXXXX END ess.r XXXXXXXXXXXXXX.



XXXXXXXXXXXXXXXXXXXXX BEGIN SS.splx.output.r XXXXXXXXXXXXXX.

dput(beta,file="beta.dat")
dput(elapse,file="elapse.dat")
dput(gamma,file="gamma.dat")
dput(U,file="U.dat")

select = (burn+1):niter
out<-list()

out$ESS.g<-ESS(gamma[select,1])

out$ESS1<-ESS(beta[select,1,1])
out$ESS2<-ESS(beta[select,2,1])
out$ESS3<-ESS(beta[select,3,1])

out$ESSpersec.g<-out$ESS.g/elapse[1]
out$ESSpersec1<-out$ESS1/elapse[1]
out$ESSpersec2<-out$ESS2/elapse[1]
out$ESSpersec3<-out$ESS3/elapse[1]

out$acf.g<-acf(gamma[select,1],plot=F)[c(1,5,10)]
out$acf.1<-acf(beta[select,1,1],plot=F)[c(1,5,10)]
out$acf.2<-acf(beta[select,2,1],plot=F)[c(1,5,10)]
out$acf.3<-acf(beta[select,3,1],plot=F)[c(1,5,10)]

dput(out, "ESS.dat")


XXXXXXXXXXXXXXXXXXXXX END ess.r XXXXXXXXXXXXXX.