golfmcmc<-function(niter=200){

#priors
m.mu<-72
s2.mu<-2
a.tau<-18
b.tau<-.015
a.sig<-83
b.sig<-.0014

#reading in the data
golf04<-read.table("golfdataR.dat",header=F)
names(golf04)<-c("id","shots","tourney")
y<-golf04$shots
players<-unique(golf04$id)
np<-length(unique(golf04$id))
nround<-table(golf04$id)

#starting values
theta<-matrix(72,ncol=np,nrow=niter)
mu<-rep(72,niter)
tau2<-rep(1,niter)
sig2<-rep(1,niter)

for(i in 2:niter){
      #generate theta's
      ssq<-0
      for(j in 1:np){
      	    mustar<-(sum(y[golf04$id==players[j]])*tau2[i-1]+mu[i-1]*sig2[i-1])/
	            (nround[players[j]]*tau2[i-1]+sig2[i-1])
            sigstar<-tau2[i-1]*sig2[i-1]/(nround[players[j]]*tau2[i-1]+sig2[i-1])
	    theta[i,j]<-rnorm(1,mustar,sqrt(sigstar))
	    ssq<-ssq+sum((y[ golf04$id==players[j]]-theta[i,j])^2)   
	    }

      #generate mu
      mustar<-(mean(theta[i,])*np*s2.mu+m.mu*tau2[i-1])/(np*s2.mu+tau2[i-1])
      sigstar<-(s2.mu*tau2[i-1])/(np*s2.mu+tau2[i-1])
      mu[i]<-rnorm(1,mustar,sqrt(sigstar))

      #generate sig2
      astar<-a.sig+sum(nround)/2
      bstar<-(1/b.sig+ssq/2)^(-1)
      sig2[i]<-1/rgamma(1,astar,scale=bstar)

      #generate tau2
      astar<-a.tau+np/2
      bstar<-(1/b.tau+.5*sum(theta[i,]-mu[i])^2)^(-1)
      tau2[i]<-1/rgamma(1,astar,scale=bstar)
      out<-cbind(theta,mu,tau2,sig2)
      cat("iter =",i,"\n")
 
}
     return(out)
}