# data 
ybar_c(2.145,7.59)
sighat_matrix(c(0.45,0.275,0.275,0.981),2,2)
n_50



#GENERATE INVERSE WISHART----------------------------------------------------
"riwish" <- function(k,df,S)
{
#
# Delivers a draw from the IW distribution in k.k dimensions, with degrees of 
# freedom df and scale matrix S 
#
# For example, in the normal random sampling model with Jeffreys' prior,
#  the IW posterior for the variance matrix is sampled by the call 
#  riwish(2,n-1,sighat*(n-k)) where sighat is the sample variance matrix
#  based on n observations
# 
R <- diag(sqrt(2*rgamma(k,(df + k  - 1:k)/2)))
R[outer(1:k, 1:k,  "<")] <- rnorm (k*(k-1)/2)
R <- t(solve(R))%*% chol(S)
return(t(R)%*%R)
}

	

nmc_2000; yn_mu_matrix(0,2,nmc)

# direct posterior and predictive samples ... 
for (i in 1:nmc) { 
  tcholsig_t(chol(  riwish(2,n-1,sighat*(n-2))  ))
  mu[,i]_ybar+(tcholsig/sqrt(n))     %*% rnorm(2)	
  yn[,i]_ybar+(tcholsig*(1+1/n)) %*% rnorm(2)	 
# n.b., you can sample the predictive alternatively via 
#   yn[,i]_mu[,i]+tcholsig %*% rnorm(2)	 
# as some of you did, very neatly 
} 

	
par(mfrow=c(4,1),bty='n')

hist(mu[1,],nclass=30,probability=T)
v_sqrt(sighat[1,1]/n)
x_seq(ybar[1]-3*v,ybar[1]+3*v,length=100)
lines(x,dt((x-ybar[1])/v,n-1)/v,col=3) 

hist(mu[2,],nclass=30,probability=T)
v_sqrt(sighat[2,2]/n)
x_seq(ybar[2]-3*v,ybar[2]+3*v,length=100)
lines(x,dt((x-ybar[2])/v,n-1)/v,col=3) 

hist(yn[1,],nclass=30,probability=T)
v_sqrt(sighat[1,1]*(1+1/n))
x_seq(ybar[1]-3*v,ybar[1]+3*v,length=100)
lines(x,dt((x-ybar[1])/v,n-1)/v,col=3) 

hist(yn[2,],nclass=30,probability=T)
v_sqrt(sighat[2,2]*(1+1/n))
x_seq(ybar[2]-3*v,ybar[2]+3*v,length=100)
lines(x,dt((x-ybar[2])/v,n-1)/v,col=3)