MCMC linear model template

Generating simulated data with known parameters to recapture

nReps<-10
x<-seq(from=0,to=100,by=1)
x<-rep(x,nReps)
b0<-1.0; b1<-0.2
ymean<-b0+b1*x
y<-rnorm(n=length(ymean),mean=ymean,sd=1.5)
N<-length(y)

wtData<-data.frame(weight_lbs=y, height_in=x)

Metropolis Algorithym

1. Intialize the chain at some initial value of p

2. Generate a proposed value of p

3. Calculate acceptabce probability r

4. Accept proposed value of p with probability r

5. Repeat n times

nIter<-...
proposalSD<-...
priorMean<-...; priorSD<-...; alphaPrior<-...; betaPrior<-...
b0<-b1<-sd<-vector()
b0[1]<-...; b1[1]<-...; sd[1]<-... # Initial values



for(i in 1:nIter){ # Iterations in sample: i loop
  

    # first sample b0 using a Metropolis step
    b0Star<-rnorm(n=1,...)
    r<-min(1, exp( posteriorLik(...) - posteriorLik(...) ) )
    if(runif(1)<r) {
      b0[i+1]<-b0Star
    } else {
      b0[i+1]<-b0[i]
    } 
    
    # then sample b1 using a Metropolis step
    


   # then sample sd^2 using a Gibbs step or a Metropolis step
   wtPred<-...
   myVar<-sampNormVar(...)
   sd[i+1]<-sqrt(myVar)
  
} 

Examining results using Coda

library(coda)
burnIn<-...
dataMCMCburnIn<-as.mcmc(cbind(b0[-c(1:burnIn)],b1[-c(1:burnIn)],sd[-c(1:burnIn)]))
densplot(dataMCMCburnIn)
HPDinterval(dataMCMCburnIn)
autocorr.plot(dataMCMCburnIn)
plot(dataMCMCburnIn)

Function defintions

sampNormVar<-function(yPred,alpha,beta,data){
    require(invgamma)
    yObs<-data$weight_lbs
    n<-length(yPred)
    alphaPrime<-...
    betaPrime<-...
    rinvgamma(n=1,shape=..., rate=...)
  }

  

Data Likelihood

dataLik<-function(b0,b1,sd,data){
    wt<-...
    ht<-...
    wtPred<-b0+b1*...
    llik<- sum(dnorm(...,log=TRUE))
    return(llik)
  }

Normal Prior on B0 and B1

normPriorL<-function(beta,priorMean,priorSD){
    lik<-dnorm(...,log=TRUE)
    return(...)
  }

Inverse Gamma prior

invGamma<-function(sd,alphaPrior,betaPrior){
    require(invgamma)
    var<-...
    lik<-dinvgamma(var...)
    return(log(...))
  }

Posterior Likelihood

posteriorLik<-function(...){

llik<-dataLik(...) + normPriorB0(...) + normPriorB1(...) + invGammaL(...)

return(llik)
}