golf<-read.table("c:/teaching/golf2.txt")
names(golf)<-c("yr","name","place","score","score2","y","money")

attach(golf)

golf$y<-golf$y/nround

burn <- 1000
length <-5000
k<-length(unique(golf$name))

theta <- matrix(72,ncol=k,nrow=(length+burn))
mu <- rep(72,length+burn)
sig2 <- rep(1,length+burn)
tau2 <- rep(1,length+burn)
ybar<-tapply(golf$y,golf$name,mean)

#k<-length(unique(golf$name))
namevec<-unique(golf$name)
ni <- numeric(k)
for(i in 1:k){
	ni[i]<-length(golf$y[golf$name==namevec[i]])
}

m.mu <- 72 
s2.mu <- 3.0
a.tau <- 3
b.tau <- .15

a.sig <- 3
b.sig <- .2

for(i in 2:(length+burn)){

# update theta (i = 1,...,k)
	for (j in 1:k){
		mustar <- (sum(golf$y[golf$name==namevec[j]])*tau2[i-1]+mu[i-1]*sig2[i-1])/(ni[j]*tau2[i-1]+sig2[i-1])
		sigstar <- sqrt((tau2[i-1]*sig2[i-1])/(ni[j]*tau2[i-1]+sig2[i-1]))
		theta[i,j] <- rnorm(1,mustar,sigstar)
	}

# update mu
	mustar <- (sum(theta[i,])*s2.mu+m.mu*tau2[i-1])/(k*s2.mu+tau2[i-1])
	sigstar <- sqrt((s2.mu*tau2[i-1])/(k*s2.mu+tau2[i-1]))
	mu[i] <- rnorm(1,mustar,sigstar)
 		
# update sigma2
	astar <- a.sig + sum(ni)/2
	bstar <- 0
	for(j in 1:k){
	   bstar <- bstar+sum((golf$y[golf$name==namevec[j]]-theta[i,j])^2)/2
	}
	bstar <- bstar+1/b.sig
	bstar <- 1/bstar
	sig2[i] <- 1/rgamma(1,astar,scale=bstar)

# update tau2
	astar <- a.tau + k/2
	bstar <- (1/b.tau + sum((theta[i,]-mu[i])^2)/2)^(-1)
	tau2[i] <- 1/rgamma(1,astar,scale=bstar)
 cat("iteration",i,"\n")
}

# plot prior predictive distribution

#plot(density(rnorm(10000,rnorm(10000,m.mu,sqrt(s2.mu)),sqrt(1/rgamma(10000,a.tau,scale=b.tau)))))
#plot posteriors for individual players
# Mike Weir
# Tiger Woods
# David Duval
# Phil Mickelson
#plot posterior predictive
# plot(density(rnorm(length-burn,mu[(burn+1):(length+burn)],sqrt(tau2[(burn+1):(length+burn)]))))
#calc posterior means
#calc posterior cov matrix