-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path3 Single parameter models (pp.39-44).txt
121 lines (83 loc) · 3.58 KB
/
3 Single parameter models (pp.39-44).txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#Chapter 3 (pp. 39-44)
#Single parameter models
#Normal Distribution with Known Mean (mean=0) but unknown Variance
library(LearnBayes)
data(footballscores)
head(footballscores)
attach(footballscores)
#the goal is to estimate variance (sigma^2)
#nonnformative prior density is proportional to 1/sigma^2
#this is assigned to the variance
#the technique is to define precision parameter P=1/sigma^2
#P is distributed as U/v, where U has a chi-squared distribution with n degrees of freedom
#we are interested in a point estimate with 95% probability interval
# d is the difference variable
d=favorite-underdog-spread
# n is the sample size
n=length(d)
# v is the sum of squares of differences
v=sum(d^2)
#We simulate 1000 values from the posterior distribution
#of the standard deviation in 2 steps
#Step 1: simulate values of the precision parameter P from the scaled chi-square(n)
p=rchisq(1000,n)/v
#Step 2: Peform transformation
s=sqrt(1/p)
#Get histogram of draws of sigma
hist(s,main="")
#quantile command is used to extract the 2.5%, 50%, and 97.5%
#percentiles of this simulated sample
#A point estimate for sigma is the posterior median (50%)
#The 95%probability bounds are given by 2.5% and 97.5%
quantile(s,probs=c(0.025, 0.5, 0.975))
#Estimating a heart transplant mortality rate (lambda)
#estimate rate of success of hearth transplant
#n is number of surgeries
#y is number of deaths within 30 days
#standard MLE is Poisson distribution but this is poor when deaths is close to zero
#in thiss situation we use prior gamma(alpha, beta)
#Source of prior information is heart transplant data from a small group of hospitals
#that we believe has the same rate of mortality as the rate from the hospital of interest.
#Suppose we observe number of deaths z(j) and exposure o(j) for ten hospitals (j=1,...10)
#z(j) is poisson
#standard non-informative prior for lambda is proportional to 1/lambda
#prior distribution is then gamma(alpha = summation of z(j), beta = summation of o(j))
#or gamma(16, 15174)
#posterior is gamma(alpha+y(obs), Beta+exposure)
#y(obs) and exposure are the observed data
#model checking
#y(obs) = 1 only and exposure = 66
alpha=16;beta=15174
yobs=1;ex=66
y=0:10
#prior mean value
lam=alpha/beta
#we use dpois and dgamma to get (prior) predictive distibution
#prior mean value was used to ensure no "underflow" in computaions
#see formula in p.42
#dgamma in numerator is prior
#dgamma in denominator is posterior
py=dpois(y,lam*ex)*dgamma(lam,shape=alpha,
rate=beta)/dgamma(lam,shape=alpha+y,
rate=beta+ex)
#3 decimal places
cbind(y,round(py,3))
#draw 1000 values from posterior density of lambda
lambdaA=rgamma(1000,shape=alpha+yobs,rate=beta+ex)
#model checking
#y(obs) = 4 (NOT 1 as in previous example)
#but the rate is still low as exposure is also higher
e = 1767; yobs=4
y=0:10
#Again we see that the observed number of deaths seems consistent with this model
#since yobs = 4 is not in the extreme tails of this distribution.
py = dpois(y,lam*ex)*dgamma(lam,shape=alpha,rate=beta)/dgamma(lam, shape=alpha+y,rate=beta+ex)
cbind(y,round(py,3))
lambdaB=rgamma(1000,shape=alpha+yobs,rate=beta+ex)
par(mfrow=c(2,1))
plot(density(lambdaA), main="HOSPITAL A", xlab="lambdaA", lwd=3)
curve(dgamma(x, shape=alpha, rate=beta), add=TRUE)
legend("topright", legend=c("prior", "posterior"), lwd=c(1,3))
plot(density(lambdaB), main="Hospital B", xlab="lambdaB", lwd=3)
curve(dgamma(x, shape=alpha, rate=beta), add=TRUE)
legend("topright", legend=c("prior", "posterior"), lwd=c(1,3))