## com27a of september 2021,
#  stk 4021, nils exercise 12, alarm or not

#  risk0 = 5000*Pr(theta > 0.15 | data)
#  risk1 = 1000*Pr(theta < 0.15 | data) 

#  so i decide "alarm" if 1000 G(0.15) < 5000 (1-G(0.15))
#  or G(0.15) < 5/6 = 0.8333
#  with G(0.15) = Pr(theta < 0.15 | data) 

nn <- 50
yval <- seq(0,nn,by=1) 

## (a) theta is unif: 
#  theta | y is Beta(1 + y, 1 + n-y) 

## (b) theta is Beta(2,8) 
#  theta | y is Beta(2 + y, 8 + n-y) 

G1val <- pbeta(0.15, 1+yval, 1+nn-yval)
G2val <- pbeta(0.15, 2+yval, 8+nn-yval)

cbind(yval,G1val,G2val)[1:10, ]

# read off from this: 
# prior1: alarm if y \ge 5
# prior2: alarm if y \ge 5 

## (c) now more complicated, prior is 0.5 Beta(2,8) + 0.5 Beta(8,2)  
#  posterior is r_1(y) be(theta,2+y,8+n-y) + r_2(y) be(theta,8+y,2+n-y)
#  where r_1(y) = f_1(y) / sum and r_2(y) = f_2(y) / sum 

logf1 <- 0*yval
logf2 <- 0*yval

for (j in 1:length(yval))
{
y <- yval[j]
# 
aux1 <- lgamma(nn+1)-lgamma(y+1)-lgamma(nn-y+1)
aux2 <- lgamma(2+y)+lgamma(8+nn-y)-lgamma(nn+2+8) + lgamma(2+8)-lgamma(2)-lgamma(8)
aux3 <- lgamma(8+y)+lgamma(2+nn-y)-lgamma(nn+2+8) + lgamma(2+8)-lgamma(8)-lgamma(2)
logf1[j] <- aux1 + aux2 
logf2[j] <- aux1 + aux3
} 

# checking: 
sum(exp(logf1)) # 1
sum(exp(logf2)) # 1 

# probabilities for Type 1 and Type 2 given data y: 
ww1 <- exp(logf1) / ( exp(logf1) + exp(logf2) ) 
ww2 <- exp(logf2) / ( exp(logf1) + exp(logf2) ) 

matplot(yval,cbind(ww1,ww2),type="p")

G3val <- ww1*pbeta(0.15,2+yval,8+nn-yval) + ww2*pbeta(0.15,8+yval,2+nn-yval) 

letscheck <- cbind(yval,cbind(G1val,G2val,G3val))
print(round(letscheck[1:10, ],4))

#      yval  G1val  G2val  G3val
# [1,]    0 0.9997 0.9992 0.9992
# [2,]    1 0.9975 0.9956 0.9956
# [3,]    2 0.9875 0.9833 0.9833
# [4,]    3 0.9587 0.9531 0.9531
# [5,]    4 0.8978 0.8944 0.8944
# [6,]    5 0.7967 0.8012 0.8011 so from here on, G(0.15 | data) < 0.8333 
# [7,]    6 0.6600 0.6766 0.6766
# [8,]    7 0.5049 0.5337 0.5336
# [9,]    8 0.3543 0.3909 0.3907
#[10,]    9 0.2274 0.2648 0.2646

# boring answer: for each of priors a, b, c,
# take "alarm" for y \ge 5