Bayesian perspective Analysis.rmd

---
title: "STAT 8150/7150, Bayesian Data Analysis"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

1. Question 1 (18 marks)
     
PART 1.

Suppose two people, Benoit and Houying, have different prior beliefs about the average number of ER (emergency room) visits during the 10 pm - 11 pm time period. Benoit’s prior information is matched to a Gamma distribution with parameters α = 70 and β = 10, and Houying’s beliefs are matched to a Gamma distribution with α = 33.3 and β = 3.3. Let λ be the average number of visits to ER during the particular hour in the evening.

(a) (2 marks) Provide a plot representing two Gamma priors
    
    
    
```{r, message=FALSE, warning=FALSE}
library(ggplot2)
library(scales)

# Define the range of possible values for lambda
lambda <- seq(0, 20, 0.01)

# Define the parameters for Benoit's Gamma distribution
alpha_benoit <- 70
beta_benoit <- 10

# Define the parameters for Houying's Gamma distribution
alpha_houying <- 33.3
beta_houying <- 3.3

# Compute the densities for the two Gamma distributions
dens_benoit <- dgamma(lambda, shape = alpha_benoit, rate = beta_benoit)
dens_houying <- dgamma(lambda, shape = alpha_houying, rate = beta_houying)

# Combine the densities into a data frame
df_densities <- data.frame(lambda, dens_benoit, dens_houying)

# Plot the densities using ggplot2
ggplot(df_densities, aes(x = lambda)) +
  geom_line(aes(y = dens_benoit, color = "Benoit's Prior")) +
  geom_line(aes(y = dens_houying, color = "Houying's Prior")) +
  scale_color_manual(values = c("blue", "red")) +
  xlab(expression(lambda)) +
  ylab("Density") +
  ggtitle("Gamma Prior Distributions for Benoit and Houying") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 12))
```

(b) (2 marks) Compare the priors of Benoit and Houying with respect to the average value and spread. Which person believes that there will be more ER visits, on average? Which person is more confident of his/her best guess” at the average number of ER visits? Provide your reasoning
    
    The mean of a Gamma distribution with parameters alpha and beta is given by:

E(X) = alpha / beta

For Benoit's prior, we have:

     E(X) = alpha_benoit / beta_benoit 
     = 70 / 10 = 7

For Houying's prior, we have:

    E(X) = alpha_houying / beta_houying 
    = 33.3 / 3.3 = 10

    The variance of a Gamma distribution with parameters alpha and beta is given by:

Var(X) = alpha / beta^2

For Benoit's prior, we have:

    Var(X) = alpha_benoit / beta_benoit^2 
    = 70 / 10^2 = 0.7

For Houying's prior, we have:

    Var(X) = alpha_houying / beta_houying^2 
    = 33.3 / 3.3^2 = 3.03

Houying's prior has a higher mean than Benoit's prior. A higher prior mean implies that, on average, Houying believes there will be more ER visits during the 10 pm - 11 pm time period than Benoit does. On the other hand, Benoit is more confident in her best guess at the average number of ER visits, since her prior distribution has a lower variance. 

(c) (2 marks) Construct 90% interval estimates for λ using Benoit’s prior and Houying’s prior.
    

```{r}
alpha_benoit <- 70
beta_benoit <- 10

lower_benoit <- qgamma(0.05, shape = alpha_benoit, rate = beta_benoit)
upper_benoit <- qgamma(0.95, shape = alpha_benoit, rate = beta_benoit)

cat("Benoit's 90% interval estimate for λ is (", round(lower_benoit, 2), ", ", round(upper_benoit, 2), ")\n")

```

