########################################################
#           MINITAB MACRO robust.mtb                   #
# -----------------------------------------------------#
# Robust analysis of a mean using a t(v) likelihood    #
# ######################################################

# ----------------------- SETUP ------------------------
# - name five columns 'y' 'lam' 'theta' 's_lam' 's_lam2'
# - put observed data into 'y'
# - let k4 be equal to sampling variance
# - let k5 be equal to sample size
# - let k10 be equal to degrees of freedom v
# - let k11 be equal to (v+1)/2
# - initialize lam to be 1 (let 'lam'=1+0*'y')
# - initialize s_lam to be zeros (let 's_lam'=0*'y')
# - initialize s_lam2 to be zeros
# - initialize counter k1 to be equal to 0
#-------------------------------------------------------

#----------------- TO RUN GIBBS SAMPLER ----------------
# if you wish to do 500 iterations of the sampler, type
# exec 'robust' 500

# simulated stream of theta in 'theta'
# sum of lambda in 's_lam', sum of lambda**2 in 's_lam2'

let k1=k1+1

#-------------------------------- simulate lambda
let k2=sum('y'*'lam')/sum('lam')
let k3=sqrt(k4/sum('lam'))
rand 1 c100;
  normal k2 k3.
let 'theta'(k1)=c100

#-------------------------------- simulate theta
rand k5 'lam';
gamma k11 1.
let 'lam'='lam'/(.5*(('y'-'theta'(k1))**2/k4+k10))
let 's_lam'='s_lam'+'lam'
let 's_lam2'='s_lam2'+'lam'**2

prin k1

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