# com32b of September 2021, Gott würfelt nicht,
# with prior being a mixture of point mass and Dirichlet.

# below I resort for simplicity to the approximation trick
# portraying "degenerate prior at p0 = (1/6,...,1/6)"
# as almost the same as Dir(cc/6,...,cc/6), with cc = big.

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

# counts observed: 
yy <- c(2, 5, 3, 7, 5, 8) 

# may also try out this with other data,
# and e.g. bigger sample size: 
# yy <- c(2, 2, 2, 8, 8, 8)

nn <- sum(yy) 

# prior mix: 
# w0 pi0 + w1 pi1 with Dir(cc p0) and Dir(6 p0), 
# p0 = (1/6, ..., 1/6) 
# when cc is Inf this means pi0 is totally sure of p0,
# i.e. the die is a fair one. 
# I write this as w0 Dir(cc0 p0) + w1 Dir(cc1 p0),
# where cc0 is big and cc1 = 6

cc0 <- 10^4
cc1 <- 6 

p0 <- 1/6 + 0*(1:6) 
w0 <- 0.5
w1 <- 1 - w0 

hei1 <- lgamma(cc0) - 6*lgamma(cc0/6) 
hei2 <- sum(lgamma(cc0/6+yy)) - lgamma(cc0 + nn)
logf0 <- hei1 + hei2 

hei3 <- lgamma(cc1) - 6*lgamma(cc1/6) 
hei4 <- sum(lgamma(cc1/6+yy)) - lgamma(cc1 + nn)
logf1 <- hei3 + hei4

aux0 <- logf0 - mean(c(logf0,logf1))
aux1 <- logf1 - mean(c(logf0,logf1))

# print(c(logf0, logf1))
# print(c(aux0,aux1))
# these have identical differences 

post0 <- w0*exp(aux0) / (w0*exp(aux0) + w1*exp(aux1))
post1 <- w1*exp(aux1) / (w0*exp(aux0) + w1*exp(aux1))

print( round(c(post0,post1),4) ) 
# 0.9173 and 0.0827, posterior probabilities 
# associated with priorpart0 and priorpart1 
# with the original data set given 

# all later posterior distributions: 
# post0 * pi0(p | y) + post1 * pi1(p | y) 
# and under pi0(p | y), the p is simply fixed at p0.

# I first create posterior realisations from the 
# pi1(p | y) part: 

sim <- 10^4
# may use sim = 10^5 if seriously needed 

# first the gamma variables, column by column: 
zz <- 0*(1:sim) %*% t(1:6)
for (j in 1:6)
{ zz[ ,j] <- rgamma(sim,1 + yy[j]) }

# then using these to create the pp, row by row: 
pp <- 0*zz
for (ss in 1:sim)
{ pp[ss, ] <- zz[ss, ] / sum(zz[ss, ]) }

# 
rhoprelimsim <- 0*(1:sim)
alphaprelimsim <- 0*(1:sim) 
betaprelimsim <- 0*(1:sim) 
gammaprelimsim <- 0*(1:sim) 

for (ss in 1:sim)
{
rhoprelimsim[ss] <- pp[ss,6]/pp[ss,1]
# 
aux11 <- pp[ss, ] 
alphaprelimsim[ss] <- mean( (aux11-1/6)^2 )
betaprelimsim[ss] <- mean( abs(aux11-1/6) )
oi1 <- (aux11[4]*aux11[5]*aux11[6])^(1/3) 
oi2 <- (aux11[1]*aux11[2]*aux11[3])^(1/3)
gammaprelimsim[ss] <- oi1/oi2 
}

# then to the final thing: 
rhosim <- 0*(1:sim)
alphasim <- 0*(1:sim) 
betasim <- 0*(1:sim) 
gammasim <- 0*(1:sim) 

ok <- 0*(1:sim)
uu <- runif(sim) 

for (ss in 1:sim)
{
ok[ss] <- 1*(uu[ss] <= post0)
rhosim[ss] <- ok[ss]*1 + (1-ok[ss])*rhoprelimsim[ss]
alphasim[ss] <- ok[ss]*0 + (1-ok[ss])*alphaprelimsim[ss]
betasim[ss] <- ok[ss]*0 + (1-ok[ss])*betaprelimsim[ss]
gammasim[ss] <- ok[ss]*1 + (1-ok[ss])*gammaprelimsim[ss]
}

# can then investigate these posterior distributions,
# finding posterior means, medians, quantiles, etc. 

# plot(density(gammasim))
plot(gammasim)

# 
print( round(describeme(rhosim),4) )
print( round(describeme(alphasim),4) )
print( round(describeme(betasim),4) )
print( round(describeme(gammasim),4) )

# with sim = 10^5 I find: 
# mean   sd     5%     50%    95% 
# 1.2914 1.7046 1.0000 1.0000 2.7404 rho 
# 0.0006 0.0022 0.0000 0.0000 0.0056 alpha 
# 0.0057 0.0196 0.0000 0.0000 0.0631 beta
# 1.0902 0.3848 1.0000 1.0000 1.7511 gamma