```{r}
alpha_houying <- 33.3
beta_houying <- 3.3

lower_houying <- qgamma(0.05, shape = alpha_houying, rate = beta_houying)
upper_houying <- qgamma(0.95, shape = alpha_houying, rate = beta_houying)

cat("Houying's 90% interval estimate for λ is (", round(lower_houying, 2), ", ", round(upper_houying, 2), ")\n")

```
(d) (2 marks) After some thought, Benoit believes that his best prior guess at λ is correct, but he is less confident in this guess. Explain how Benoit can adjust the parameters of his Gamma prior to reflect this new prior belief.

To adjust the parameters of Benoit's Gamma prior to reflect his new prior belief, Benoit can use a method called "Bayesian updating" which involves combining his prior belief (represented by the Gamma distribution) with new prior belief .Since Benoit now believes that his best prior guess at λ is correct, he can adjust his prior to have a lower variance, reflecting his decreased uncertainty in his estimate. To achieve this, he can increase the value of the shape parameter α, which will increase the weight placed on the mean of the distribution and decrease the spread. Alternatively, he can decrease the value of the rate parameter β, which will also decrease the spread of the distribution.


(e) (2 marks) Houying also revisits her prior. Her best guess at the average number of ER visits is now 3 larger than her previous best guess, but the degree of confidence in this guess hasn’t changed. Explain how Houying can adjust the parameters of her Gamma prior to reflect this new prior belief.

To adjust the parameters of Houying's Gamma prior to reflect her new prior belief, she can also use Bayesian updating to combine her prior belief with new data to obtain a new posterior distribution. Since Houying's best guess at the average number of ER visits is now 3 larger than her previous best guess, she can adjust the mean of her prior by increasing it by 3. However, she wants to keep the degree of confidence in her prior guess the same, which means that the variance of her prior should remain the same.One way to achieve this is to choose new values of the shape parameter α and rate parameter β that will yield a new Gamma distribution with the same variance but with a mean that is 3 larger than the original mean.

PART 2.

A hospital collects the number of patients in the emergency room admitted between 10 pm and 11 pm for each day of a week. The following table records the day and the number of ER visits for the given day.

  Day          Number of ER visits
  Sunday         8
  Monday         6
  Tuesday        6
  Wednesday      9
  Thursday       8
  Friday         9
  Saturday       7
    
Suppose one assumes Poisson sampling for the counts and a conjugate Gamma prior with parameters α = 70 and β = 10 for the Poisson rate parameter λ.

(f) (4 marks) Given the sample shown in the Table, obtain the posterior distribution for λ through the Gamma-Poisson conjugacy. Obtain a 95% posterior credible interval for λ.

To obtain the posterior distribution for λ using the Gamma-Poisson conjugacy, we need to use Bayes' theorem and the likelihood and prior distributions.

Let Y denote the number of ER visits in a particular day, and let λ be the Poisson rate parameter. Assuming Poisson sampling, the likelihood function can be written as:

    P(Y|λ) = λ^y * exp(-λ) / y!

Using the conjugate Gamma prior with parameters α = 70 and β = 10, the prior distribution for λ can be written as:

    P(λ) = β^α / Γ(α) * λ^(α-1) * exp(-βλ), where Γ is the Gamma function.

Using Bayes' theorem, we can therefore obtain the posterior distribution for λ as:

    P(λ|Y) ∝ P(Y|λ) * P(λ)

R code for posterior distribution and the 95% posterior credible interval:

```{r}
# Define the parameters for the Gamma prior distribution
alpha <- 70
beta <- 10

# Define the data
Y <- c(8, 6, 6, 9, 8, 9, 7)

# Calculate the total number of ER visits for the week
n <- sum(Y)

# Calculate the posterior parameters
alpha_post <- alpha + n
beta_post <- beta + sum(Y)

# Define a range of values for lambda
lambda_range <- seq(0, 30, 0.01)

# Calculate the posterior distribution
posterior <- dgamma(lambda_range, shape = alpha_post, rate = beta_post)

# Find the 95% posterior credible interval
lower <- qgamma(0.025, shape = alpha_post, rate = beta_post)
upper <- qgamma(0.975, shape = alpha_post, rate = beta_post)

# Print the results
cat("Posterior distribution parameters: alpha =", alpha_post, "beta =", beta_post, "\n")
cat("95% posterior credible interval for lambda: [", lower, ",", upper, "]")

```
(g) (2 marks) Suppose a hospital administrator states that the average number of ER visits during any evening hour does not exceed 6. By computing a posterior probability, evaluate the validity of the administrators statement.


