Code
Code
2 Aug 2021
To evaluate the posterior probability
Pr(µTue > µMon|x, α0, β0) =
∫ +∞
0
∫ +∞
µMon
f(µTue|x, α0, β0)f(µMon|x, α0, β0)dµTuedµMon.
### 2 Aug 2021
# By EM
# evaluating posterior probability
rm(list=ls())
mu.mon <- rgamma(10^4, 750, 11) mu.tue <- rgamma(10^4, 765, 11) hist(mu.mon) ## Pr(mu.tue > mu.mon|data)
mean(mu.tue>mu.mon)
To evaluate the posterior predictive probability:
Pr(x̃ > 80|x, α0, β0)
### Pr(x.tilde > 80 | data)
# 1. simulate 10^4 mus from the post. dist
# 2. simulate 10^4 new x from the likelihood
mu.mon <- rgamma(10^4, 750, 11)
x.tilde <- rpois(10^4,mu.mon)
mean(x.tilde>80)
quantile(x.tilde,.99)
Which day-of-the-week has the highest average demand?
### which DoW has the highest demand
mu.wed <- rgamma(10^4, 772, 11)
mu.thu <- rgamma(10^4, 763, 11)
mu.fri <- rgamma(10^4, 770, 11)
mu.table <- cbind(mu.mon,mu.tue,mu.wed,mu.thu,mu.fri)
1
DoW <- apply(mu.table,1,which.max)
table(DoW)/length(DoW)
## ratio tue/mon
ratio.t.m <- mu.tue/mu.mon
mean(ratio.t.m)
quantile(ratio.t.m,c(.025,.975))
2