--- title: "Simple MCMC chains" author: Jeff Rouder date: Sept. 28th, 2018 output: html_document --- # Motivating Problem Same problem as before: Model: \begin{aligned} Y|\mu,\sigma^2 &\sim \mbox{Normal}(\mu,\sigma^2)\\ \mu &\sim \mbox{Normal}(a,b)\\ \sigma^2 &\sim \mbox{Inverse Gamma}(r,s) \end{aligned} Data: $y=105,100,124,104$ Desire: $\mu|Y$, $\sigma^2|Y$ (marginals) Know: $\mu|Y,\sigma^2$, $\sigma^2|Y,\mu$ Example: {r} library(MCMCpack) dat=c(105,100,124,104) N=length(dat) mean.dat=mean(dat) #prior settings a=100 b=5^2 r=2 s=100 #graph of priors mu=seq(70,130,.1) plot(mu,dnorm(mu,a,sqrt(b)),typ='l') s2=seq(1,500,1) plot(s2,dinvgamma(s2,shape=r,scale=s),typ='l')  {r} M=1000 mu=1:M sig2=1:M #really bad starting values mu[1]=100 sig2[1]=100 for (m in 2:M) { v=1/((N/sig2[m-1]+1/b)) c=N*mean.dat/sig2[m-1]+a/b mu[m]=rnorm(1,c*v,sqrt(v)) s2.scale=sum((dat-mu[m])^2)/2+s sig2[m]=rinvgamma(1,shape=r+N/2,scale=s2.scale) } plot(mu,typ='l') abline(h=mean(dat),col='red',lwd=2) abline(h=mean(mu),col='blue',lwd=2) plot(sig2,typ='l') plot(sqrt(sig2),typ='l') abline(h=sd(dat),col='red',lwd=2) abline(h=mean(sqrt(sig2)),col='blue',lwd=2)  Your Turn (but wait) * Group 1, Explore the effect of different priors for this data as well as from data created from {r} seed=123 dat=rnorm(20,106,10)  * Group 2, Explore the effect of different starting values. Choose some really bad ones * Group 3, Figure out this model in JAGS * Group 4, Figure out this model in STAN # Two groups Cell Model: \begin{aligned} X_i|\mu_x,\sigma^2 &\sim \mbox{Normal}(\mu_x,\sigma^2)\\ Y_i|\mu_y,\sigma^2 &\sim \mbox{Normal}(\mu_y,\sigma^2)\\ \mu_x &\sim \mbox{Normal}(a,b)\\ \mu_y &\sim \mbox{Normal}(a,b)\\ \sigma^2 &\sim \mbox{Inverse Gamma}(r,s) \end{aligned} Effects Model: \begin{aligned} X_i|\mu,\alpha,\sigma^2 &\sim \mbox{Normal}(\mu+\alpha/2,\sigma^2)\\ Y_i|\mu_,\alpha,\sigma^2 &\sim \mbox{Normal}(\mu-\alpha/2,\sigma^2)\\ \mu &\sim \mbox{Normal}(a_1,b_1)\\ \alpha &\sim \mbox{Normal}(a_2,b_2)\\ \sigma^2 &\sim \mbox{Inverse Gamma}(r,s) \end{aligned} Here is some data for response times in a Stroop task. {r} set.seed(124) y=rnorm(35,.5,.15) x=rnorm(35,.54,.15)  * Group 5, Implement these model in STAN. Which do you prefer in real life? * Group 6, Derive full conditionals for effects model. Implement them in a chain in R.