y <- scan("faculty.dat")

# priors

m.mu <- 5
s2.mu <- 1
a.sig <- 6
b.sig <- 1/10


# starting values

M<-10000
mu<-sig2<-rep(0,M)
mu[1] <- 5
sig2[1] <- 2

n <- length(y)
ybar <- mean(y)

for(i in 2:M){

# generate a value from the complete conditional
# for mu

mustar <- (n*ybar*s2.mu+m.mu*sig2[i-1])/(sig2[i-1]+n*s2.mu)
sigstar <- sqrt(sig2[i-1]*s2.mu/(n*s2.mu+sig2[i-1]))
mu[i] <- rnorm(1,mustar,sigstar)

# generate a value from the complete conditional
# for sigma2

astar <- n/2+a.sig
bstar <- 1/(sum((y-mu[i])^2)/2 + 1/b.sig)
sig2[i] <- 1/rgamma(1,astar,scale=bstar)
}
output<-cbind(mu,sig2)
return(output)