---
title: "MCMC linear model"
output: html_notebook
---

```{r}
 setwd("~/Documents/Web/Teaching/PBIO_294/Homework/10.25.2017")
```



```{r}
wtData<-read.table("http://www.uvm.edu/~bbeckage/Teaching/DataAnalysis/Data/wts.txt",header=TRUE)
head(wtData)

summary(lm(weight_lbs~height_in,data=wtData))
```

Generating simulated data with known parameters to recapture

```{r}
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)

summary(lm(weight_lbs~height_in,data=wtData))
```


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



```{r}
nIter<-5000
proposalSD<-5
priorMean<-0; priorSD<-100; alphaPrior<-0.001; betaPrior<-0.001 # Prior parameters: priorMean,priorSD, alphaPrior, betaPrior
b0<-b1<-sd<-vector()
b0[1]<-1; b1[1]<-1; sd[1]<-1 # Initial values

# Create likelihood functions
dataLik<-makeDataLik(wtData)
normPriorB0<-makeNormPriorL(priorMean,priorSD)
normPriorB1<-makeNormPriorL(priorMean,priorSD)
invGammaL<-makeInvGamma(alphaPrior,betaPrior)
sampNormSD<-makeSampNormSD(alphaPrior,betaPrior,wtData)

for(i in 1:nIter){ # Iterations in sample: i loop
  
    if(identical(i %% 1000,0)) cat('iter: ',i,fill=TRUE) # keep track of progress
    # first sample b0 using a Metropolis step
    b0Star<-rnorm(n=1,mean=b0[i],sd=proposalSD)
    
    #r<-min(1, ( posteriorLik(b0Star,b1[i],sd[i]) ) / ( posteriorLik(b0[i],b1[i],sd[i]) ) )
    r<-min(1, exp( posteriorLik(b0Star,b1[i],sd[i]) - posteriorLik(b0[i],b1[i],sd[i]) ) )
    if(runif(1)<r) {
      b0[i+1]<-b0Star
    } else {
      b0[i+1]<-b0[i]
    } 
    
    # then sample b1 using a Metropolis step
    
    b1Star<-rnorm(n=1,mean=b1[i],sd=proposalSD)
    
    # r<-min(1, ( posteriorLik(b0[i],b1Star,sd[i]) ) / ( posteriorLik(b0[i],b1[i],sd[i]) ) )
    r<-min(1, exp( posteriorLik(b0[i],b1Star,sd[i]) - posteriorLik(b0[i],b1[i],sd[i]) )  )  
    if(runif(1)<r) {
      b1[i+1]<-b1Star
    } else {
      b1[i+1]<-b1[i]
    } 

   # then sample sd^2 using a Gibbs step
   wtPred<-b0[i]+b1[i]*wtData$height_in
   myVar<-sampNormSD(wtPred)
   sd[i+1]<-sqrt(myVar)
  
} # end of i for loop


```

### Examining results using Coda

```{r}
library(coda)
burnIn<-1000
# dataMCMC<-as.mcmc(cbind(b0,b1,sd))
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 

# Inverse Gamma: Varaince Conjugate Sampling
from Clark, Ecologial Modeling, page 89

```{r}

makeSampNormSD<-function(alpha,beta,data){
  yObs<-data$weight_lbs
  function(yPred){
    require(invgamma)
    # alpha = prior parm; beta = prior parm
    n<-length(yPred)
    alphaPrime<-alpha + 0.5*n
    betaPrime<-beta + 0.5*sum( (yObs-yPred)^2 )
    rinvgamma(n=1,shape=alphaPrime, rate=betaPrime)
  }
}
  
```

### Data Likelihood

```{r}
makeDataLik<-function(data){

  function(b0,b1,sd){
    wt<-data$weight_lbs
    ht<-data$height
    wtPred<-b0+b1*ht
    llik<- sum(dnorm(x=wt,mean=wtPred,sd=sd,log=TRUE))
    #browser()
    # exp(llik)
    return(llik)
  }
  
}

 # dataLik<-makeNormPriorL(wtData)
# ls.str(environment(dataLik))

```

### Normal Prior on B0 and B1

```{r}

makeNormPriorL<-function(priorMean,priorSD){
  function(beta){
    ## Prior for Beta0 and Beta1
    lik<-dnorm(x=beta,mean=priorMean,sd=priorSD,log=FALSE)
    return(log(lik))
  }
}

# normPriorB0<-makeNormPriorL(0,100)
# normPriorB1<-makeNormPriorL(0,100)
# ls.str(environment(normPriorB0))

```

### Inverse Gamma prior

```{r}
makeInvGamma<-function(alphaPrior,betaPrior){
  
  function(sd){
    require(invgamma)
    # alpha = prior parm; beta = prior parm
    var<-sd^2
    lik<-dinvgamma(var,shape=alphaPrior, rate=betaPrior)
    return(log(lik))
  }
}

# invGammaL<-makeInvGamma(0.001,0.001)
# ls.str(environment(invGammaL))

```

### Posterior Likelihood
```{r}

posteriorLik<-function(bo,b1,sd){
# p(data|bo,b1,sd) * p(bo|priorMean,priorSD) * p(b1|priorMean,priorSD) * p(sd|alpha,beta)
# lik<-dataLik(bo,b1,sd) * normPriorB0(b0) * normPriorB1(b1) * invGammaL(sd)
llik<-dataLik(bo,b1,sd) + normPriorB0(b0) + normPriorB1(b1) + invGammaL(sd)
# browser()
return(llik)
}


```


Testing and debugging

```{r}
posteriorLik(b0[i],b1[i],sd[i])
dataLik(b0[i],b1[i],sd[i])
normPriorB0(b0[i])
normPriorB1(b1[i])
invGammaL(sd[i])

dataLik(b0[i],b1Star,sd[i])
normPriorB0(b0[i])
normPriorB1(b1Star)
invGammaL(sd[i])
posteriorLik(b0[i],b1Star,sd[i])

```



### Exploring priors

```{r}
x<-seq(from=0,to=100,length=1000)
y<-dgamma(x,shape=0.001,rate=0.001)
plot(x,y)
```


```{r}
library(invgamma)
x<-seq(from=0,to=100,length=1000)
y<-dinvgamma(x, shape=0.001, rate = 0.001, log = FALSE) # shape alpha; rate beta
plot(x,y)
```