```{r}
# Poisson counts of ER visits
y <- c(8, 6, 6, 9, 8, 9, 7)

# Prior hyperparameters for Gamma distribution
alpha <- 70
beta <- 10

# Compute posterior hyperparameters for Gamma distribution
alpha_post <- alpha + sum(y)
beta_post <- beta + length(y)

# Posterior distribution for lambda (Poisson rate parameter)
lambda_post <- rgamma(10000, alpha_post, beta_post)

# Posterior probability that lambda <= 6
prob_lambda_le_6 <- pgamma(6, alpha_post, beta_post)

# Print the posterior probability
print(paste0("Posterior probability that the average number of ER visits during any evening hour is less than or equal to 6 is ", round(prob_lambda_le_6, 3)))

```

A posterior probability of 0.024 means that, given the data and the prior, there is only a 2.4% chance that the average number of ER visits during any evening hour is less than or equal to 6. This is a low probability, and it suggests that the administrator's statement is unlikely to be valid based on the observed data.


(h) (2 marks) The hospital is interested in predicting the number of ER visits between 10pm and 11pm for another week. Use simulations to generate posterior predictions of the number of ER visits for another week (seven days).


```{r}
# Poisson counts of ER visits
y <- c(8, 6, 6, 9, 8, 9, 7)

# Prior hyperparameters for Gamma distribution
alpha <- 70
beta <- 10

# Compute posterior hyperparameters for Gamma distribution
alpha_post <- alpha + sum(y)
beta_post <- beta + length(y)

# Posterior distribution for lambda (Poisson rate parameter)
lambda_post <- rgamma(10000, alpha_post, beta_post)

# Generate posterior predictions for another week
predicted_ER_visits <- rpois(7, lambda_post)
days <- c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")
ER_visits_pred <- data.frame(days, predicted_ER_visits)
ER_visits_pred
```

2. Question 2 (22 marks)

PART 1. 

