# Lab.2 -- TSdD 2024-2025 (SOLUZIONI) -- 02 ottobre 2024

# Esercizio 1

# consultare help(dbeta)

# iperparametri
al=9.2
be=13.8

# dati
sn=15
n=20

# (a)  grafico d. a priori
# d. a priori
curve(dbeta(x,shape1=al,shape2=be), ylim=c(0,6),ylab="", xlab=expression(theta))

# (b)  probab a priori
1-pbeta(0.2,al,be)
pbeta(0.6,al,be)-pbeta(0.2,al,be)    #  0.957 

# (c) grafici
sn=15
n=20
# funzione di verosimiglianza
curve(dbeta(x,sn+1,n-sn+1),lty=2, add=TRUE)
# oppure
# ver.fun=function(theta){
#  (theta^sn)*((1-theta)^(n-sn))
# }
# curve(ver.fun(x),from=0,to=1)     # NB: non è relativa

# d. a posteriori
#curve(dbeta(x,al+sn,be+n-sn),lty=4, add=TRUE)
al.p=al+sn
be.p=be+n-sn
curve(dbeta(x,al.p,be.p),lty=4, add=TRUE)

# (d) stime puntuali e probab a posteriori
# stime puntuali
E.p=al.p/(al.p+be.p)
E.p
Mo.p=(al.p-1)/(al.p+be.p-2)
Mo.p
Me.p=qbeta(0.5,al.p,be.p)
Me.p

# probab a posteriori
al.p=al+sn
be.p=be+n-sn
1-pbeta(0.2,al.p,be.p)                            # 0.99
pbeta(0.6,al.p,be.p)-pbeta(0.2,al.p,be.p)   # 0.68

# (e) insiemi ET
L.ET=qbeta(0.025,al.p,be.p)
U.ET=qbeta(0.975,al.p,be.p)
print(c(L.ET,U.ET))

# (f) insiemi HPD
library(HDInterval)
help(hdi)

hdi(qbeta,0.95,shape1=al.p,shape2=be.p)

# per ottenere solo L.hpd e U.hpd
L.hpd=hdi(qbeta,0.95,shape1=al.p,shape2=be.p)[1]
U.hpd=hdi(qbeta,0.95,shape1=al.p,shape2=be.p)[2]
print(c(L.hpd,U.hpd))
# confronto lunghezze ET e HPD
U.ET-L.ET
U.hpd-L.hpd

# (g) Analisi non informative
al=1
be=1
al.p=al+sn
be.p=be+n-sn
al.p
be.p
# grafici
curve(dbeta(x,shape1=al,shape2=be), ylim=c(0,6),ylab="", xlab=expression(theta))
curve(dbeta(x,sn+1,n-sn+1),lty=2, add=TRUE,lwd=2)
curve(dbeta(x,sn+1,n-sn+1),lty=4, add=TRUE)
# etc etc per casa (commentare)
# analogo assumendo alpha=beta=0 
# (per casa; commentare)

# (h) maggiore informazione dai dati
# iperparametri
al=9.2
be=13.8
# dati
sn=30
n=40
# per casa: ripetere e commentare
# prior
curve(dbeta(x,shape1=al,shape2=be), ylim=c(0,6.5),ylab="", xlab=expression(theta))
# Lik
curve(dbeta(x,sn+1,n-sn+1),lty=2, add=TRUE)
# Posterior
al.p=al+sn
be.p=be+n-sn
curve(dbeta(x,al.p,be.p),lty=4, add=TRUE)



# ESERCIZIO 2

al=2
be=3
s.n=15
n=5

# (a)
# prior
curve(dgamma(x,shape=al,rate=be),xlim=c(0,5),ylab="prior density", xlab="theta",lty=3,ylim=c(0,1.5))

# lik
curve(dgamma(x,shape=s.n+1,rate=n),add=T,lty=2)

# (b)
# posterior
al.p=al+s.n
be.p=b+n

# (c)
curve(dgamma(x,shape=al.p,rate=be.p),add=T, lty=1)
E.p=al.p/be.p
E.p
V.p=al.p/be.p^2
V.p
Mo.p=(al.p-1)/be.p
Mo.p
Me.p=qgamma(0.5,al.p,rate=be.p)
Me.p

# (d)
# ET-intervals
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))

# (e) HPD
hdi(qgamma,0.95,al.p,rate=be.p)

# per ottenere solo L.hpd e U.hpd
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))
U.n-L.n
U.hpd-L.hpd

# (f)
# posterior prob's
th.t=2
pr.H0=pgamma(th.t,shape=al.p,rate=be.p) 	# Prob (Th<2|z)
pr.H1=1-pgamma(th.t,shape=al.p,rate=be.p)	# Prob(Th>2|z)
print(c(pr.H0,pr.H1))

# (g) comparisons histogram/density
M=10000
theta.simul=rgamma(M,shape=al.p, rate=be.p)
hist(theta.simul,freq=F,ylim=c(0,1.00),nclass=20)
curve(dgamma(x,shape=al.p,rate=be.p),add=T)

# I restanti punti come esercizio per casa