# this is com37a of a busy October 2021
# re analysis of a smoking-related data set from DASL;
# see Nils Exercise 22 

# note: there are various different ways of 
# parametrising these kinds of priors; 
# I follow the "most natural" of these,
# namely my own, as per Nils Exercises 22 and 23; 
# this is "nicer" and more logical than in BDA 

# this part: only y, death rates per 100,000 by lung cancer,
# for 44 states. 

#   1. CIG = Number of cigarettes smoked (hds per capita)
#   2. BLAD = Deaths per 100K population from bladder cancer
#   3. LUNG = Deathes per 100K population from lung cancer
#   4. KID = Deaths per 100K population from bladder cancer
#   5. LEUK = Deaths per 100 K population from leukemia 

data <- matrix(scan("cigarette_data2", skip=15), ncol=5, byrow=T)

# here: x = CIG, y = LUNG (death rate per thousand from lung cancer)

xx <- data[ ,1] # CIG
yy <- data[ ,3] # LUNG 

plot(xx,yy,
   xlab="cigarette consumption",ylab="lung cancer rate per 100K")

nn <- length(yy)
ybar <- mean(yy)
Q0 <- sum( (yy-ybar)^2 ) 
shat <- sqrt(Q0/(nn-1))

low0 <- sqrt((nn-1)/qchisq(0.95,nn-2))
high0 <- sqrt((nn-1)/qchisq(0.05,nn-2))

low <- shat*low0
high <- shat*high0
print(c(low,high))
# 3.636 5.226

## building the GN_2(a, b, mu0, v0) prior: 

#  first the gamma prior: 

sigmalow <- 2.2
sigmaup <- 7.7

Q <- function(para)
{
a <- para[1]
b <- para[2]
# 
hei1 <- 1 - pgamma(1/sigmaup^2, a, b) - 0.90
hei2 <- 1 - pgamma(1/sigmalow^2, a, b) - 0.10
# 
hei1^2 + hei2^2
}

nils <- nlm(Q, c(4.0, 4.0))
a0 <- nils$estimate[1]
b0 <- nils$estimate[2]

# a0 = 1.3730, b0 = 14.1507
# checking: 

# 1 - pgamma(1/sigmaup^2, a0, b0) = 0.90 
# 1 - pgamma(1/sigmalow^2, a0, b0) = 0.10
 
# then the mu given sigma part: 

# prior mean: 
mu0 <- 15.0
v0 <- 1/3.3^2 

mustar <- (v0*mu0 + nn*ybar)/(v0 + nn) # 19.6435 
dn <- 1 / (1/v0 + 1/nn)                # 0.0916   

# updated parameters: 
# GN(an, bn, mustar, vn)
# which is GN(23.3730,399.4985,19.6435,44.0918) 
# coming from GN(1.3730,14.1507,15.00,3.3) 

an <- a0 + 0.5*nn
bn <- b0 + 0.5*Q0 + 0.5*dn*(ybar-mu0)^2 
vn <- v0 + nn 

# now for simulations: 
# first sigma from sqrt of inverse gamma;
# then mu from N(mustar, sigma^2/vn), 
# then alpha = Pr( y \ge y0), y0 = 25.00
y0 <- 25.0

sim <- 10^6
# first sigma:
sigmasim <- 1/sqrt(rgamma(sim, an, bn))
# then mu for given sigma: 
musim <- mustar + (sigmasim/sqrt(vn)) * rnorm(sim)
# then alpha: 
alphasim <- 1 - pnorm((y0 - musim)/sigmasim)

keepI <- 0*(1:3) %*% t(1:2) 
keepI[1, ] <- quantile(sigmasim,c(0.025,0.975))
keepI[2, ] <- quantile(musim,c(0.025,0.975))
keepI[3, ] <- quantile(alphasim,c(0.025,0.975))

print(round(keepI,4))

#  3.4347  5.1880 sigma
# 18.3854 20.9018 mu
#  0.0458  0.1848 alpha 

muclassicL <- ybar - qt(0.975,nn-1)*shat/sqrt(nn)
muclassicR <- ybar + qt(0.975,nn-1)*shat/sqrt(nn)

#  3.6366  5.2262 for sigma, classic 
# 18.3672 20.9387 for mu, classic 

#  3.4402  5.1673 for sigma
# 18.3918 20.8939 for mu 
#  0.0459  0.1847 for alpha 

# then redoing the Bayesian analysis
# with the simpler prior that takes mu given sigma flat: 
# this corresponds to v0 = 0 above, 
# and posterior distribution GN(a0+0.5*n,b0+0.5*Q0,ybar,n).

an <- a0 + 0.5*nn
bn <- b0 + 0.5*Q0
mustar <- ybar
vn <- nn 

# sim : as above 
# first sigma:
sigmasim <- 1/sqrt(rgamma(sim, an, bn))
# then mu for given sigma: 
musim <- mustar + (sigmasim/sqrt(vn)) * rnorm(sim)
# then alpha: 
alphasim <- 1 - pnorm((y0 - musim)/sigmasim)

# then reading off credibility intervals: 

keepII <- 0*(1:3) %*% t(1:2) 
keepII[1, ] <- quantile(sigmasim,c(0.025,0.975))
keepII[2, ] <- quantile(musim,c(0.025,0.975))
keepII[3, ] <- quantile(alphasim,c(0.025,0.975))

print(round(keepII,4))

#  3.6366  5.2262 for sigma, classic 
# 18.3672 20.9387 for mu, classic 

# with informative prior for mu: 
#  3.4402  5.1673 for sigma
# 18.3918 20.8939 for mu 
#  0.0459  0.1847 for alpha 

# with flat prior for mu -- very tiny differences here: 
#  3.4358  5.1802 for sigma
# 18.3978 20.8972 for mu 
#  0.0460  0.1851 for alpha