MACRO
PoisGamP.mac y;
gamma r v;
prior mu pri;
likelihood like;
posterior post;
constants r1 v1.

# PoisGamP.mac calculates the posterior probability 
# distribution when the prior density for  mu is gamma(r,v)

# Instructions for using this macro 
# are contained in  "Introduction to Bayesian Statistics, 2nd Edition" by 
# William M. Bolstad (John Wiley & Sons)

############################################################### 
#  
#  Neither Minitab Inc. nor the author of this MACRO makes 
#  any claim or offers any warranty whatsoever with regard to 
#  the accuracy of this MACRO or its suitability for use, and 
#  Minitab Inc. and the author of this MACRO each disclaims     
#  any liability with respect thereto.            
# 
###############################################################


# MACRO NAME: PoisGamP.mac
# AUTHORS NAME AND ADDRESS 
# William M (Bill) Bolstad 
# Statistics Department 
# Waikato UNIVERSITY 
# Hamilton, New Zealand
# 
# DATE: March 2007 

mconstant n  r v vinv r1 v1 rl vl v1inv  i k1 j k2 k3 k4 m1  s2 s1 ku ku1
mconstant kalpha ysum kint
mcolumn y theta pri like  post dist mu sub y1 invf
mmatrix 

# y is vector of observations 
# n is number of observations
# mu and theta contains possible values.  pri is the prior probability of each value.
# like is the likelihood function given the observed value k
name mu "mu" pri "prior" like "likelihood" post "posterior" y "y"
let n=n(y)
let ysum=sum(y)

note "sample size and sample mean"
print n ysum
let rl=ysum+1
let vl=n
if (v > 0)		#Proper gamma prior
Note "Minitab uses gamma(r, 1/v) prior which is gamma(r,v) in our notation"

let vinv=1/v
invcdf .9999 ku;
gamma r  vinv.
print ku

note "99.99 th percentile of prior distribution"
let k2=ku/10000

set mu
k2:ku/k2
end

pdf mu like;
gamma rl vl.


let r1=r+ysum
let v1=v+n
let v1inv=1/v1

pdf mu pri;
gamma r vinv.
pdf mu post;
gamma r1 v1inv.

invcdf .9999 ku;
gamma r1 v1inv.

note "gamma posterior"
print r1 v1 

print ku
note "99.99th percentile of posterior distribution"



plot pri*mu post*mu;
overlay;
connect;
 Title "Shape of gamma prior and posterior for Poisson mean";
scale 1;
max ku;
scale 2;
ldisplay 1 0 0 0.

elseif (v=0)	#improper gamma prior
let r1=r+ysum
let v1=v+n
note "gamma posterior"
print r1 v1
let v1inv=1/v1
invcdf .9999 ku;
gamma r1 v1inv.
print ku
note "99.99th percentile of posterior distribution"
let k4=ku/10000
set mu
k4:ku/k4
end
let pri=mu**(r-1)* exp(-v*mu)
let kint=(2*sum(pri)-pri(1)-pri(10000))*k4/2

pdf mu like;
gamma rl vl.


let pri=pri/kint
pdf mu post;
gamma r1 v1inv.

let ku1=max(post)/.9

plot pri*mu post*mu;
overlay;
connect;
 Title "Shape of improper gamma prior and posterior for Poisson mean";
scale 1;
max ku;
scale 2;
max ku1;
ldisplay 1 0 0 0.


endif


endmacro