R Code 3: Normal + Inv Gamma


#posterior of mu given known sig2

t.mu=8
t.s2=4
N=6
y=rnorm(N,t.mu,sqrt(t.s2))

post.var=function(s2,s2.0) 1/((N/s2)+(1/s2.0))
post.mean=function(y,mu.0,s2,s2.0)
{
pv=post.var(s2,s2.0)
pm=pv*((sum(y)/s2)+(mu.0/s2.0))
return(pm)
}

#####################################
#analysis

s2.0=2
mu.0=0

z=seq(-10,15,.1)
plot(z,dnorm(z,mu.0,sqrt(s2.0)),typ='l',ylim=c(0,.8),col='darkgreen',lwd=2)
lines(z,dnorm(z,mean(y),sqrt(t.s2/length(y))),col='darkred',lwd=2)
lines(z,dnorm(z,post.mean(y,mu.0,t.s2,s2.0),sqrt(post.var(t.s2,s2.0))),lwd=2,col='darkblue')

##########################################
#Prior on s2

s2=seq(.1,10,.1)
plot(s2,1/s2,typ='l')

s2=seq(.01,100,.01)
plot(s2,1/s2,typ='l')

#########################################

#lets plot s2 on a log scale.

logs2=rnorm(1000,0,1)
hist(logs2)
hist(logs2,axes=F,xlab="Variance",main="")
axis(1,at=-3:3,label=10^(-3:3))
var=10^logs2
hist(var,breaks=100)

############################################
#Inverse Gamma
library('MCMCpack')

a=.01
b=.01
s2=seq(.01,10,.01)
plot(s2,dinvgamma(s2,a,b),typ='l')

post.shape=N/2+a
post.scale=sum((y-t.mu)^2)/2+b

lines(s2,dinvgamma(s2,post.shape,post.scale),col='darkblue')

###########################################

#joint mu and sigma^2

dat=rnorm(4,110,15)

post=function(mu,s2,mu0,s20,a,b) sum(dnorm(dat,mu,sqrt(s2),log=T))+dnorm(mu,mu0,sqrt(s20),log=T)+log(dinvgamma(s2,shape=1,scale=b))

mu=seq(90,110,.1)
s2=seq(40,500,2)
M=length(mu)
N=length(s2)

ans=matrix(nrow=M,ncol=N)
for (m in 1:M)
for (n in 1:N)
{
ans[m,n]=post(mu[m],s2[n],mu0=100,s20=5,a=1,b=1)
}

image(mu,sqrt(s2),exp(ans))

AttachmentSize
byNorm.pdf86.07 KB

Research Methods 3020