# Lab. 3 TSdD 2024 -- 9 ottobre 2024


# Esercizio 1 (modello normale-normale)

# (a)
# dati
x.med=130
sig=5
n=2
# parametri a priori
mu.0=120
n.0=1/4
# parametri a posteriori
mu.p=(n*x.med+n.0*mu.0)/(n.0+n)
sig.p=sig/sqrt(n.0+n)

#(b) n.0 vale 1/4 di una osservazione
n/n.0

# (c) e (d) 
prior=function(x){dnorm(x,mu.0,sig/sqrt(n.0))}
posterior=function(x){dnorm(x,mu.p,sig.p)}
lik=function(x){dnorm(x,x.med,sig/sqrt(n))}

# (e)
curve(prior(x),from=40,to=200,ylim=c(0,0.15))
curve(lik,add=T,lty=2)
curve(posterior,add=T,lty=3)

# (f)
# stima puntuale
mu.p
# stima intervallare ET=HPD
L=qnorm(0.025,mu.p, sig.p)
U=qnorm(0.975,mu.p, sig.p)
c(L,U)

# (g)
1-pnorm(135,mu.p,sig.p)

# (h)
smv=x.med
smv
IC=c(x.med-qnorm(0.975)*sig/sqrt(n), x.med+qnorm(0.975)*sig/sqrt(n))
IC
# p.value per theta <=135 vs. theta>135
W=sqrt(n)*(x.med-135)/sig
p.value=1-pnorm(W)
p.value


# Esercizio 2 (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)