The Exponential distribution is often used as a model to describe the time between events, such as traffic accidents. A random variable Y has an Exponential distribution if its pdf is as follows.

    f(y | λ) = (λ exp(−λy), if y ≥ 0
               0, if y < 0
               
Here, the parameter λ > 0, considered as the rate of event occurrences. This is a oneparameter model.

(a) (4 marks) Use the prior distribution λ ∼ Gamma(a, b), and find its posterior distribution π (λ | y1, · · · , yn), where yi i.i.d. ∼ Exponential(λ) for i = 1, · · · , n. Do you recognise a known distribution? Is the Gamma prior a conjugate prior for this model?

The prior distribution is λ ~ Gamma(a, b), with pdf: 

     π(λ) = b^a λ^(a-1) exp(-bλ) / Γ(a) where Γ is the Gamma function.

The likelihood function for a sample of n independent and identically distributed Exponential random variables with parameter λ is given by:
 
    L(λ|y1, …, yn) = λ^n exp(-λ(y1 + … + yn))

Using Bayes' theorem, the posterior distribution is proportional to the product of the prior and likelihood:

    π(λ | y1, …, yn) ∝ π(λ) L(λ|y1, …, yn)

Substituting the expressions for the prior and likelihood, we get:

    π(λ | y1, …, yn) ∝ λ^(a+n-1) exp(-λ(b+y1+…+yn))

This is the kernel of a Gamma distribution with parameters a' = a + n and b' = b + y1 + … + yn. Therefore, the posterior distribution is:

    π(λ | y1, …, yn) = Gamma(a', b')

So, the posterior distribution is also a Gamma distribution with updated parameters a' = a + n and b' = b + y1 + … + yn. This is a known distribution, and it is the conjugate prior for the Exponential distribution. Therefore, the posterior distribution is:

    π(λ | y1, …, yn) = Gamma(a + n, b + y1 + … + yn)
    
The Gamma distribution is a known conjugate prior distribution for the Exponential likelihood, hence the posterior distribution is also a Gamma distribution with updated parameters.

(b) (3 marks) Define the Jeffreys prior for the parameter λ. Is it proper or improper prior? (Justify your claim)

The Jeffreys prior for the parameter λ is a non-informative prior that is invariant under reparameterization. For the Exponential distribution, the Jeffreys prior is proportional to the square root of the Fisher information:

    π_J(λ) ∝ √I(λ) = √(E[(d/dλ) ln f(Y | λ)]^2)
    = √(E[(d/dλ) (n ln λ - λ ∑ yi)]^2)
    = √(n/λ^2)

where Y = (y1, ..., yn) is the observed data. Therefore, the Jeffreys prior for λ is:

    π_J(λ) ∝ 1/λ

This is an improper prior, because it does not integrate to a finite value over the entire parameter space. Specifically, the integral of 1/λ over the positive real line diverges:

(c) (3 marks) Suppose 10 times between traffic accidents are collected: 1.5, 15, 60.3,30.5, 2.8, 56.4, 27, 6.4, 110.7, 25.4 (in minutes). With the posterior distribution derived in part (a), use Monte Carlo approximation to calculate the posterior mean,median, and middle 95% credible interval for the rate λ.

We know from part (a) that the posterior distribution of λ is a Gamma distribution with shape parameter a' = a + n = 10 and scale parameter b' = b + ∑ i=1 to n yi = 334.6. Therefore, the posterior distribution is:

     π(λ | y1, ..., yn) = Gamma(λ | 10, 334.6)

```{r, message=FALSE, warning=FALSE}
# Observed  data
y <- c(1.5, 15, 60.3, 30.5, 2.8, 56.4, 27, 6.4, 110.7, 25.4)

# Prior parameters
a <- 70
b <- 10

# Posterior parameters
n <- length(y)
a_post <- a + n
b_post <- b + sum(y)

# Draw samples from the posterior distribution
library(MCMCpack)
set.seed(123)
samples <- rgamma(100000, shape = a_post, rate = b_post)

# Posterior mean
post_mean <- mean(samples)
cat("Posterior mean: ", post_mean, "\n")

# Posterior median
post_median <- median(samples)
cat("Posterior median: ", post_median, "\n")

# Middle 95% credible interval
lower <- quantile(samples, probs = 0.025)
upper <- quantile(samples, probs = 0.975)
cat("Middle 95% credible interval: (", lower, ", ", upper, ")\n")
```

(d) (2 marks)  Use Monte Carlo approximation to generate another set of 10 predicted times between events.


```{r}
# posterior distribution parameters (from part c)
a_post <- 11  # shape parameter
b_post <- 337  # scale parameter

# generate 10,000 random samples from the posterior distribution
samples <- rgamma(10000, a_post, scale = 1/b_post)

# generate 10 predicted times between events
set.seed(123)  # for reproducibility
predicted_times <- rexp(10, rate = sample(samples, 1))

# Results in a table
results <- data.frame(Predicted_Time_Between_Events = predicted_times)
rownames(results) <- NULL
print(results)
```

PART 2. 
      
The Weibull distribution is often used as a model for survival times in biomedical, demographic, and engineering analyses. A random variable Y has a Weibull distribution if its pdf is as follows.

    f(y | α, λ) = λαyα−1 exp (−λyα) for y > 0.

Here, α > 0 and λ > 0 are parameters of the distribution. For this problem, assume that α = α0 is known, but λ is not known, i.e. a simplified case of a one-parameter model. Also, assume that software routines for simulating from Weibull distributions are available (e.g., rweibull())

(e) (3 marks) Assuming a prior distribution π (λ | α = α0) ∝ 1, find its posterior π (λ | y1, . . . , yn, α = α0), where yi
i.i.d. ∼ Weibull (λ, α = α0) for i = 1, · · · , n. Write the name of the distribution and expressions for its parameter values.

The posterior distribution of λ can be obtained by applying Bayes' theorem:

    π (λ | y1, . . . , yn, α = α0) ∝ f(y1, . . . , yn | λ, α = α0) π(λ | α = α0)

where f(y1, . . . , yn | λ, α = α0) is the likelihood function of λ given the observed data y1, . . . , yn, and α = α0 is the known value of α.

Since the observations are independent and identically distributed (i.i.d.), the likelihood function is the product of the individual likelihoods:

    f(y1, . . . , yn | λ, α = α0) = ∏i=1n λα0yiα0−1 exp(−λyiα0)

Using the prior distribution π(λ | α = α0) ∝ 1, we can write:

    π (λ | α = α0) = 1/K where K is a normalizing constant.

Then, the posterior distribution is proportional to the product of the likelihood and prior:

     π (λ | y1, . . . , yn, α = α0) ∝ ∏i=1n λα0yiα0−1 exp(−λyiα0) * 1

    π (λ | y1, . . . , yn, α = α0) ∝ λ^{nα0} * exp(-λsum(yi^α0))

This is a Gamma distribution with shape parameter k = nα0 and rate parameter θ = sum(yi^α0). So the posterior distribution is:

    Gamma(λ | k = nα0, θ = sum(yi^α0))

Therefore, the name of the distribution is Gamma distribution with parameters k = nα0 and θ = sum(yi^α0).


(f) (5 marks) Using the posterior distribution derived in part (e), explain step-by-step how you would use Monte Carlo simulation to approximate the posterior median survival time, assuming that α = α0. Write also a pseudo-code in R.

To use Monte Carlo simulation to approximate the posterior median survival time, we can follow these steps:

1.Draw a large number of samples from the posterior distribution of λ given the data y1, . . . , yn and α = α0, which is a Gamma distribution with shape parameter a = n and rate parameter b = ∑i=1^n yiα0.

2For each sample, simulate a single value of Y from the Weibull distribution with the corresponding value of λ and α0.

Calculate the empirical median of the simulated Y values across all the samples drawn in step 1.


```{r}
# Collected data
y <- c(1.5, 15, 60.3, 30.5, 2.8, 56.4, 27, 6.4, 110.7, 25.4)

# Posterior parameters
n <- length(y)
a_post <- n
b_post <- sum(y)

# Draw samples from the posterior distribution
set.seed(123)
samples <- rgamma(10000, shape = a_post, rate = b_post)

# Simulate survival times from the Weibull distribution
weibull_sim <- function(lambda, alpha) {
  rweibull(1, shape = alpha, scale = 1/lambda)
}
simulated_times <- sapply(samples, weibull_sim, alpha = a_post)

# Calculate the posterior median survival time
posterior_median <- median(simulated_times)
posterior_median
```

(g) (2 marks) What family of distributions represents the conjugate prior distributions for λ, assuming that α = α0.

The family of gamma distributions represents the conjugate prior distributions for λ, assuming that α = α0 in the Weibull distribution. Specifically, if we assume a gamma prior distribution for λ, with shape parameter a and rate parameter b, then the posterior distribution of λ given the data y has a gamma distribution with updated shape and rate parameters:

     shape parameter = a + n*α0
     rate parameter = b + sum(y^α0)

This means that if we use a gamma distribution as our prior distribution for λ, the resulting posterior distribution will also be a gamma distribution. This makes the gamma distribution a conjugate prior for λ in the Weibull distribution with fixed shape parameter α0.