# 1= no wing cracks
# 2= detectable wing cracks
# 3= critical wing cracks

aircraft.data <- c(1,3,1,1,1,1,1,1,1,1,2,1,1,1,2,2,2,2,2,1,2,3,2,1,1,1,1,2,1,1)

table(aircraft.data)

x1<-table(aircraft.data)[1]
x2<-table(aircraft.data)[2]
x3<-table(aircraft.data)[3]

alpha1 <- 3
alpha2 <- 2
alpha3 <- 1

# This implies a prior expected value for p1 of 3/6.
# This implies a prior expected value for p2 of 2/6.
# This implies a prior expected value for p3 of 1/6.

posterior.params <- c(x1+alpha1,x2+alpha2,x3+alpha3)

# exact posterior mean for each of p1, p2, p3:

pm.alpha1 <- (x1+alpha1)/sum(posterior.params)
pm.alpha2 <- (x2+alpha2)/sum(posterior.params)
pm.alpha3 <- (x3+alpha3)/sum(posterior.params)

print(c(pm.alpha1,pm.alpha2,pm.alpha3))

# Monte Carlo approximation to help get posterior medians and credible intervals:

library(gtools)

MC.sample <- rdirichlet(n=10000, alpha=posterior.params)

apply(MC.sample, 2, mean)  # approximate posterior mean for each of p1, p2, p3

# approximate posterior medians and 95% credible intervals for each of p1, p2, p3:
quantile(MC.sample[,1], probs=c(.025, 0.5, .975))
quantile(MC.sample[,2], probs=c(.025, 0.5, .975))
quantile(MC.sample[,3], probs=c(.025, 0.5, .975))

# Sample proportions:

x1/(x1+x2+x3)
x2/(x1+x2+x3)
x3/(x1+x2+x3)


##########################################################################################

#### Redoing with a uniform prior on (p1,p2,p3):


alpha1 <- 1
alpha2 <- 1
alpha3 <- 1

# This implies a prior expected value for p1 of 1/3.
# This implies a prior expected value for p2 of 1/3.
# This implies a prior expected value for p3 of 1/3.

posterior.params <- c(x1+alpha1,x2+alpha2,x3+alpha3)

# exact posterior mean for each of p1, p2, p3:

pm.alpha1 <- (x1+alpha1)/sum(posterior.params)
pm.alpha2 <- (x2+alpha2)/sum(posterior.params)
pm.alpha3 <- (x3+alpha3)/sum(posterior.params)

print(c(pm.alpha1,pm.alpha2,pm.alpha3))

# Monte Carlo approximation to help get posterior medians and credible intervals:

library(gtools)

MC.sample <- rdirichlet(n=10000, alpha=posterior.params)

apply(MC.sample, 2, mean)  # approximate posterior mean for each of p1, p2, p3

# approximate posterior medians and 95% credible intervals for each of p1, p2, p3:
quantile(MC.sample[,1], probs=c(.025, 0.5, .975))
quantile(MC.sample[,2], probs=c(.025, 0.5, .975))
quantile(MC.sample[,3], probs=c(.025, 0.5, .975))


## Do the substantive conclusions change much?

##########################################################################################

## What if we were more certain of our prior knowledge?  Increase alpha's in prior distribution:

alpha1 <- 30
alpha2 <- 20
alpha3 <- 10

# This implies a prior expected value for p1 of 30/60.
# This implies a prior expected value for p2 of 20/60.
# This implies a prior expected value for p3 of 10/60.
# (Same as in first example)

posterior.params <- c(x1+alpha1,x2+alpha2,x3+alpha3)

# exact posterior mean for each of p1, p2, p3:

pm.alpha1 <- (x1+alpha1)/sum(posterior.params)
pm.alpha2 <- (x2+alpha2)/sum(posterior.params)
pm.alpha3 <- (x3+alpha3)/sum(posterior.params)

print(c(pm.alpha1,pm.alpha2,pm.alpha3))

# Monte Carlo approximation to help get posterior medians and credible intervals:

library(gtools)

MC.sample <- rdirichlet(n=10000, alpha=posterior.params)

apply(MC.sample, 2, mean)  # approximate posterior mean for each of p1, p2, p3

# approximate posterior medians and 95% credible intervals for each of p1, p2, p3:
quantile(MC.sample[,1], probs=c(.025, 0.5, .975))
quantile(MC.sample[,2], probs=c(.025, 0.5, .975))
quantile(MC.sample[,3], probs=c(.025, 0.5, .975))

##########################################################################################
##########################################################################################
##########################################################################################
##########################################################################################