# com32a of September 2021,
# stk4021 exercise 14, Gott würfelt nicht: 

describe <- function(x)
{c(mean(x),sd(x),quantile(x,0.05),quantile(x,0.50),quantile(x,0.95))}

# outcomes of types 1,2,3,4,5,6: 
yy <- c(2,5,3,7,5,8) 

sim <- 10^5

keepR <- 0*(1:3) %*% t(1:5)
keepA <- 0*(1:3) %*% t(1:5)
keepB <- 0*(1:3) %*% t(1:5)
keepC <- 0*(1:3) %*% t(1:5)

## (a1)
#  then p | data is Dir[(1,1,1, 1,1,1) + yy]

xx <- 0*(1:sim) %*% t(1:6)
pp <- 0*(1:sim) %*% t(1:6)

for (j in 1:6)
{ xx[ ,j] <- rgamma(sim,1 + yy[j],1) }

for (i in 1:sim)
{ pp[i, ] <- xx[i, ] / sum(xx[i, ]) }

rhosim <- pp[ ,6] / pp[ ,1]

alphasim <- 0*(1:sim)
betasim <- 0*(1:sim) 

for (ss in 1:sim)
{
alphasim[ss] <- mean( (pp[ss, ] - 1/6)^2 )
betasim[ss] <- mean( abs(pp[ss, ] - 1/6) ) 
}

aux1 <- (1/3) * ( log(pp[ ,1])+log(pp[ ,2])+log(pp[ ,3]) )
aux2 <- (1/3) * ( log(pp[ ,4])+log(pp[ ,5])+log(pp[ ,6]) )
gammasim <- exp(aux1 - aux2) 

keepR[1, ] <- round(describe(rhosim),4)
keepA[1, ] <- round(describe(alphasim),4)
keepB[1, ] <- round(describe(betasim),4)
keepC[1, ] <- round(describe(gammasim),4)

## (a2)
#  then p | data is Dir[(3,3,3, 3,3,3) + yy]

xx <- 0*(1:sim) %*% t(1:6)
pp <- 0*(1:sim) %*% t(1:6)

for (j in 1:6)
{ xx[ ,j] <- rgamma(sim,3 + yy[j],1) }

for (i in 1:sim)
{ pp[i, ] <- xx[i, ] / sum(xx[i, ]) }

rhosim <- pp[ ,6] / pp[ ,1]

alphasim <- 0*(1:sim)
betasim <- 0*(1:sim) 

for (ss in 1:sim)
{
alphasim[ss] <- mean( (pp[ss, ] - 1/6)^2 )
betasim[ss] <- mean( abs(pp[ss, ] - 1/6) ) 
}

aux1 <- (1/3) * ( log(pp[ ,1])+log(pp[ ,2])+log(pp[ ,3]) )
aux2 <- (1/3) * ( log(pp[ ,4])+log(pp[ ,5])+log(pp[ ,6]) )
gammasim <- exp(aux1 - aux2) 

keepR[2, ] <- round(describe(rhosim),4)
keepA[2, ] <- round(describe(alphasim),4)
keepB[2, ] <- round(describe(betasim),4)
keepC[2, ] <- round(describe(gammasim),4)

## (a3)
#  then p | data is Dir[(.1,.1,.1, .1,.1,.1) + yy]

xx <- 0*(1:sim) %*% t(1:6)
pp <- 0*(1:sim) %*% t(1:6)

for (j in 1:6)
{ xx[ ,j] <- rgamma(sim,0.1 + yy[j],1) }

for (i in 1:sim)
{ pp[i, ] <- xx[i, ] / sum(xx[i, ]) }

rhosim <- pp[ ,6] / pp[ ,1]

alphasim <- 0*(1:sim)
betasim <- 0*(1:sim) 

for (ss in 1:sim)
{
alphasim[ss] <- mean( (pp[ss, ] - 1/6)^2 )
betasim[ss] <- mean( abs(pp[ss, ] - 1/6) ) 
}

aux1 <- (1/3) * ( log(pp[ ,1])+log(pp[ ,2])+log(pp[ ,3]) )
aux2 <- (1/3) * ( log(pp[ ,4])+log(pp[ ,5])+log(pp[ ,6]) )
gammasim <- exp(aux1 - aux2) 

keepR[3, ] <- round(describe(rhosim),4)
keepA[3, ] <- round(describe(alphasim),4)
keepB[3, ] <- round(describe(betasim),4)
keepC[3, ] <- round(describe(gammasim),4)

## table: 

#keepR
#     mean    sd      .05    .50     .95
#[1,] 4.5061  5.0295 1.1285 3.2448 11.7025
#[2,] 2.7429  1.8500 0.9600 2.2769  6.0343
#[3,] 7.4802 21.9690 1.2984 4.3744 21.2829

#keepA
#     mean    sd      .05    .50     .95
#[1,] 0.0070 0.0037 0.0021 0.0064 0.0140
#[2,] 0.0047 0.0027 0.0013 0.0042 0.0096
#[3,] 0.0089 0.0045 0.0029 0.0081 0.0172

#keepB
#     mean    sd      .05    .50     .95
#[1,] 0.0683 0.0190 0.0379 0.0677 0.1006
#[2,] 0.0553 0.0164 0.0296 0.0546 0.0835
#[3,] 0.0770 0.0204 0.0441 0.0766 0.1115

#keepC
#     mean    sd      .05    .50     .95
#[1,] 0.5548 0.2130 0.2738 0.5217 0.9503
#[2,] 0.6603 0.2075 0.3758 0.6320 1.0386
#[3,] 0.4811 0.2118 0.2098 0.4445 0.8756