# com38a of a busy October 2021,
# re analysis of a smoking-related data set from DASL.
# see Nils Exercises 22-23 stk 4021

#    CIG = Number of cigarettes smoked (hundreds per capita)
#    BLAD = Deaths per 100K population from bladder cancer
#    LUNG = Deaths per 100K population from lung cancer
#    KID = Deaths per 100K population from bladder cancer
#    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]
yy <- data[ ,3]

aux <- lm(yy ~ xx)$coef
xval <- seq(14.0,42.4,by=0.1)
linregval <- aux[1] + aux[2]*xval 

par(bty="l")
par(cex.lab=1.4142) 

plot(xx,yy, xlab="cigarettes", ylab="death rate lung cancer")
matlines(xval,linregval,type="l",lty=2,lwd=1.5,col=2)

## two states to the right, with high x, but not very high y:
#       x               y 
#  DC	40.46	5.60	27.27	3.13	7.08
#  NE	42.40	6.54	23.03	2.85	6.67 # outlier, good for them 

nn <- length(yy)
ML <- lm(yy ~ xx)$coef
res <- yy - ML[1] - ML[2]*xx
Q0 <- sum(res^2)
shat <- sqrt(Q0/(nn-2))

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

low <- shat*low0
high <- shat*high0
print(c(low,high))
#  2.606 3.745 is the standard 95% interval for sigma 

## building the GN_2(a, b, beta0, M0) prior: 

#  first the gamma prior: 

sigmalow <- 1.0
sigmaup <- 5.0

Q <- function(para)
{
a <- para[1]
b <- para[2]
# 
hei1 <- pgamma(1/sigmalow^2, a, b) - 0.90
hei2 <- pgamma(1/sigmaup^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 = 0.9395, b0 = 2.1967
# checking: 

# with low and big equal to 1.0 and 5.0: 
# a0 = 0.9395, b0 = 2.1967

# pgamma(1/sigmaup^2, a0, b0) = 0.10 
# pgamma(1/sigmalow^2, a0, b0) = 0.90
 
# then the (alpha,beta) part: 

# prior mean: 
beta0 <- c(15.0, 0.0)
dim(beta0) <- c(2,1)
dim(ML) <- c(2,1)

# prior precision: 
c0 <- 1/2.0^2
M0 <- c0 * rbind(c(1,0), c(0,1))

Mn <- 0*(1:2) %*% t(1:2)
Mn[1,1] <- 1
Mn[1,2] <- mean(xx)
Mn[2,1] <- mean(xx)
Mn[2,2] <- mean(xx^2)

betastar <- solve(M0 + nn*Mn) %*% (M0%*%beta0 + nn*Mn%*%ML)
Dn <- nn*M0 %*% solve(M0 + nn*Mn) %*% Mn

# updated parameters: 

an <- a0 + 0.5*(2+nn)
bn <- b0 + 0.5*Q0 + 0.5*t(ML-beta0) %*% Dn %*% (ML-beta0)
V <- solve(M0 + nn*Mn)
sigma1 <- sqrt(V[1,1])
sigma2 <- sqrt(V[2,2])
rho <- V[1,2]/(sigma1*sigma2) # a strong -0.9736

# now for simulations: 
# first sigma from sqrt of inverse gamma;
# then beta from N_2(betastar, sigma^2 V)

sim <- 10^5
sigmasim <- 1/sqrt(rgamma(sim, an, bn))

uu <- rnorm(sim)
vv <- rnorm(sim)

alphasim <- betastar[1] + sigma1*sigmasim*uu
betasim <- betastar[2] + sigma2*sigmasim*(rho*uu + sqrt(1-rho^2)*vv)

keep <- cbind(sigmasim, alphasim, betasim)

# 
print( round(c(a0,b0,an,bn),4) ) 

# a0      b0       an      bn 
# 1.3730  14.1507  24.3730 219.7017

# beta0:    (15, 0)
# betastar: (7.40,0.49)

## can then answer questions relating to posterior
#  for (alpha,beta,sigma) 

#  (i) alpha, beta, sigma: 
look1 <- c(mean(alphasim),quantile(alphasim,c(0.025,0.975)))
look2 <- c(mean(betasim),quantile(betasim,c(0.025,0.975)))
look3 <- c(mean(sigmasim),quantile(sigmasim,c(0.025,0.975)))

#  (ii) gamma = Pr(y \ge 25.0 | x = 35.0) 
#       = 1 - Phi((y0 - alpha - beta x0)/sigma)
x0 <- 35.0
y0 <- 25.0
gammasim <- 1 - pnorm((y0 - alphasim - betasim*x0)/sigmasim)
look4 <- c(mean(gammasim),quantile(gammasim,c(0.025,0.975)))

#  (iii) gamma = Pr(y \le y0 | x = x0) for Nebraska  
#       = Phi((y0 - alpha - beta x0)/sigma)
#  NE	42.40	6.54	23.03	2.85	6.67 # outlier, good for them 
x0 <- 42.40
y0 <- 23.03
deltasim <- pnorm((y0 - alphasim - betasim*x0)/sigmasim)
look5 <- c(mean(deltasim),quantile(deltasim,c(0.025,0.975)))

#   (iv) future Y0 for given x0: 
x45 <- 10.0 # a low value: 
Y45 <- alphasim + betasim*x45 + sigmasim*rnorm(sim)

x46 <- 50.0 # a low value: 
Y46 <- alphasim + betasim*x46 + sigmasim*rnorm(sim)

look45 <- c(mean(Y45),quantile(Y45,c(0.025,0.975)))
look46 <- c(mean(Y46),quantile(Y46,c(0.025,0.975)))

look <- rbind(look1,look2,look3,look4,look5,look45,look46)

print(round(look,4))

## with the informative GN prior given in Exercise 23: 

#look1   7.3973  3.4931 11.2894 alpha 
#look2   0.4936  0.3404  0.6472 beta
#look3   2.9925  2.4580  3.6852 sigma 
#look4   0.4581  0.2398  0.6871 gamma 
#look5   0.0557  0.0026  0.2161 delta, for Nebraska 
#look45 12.3345  5.8947 18.7763 Y45, for x45 = 10.0 low 
#look46 32.0686 24.9582 39.2052 Y46, for x46 = 30.0 high