# Lab 6 2024

# 1. 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)



# Predictive inference MC

# Problem 2
m=8
n.0=5
mu.0=2
sig2=1
M=10000
sd.prior=sqrt(sig2/n.0)

th.tilde=rnorm(M,mu.0,sd.prior)
x.pred.MC=rnorm(M,th.tilde,sqrt(sig2/m))

length(x.pred.MC)

# confronto istogramma e densità esatta
hist(x.pred.MC,freq=F,ylim=c(0,0.8))

sd.pred.prior=sqrt(sig2*(1/n.0+1/m))
curve(dnorm(x,mu.0,sd.pred.prior),add=T)

# confronto valori esatti e approx MC (valore atteso e varianza)

mean.pred.MC=mean(x.pred.MC)
mu.0
mean.pred.MC

sd.pred.prior.MC=sqrt(var(x.pred.MC))
sd.pred.prior.MC
sd.pred.prior

# potrei anche generare direttamente da N(mu.0,var.pred.prior)

x.pred.MC.2=rnorm(M,mu.0,sd.pred.prior)
length(x.pred.MC.2)
hist(x.pred.MC.2,freq=F,breaks=20,add=T)
curve(dnorm(x,mu.0,sd.pred.prior),add=T)
mean(x.pred.MC.2)

sd.pred.prior.MC.2=sqrt(var(x.pred.MC.2))
sd.pred.prior.MC.2
sd.pred.prior

# Problem 3

m=8
n.0=5
mu.0=2
sig2=1
n=6
x.med=0.5
M=10000

sd.post=sqrt(sig2/(n.0+n))
var.post=sig2/(n.0+n)
mean.post=(n.0*mu.0+n*x.med)/(n.0+n)
sd.pred.post=sqrt(sig2/m+var.post)

th.tilde=rnorm(M,mean.post,sd.post)
x.pred.MC=rnorm(M,th.tilde,sqrt(sig2/m))

mean.post.MC=mean(x.pred.MC)
sd.pred.post.MC=sqrt(var(x.pred.MC))

mean.post.MC
mean.post

sd.pred.post.MC
sd.pred.post

hist(x.pred.MC,freq=F,ylim=c(0,0.9))
curve(dnorm(x,mean.post,sd.pred.post),add=T)


# Problem 4

al=1/2
be=1/2
th.vero=1/4
n=5
M=5000

x.matrix=matrix(rbinom(M*n,size=1,prob=th.vero),M,n)
sn.MC=apply(x.matrix,1,sum)
al.post=al+sn.MC
be.post=be+n-sn.MC

ln.MC=qbeta(0.025,al.post,be.post)
un.MC=qbeta(0.975,al.post,be.post)

mean(ln.MC<th.vero&un.MC>th.vero)


# Problema 5 (Prob successo esperimento)

n0=5 
mu0=2 
m=6 
sig2=1 
th0=3 
alpha=0.05 
n=10 
xmed=4

# Risultati esatti
mu.prior.pred=mu0 
sig.prior.pred=sqrt(sig2*(1/m+1/n0)) 
mu.post.pred=(n0*mu0+n*xmed)/(n0+m) 
sig.post.pred=sqrt(sig2*(1/m+1/(n0+n))) 
alpha=0.05
z.al=qnorm(1-alpha/2)
z.al

# success if |Ym-theta0|>sig*z/sqrt(m)+th0
# ovvero 
# P(Ym<H1)+P(Ym>H2)  con

H2=sqrt(sig2/m)*z.al+th0
H1= - sqrt(sig2/m)*z.al+th0

# Prior prob of success 
pnorm(H1,mu.prior.pred,sig.prior.pred)+(1-pnorm(H2,mu.prior.pred,sig.prior.pred))

# Post prob prob of success 
pnorm(H1,mu.post.pred,sig.post.pred)+(1-pnorm(H2,mu.post.pred,sig.post.pred))


# Per approssimazioni MC, 
# generare da predittiva a priori e poi predittiva a posteriori di Ym e verificare

# Ad esempio, predittiva a priori

# (a) con simulazione da predittiva a priori
M=10000
th.MC=rnorm(M,mu0, sqrt(sig2/n0))
y.MC=rnorm(M,th.MC,sqrt(sig2/m))
tmp=sqrt(m)*abs(y.MC-th0)/sqrt(sig2)

mean(tmp>z.al)