# Lab. 4 TSdD 2024 -- 15 ottobre 2024


# Esercizio  (modello Poisson-gamma, approx normale)

al=2
be=3
sn=15
n=5
al.p=al+sn
be.p=be+n

# (a) densità a posteriori esatta
curve(dgamma(x,al.p,rate=be.p),from=0,to=5, xlab=expression(theta),ylab="densità a posteriori")

# (b)
th.tilde=(al.p-1)/be.p
I.B=be.p^2/(al.p-1)
sd.tilde=sqrt(1/I.B)

curve(dnorm(x,th.tilde,sd.tilde),lty=2,add=TRUE)

# (c)
# Insiemi esatti
# ET
L.n=qgamma(0.025,shape=al.p,rate=be.p)
U.n=qgamma(0.975,shape=al.p,rate=be.p)
print(c(L.n,U.n))

#  HPD
library(HDInterval)
# help(hdi)

L.hpd=hdi(qgamma,0.95,shape=al.p,rate=be.p)[1]
U.hpd=hdi(qgamma,0.95,shape=al.p,rate=be.p)[2]
print(c(L.hpd,U.hpd))


# Approssimati
gg=0.05
L.tilde=th.tilde-qnorm(1-gg/2)*sd.tilde
U.tilde=th.tilde+qnorm(1-gg/2)*sd.tilde
c(L.tilde,U.tilde)

# P(Theta<2|z.n)

# esatta
pgamma(2,al.p,rate=be.p)

# approx
pnorm(2,th.tilde,sd.tilde)


# 2. MC inferenza a posteriori (Beta-Binomial)
al=2
be=2
s.n=4
n=10

al.p=al+s.n
be.p=be+n-s.n

# (a) Exact
post.mean=al.p/(al.p+be.p)
post.prob=pbeta(0.3,al.p,be.p)
post.quant=qbeta(0.95,al.p,be.p)
print(c(post.mean,post.prob,post.quant))      #  0.4285714 0.1653975 0.6452007


# (b) MC approximations
M=10000
th.MC=rbeta(M,al.p,be.p)

post.mean.MC=mean(th.MC)
post.prob.MC=length(th.MC[th.MC<0.3])/M
post.quant.MC=quantile(th.MC,0.95)
print(c(post.mean.MC,post.prob.MC,post.quant.MC))  # 0.4300297 0.1622000 0.6453676 

# (c) Plot
hist(th.MC,freq=F,xlim=c(0,1))
curve(dbeta(x,al.p,be.p), add=T)  

# (d) ET intervals
ln=qbeta(0.025,al.p,be.p)
un=qbeta(0.975,al.p,be.p)
ln                                                 # 0.19
un                                                # 0.68

ln.MC=quantile(th.MC,0.025)
un.MC=quantile(th.MC,0.975)
ln.MC                                         # 0.19
un.MC                                        # 0.68

# (e) Invariance of ETs                 # (-0.14 0.74)
# psi=log(th/(1-th))
ln.psi.MC=log(ln.MC/(1-ln.MC))
un.psi.MC=log(un.MC/(1-un.MC))
ln.psi.MC
un.psi.MC

psi.MC=log(th.MC/(1-th.MC))
quantile(psi.MC,0.025)
quantile(psi.MC,0.975)