##################################################################################
#     MINITAB commands to fit balanced normal exchangeable model                 #
#--------------------------------------------------------------------------------#
#     data set:  Efron and Morris's baseball data                                #
#     observe (y1,..., y18), where yi is batting average for ith player in 45 ab #
#     let zi = log(yi/(1-yi)); assume that zi is N(theta(i), .11)                #
##################################################################################

# read in baseball data

 name c1 'x' c2 'y'					
 set c1							
 .267 .156 .311 .200 .400 .356 .330 .289 .178           
 .222 .378 .222 .244 .222 .311 .244 .222 .222           
 end

 let 'y'=log('x'/(1-'x'))				# convert data to logits
 let k1=sum(('y'-mean('y'))**2)				# k1=S

 name c3 't' c4 'pt'					# set up grid for hyperparameter tau
 set 't'                                                
 .001:1/.01                                             
 end                                                    

 let k2=.11						# k2 = sigma^2

 let 'pt'=(k2+'t'**2)**(-17/2)*exp(-k1/2/('t'**2+k2))   # computes p(tau|y) on grid
 Plot 'pt'*'t';						# and plots it
   Connect.						#

 let 'pt'='pt'/sum('pt')				# 
 name c5 'st'						# take a simulated sample
 rand 1000 'st';					# from p(tau | y)
 disc 't' 'pt'.						#
 dotplot 'st'						# graphs simulated values

 name c6 'smu'						# take a simulated sample
 rand 1000 'smu'					# from p(mu | tau, y)
 let 'smu'=mean('y')+sqrt((k2+'st'**2)/18)*'smu'	#

 store							# define a macro to
   let k12=k12+1					# simulate from p(theta(i) | mu, tau, y)
   let k11=k12+10					#
   rand 1000 ck11					#
   let ck11=ck11/sqrt(1/k2+1/'st'**2)+('y'(k12)/k2+'smu'/'st'**2)/(1/k2+1/'st'**2)
 end							#

 let k12=0						# initialize counter
 exec 18						# repeat macro 18 times to 
 							# simulate all theta(i)

 store							# find simulated sample from p(i)
   let k11=k11+1					# by taking inverse logit 
   let ck11=exp(ck11)/(1+exp(ck11))			# transformation of theta(i)
 end							#

 let k11=10						# initialize counter
 exec 18						# repeat macro 18 times

 describe c11-c28					# summarizes posteriors of p(i)

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