## Examples with Bayesian posteriors and estimates in R:

############### Section 16.2 examples ######################

############ Example 1:

sum.y <- 18
n <- 20      # Summary statistics

alpha <- 3
beta <- 1   # hyperparameters of our beta prior

# Plot of the prior (dashed) and the corresponding posterior (solid):

my.seq <- seq(0,1,length=300)  # a sequence of values on [0,1] for plotting

beta.prior <- dbeta(my.seq, alpha, beta) 
# Defining the prior
beta.posterior <- dbeta(my.seq, sum.y+alpha, n-sum.y+beta)
# Defining the posterior

plot(my.seq, beta.prior, ylim=c(-0.7,9.5), yaxt='n', 
     xlab="Values of p", ylab="Probability density", type='l', lty=3)
# Plotting the prior
lines(my.seq, beta.posterior)
# Plotting the posterior

### Point estimates for p:

posterior.mean <- (sum.y+alpha)/(alpha+beta+n)

posterior.median <- qbeta(0.50, sum.y+alpha, n-sum.y+beta)

posterior.mode <- my.seq[beta.posterior==max(beta.posterior)]

print(paste("posterior.mean=", round(posterior.mean,3), 
      "posterior.median=", round(posterior.median,3), "posterior.mode=", round(posterior.mode,3)))

### Change the prior parameters and see how the posterior (and the Bayesian estimates) change ...

############ Example 3:

# Our known sigma=15, or sigma^2 = 225.

sigma <- 15

sum.y <- 2850
n <- 30         # summary statistics

# Plot of the prior (dashed) and the corresponding posterior (solid):

my.seq <- seq(85,110,length=600)

# prior parameters for normal prior on mu:

my.delta <- 100; my.tau <- 5  # mean and sd of normal prior

prior.mu <- dnorm(my.seq, mean=my.delta, sd=my.tau)  # Defining the prior

posterior.mu <- dnorm(my.seq, mean = ( (my.delta/(my.tau^2) + sum.y/(sigma^2))/(1/(my.tau^2) + n/(sigma^2)) ), 
                sd = sqrt( (sigma^2)*(my.tau^2) / ((sigma^2) + n*(my.tau^2)) )  )
# Defining the posterior

plot(my.seq, prior.mu, ylim=c(-0.02,0.2), yaxt='n', 
     xlab="Values of mu", ylab="Probability density", type='l', lty=3)
# Plotting the prior
lines(my.seq, posterior.mu)
# Plotting the posterior

### Point estimates for mu:

posterior.mean <- ( (my.delta/(my.tau^2) + sum.y/(sigma^2))/(1/(my.tau^2) + n/(sigma^2)) )

posterior.median <- qnorm(0.50, mean = ( (my.delta/(my.tau^2) + sum.y/(sigma^2))/(1/(my.tau^2) + n/(sigma^2)) ), 
                sd = sqrt( (sigma^2)*(my.tau^2) / ((sigma^2) + n*(my.tau^2)) ) )

posterior.mode <- my.seq[posterior.mu==max(posterior.mu)]

print(paste("posterior.mean=", round(posterior.mean,3), 
      "posterior.median=", round(posterior.median,3), "posterior.mode=", round(posterior.mode,3)))

# Why are all three estimates the same (within roundoff error?)


############### Section 16.3 examples ######################


############ Example 1:
# 95% Credible interval for p:

sum.y <- 18
n <- 20      # Summary statistics

alpha <- 3
beta <- 1   # hyperparameters of our beta prior

qbeta(.025, shape1=sum.y+alpha, shape2=n-sum.y+beta)

qbeta(.975, shape1=sum.y+alpha, shape2=n-sum.y+beta)

############ Example 2:
# 90% Credible interval for theta:

qgamma(.05, shape=31, scale=0.00744)

qgamma(.95, shape=31, scale=0.00744)

# 90% Credible interval for 1/theta:

1/qgamma(.95, shape=31, scale=0.00744)

1/qgamma(.05, shape=31, scale=0.00744)

############### Section 16.4 examples ######################

### Example 1:

sum.y <- 18
n <- 20      # Summary statistics

alpha <- 3
beta <- 1   # hyperparameters of our beta prior

pbeta(0.8, shape1=sum.y+alpha, shape2=n-sum.y+beta)

### Example 2:

1 - pgamma(0.25, shape=31, scale=0.00744)