# Lab 09 TSdD 2024 in classe

# Es. 1 (AfE con perdita assoluta)
# (a) funzione R per perdita finale attesa (rho)
al.p=5
be.p=8
rho.abs.fun=function(a){
  integranda.fun=function(theta){abs(a-theta)*dbeta(theta,al.p,be.p)}
  return(integrate(integranda.fun,0,1)$value)
}
# valore di rho per due possibili stime (azioni)
rho.abs.fun(0.3)
rho.abs.fun(0.4)
# (b) grafico di rho al variare di a
a.val=seq(0,1,0.001)
length(a.val)
a.val=matrix(a.val,length(a.val),1)
rho.ABS.val=apply(a.val,1,rho.abs.fun)
plot(a.val,rho.ABS.val,type="l",xlab="a",ylab = "perdita finale attesa",ylim=c(0,0.5))
# (c) Calcolo di a.star
# Ricorda: a.star=mediana di pi(theta|zn)
# da risultati analitici (teoria) so che
a.ABS.star=qbeta(0.5,al.p,be.p)
a.ABS.star
abline(v=a.ABS.star)
# numericamente (con R)
a.val[rho.ABS.val==min(rho.ABS.val)]
# perdita finale minima
rho.abs.fun(a.star)
# per confrontare perdita finale attesa di azione non ottima con rho di azione ottima
rho.abs.fun(0.5)/rho.abs.fun(a.ABS.star)


# Es. 2 (AfE con perdita lineare asimmetrica)
# (a)  funzione R per perdita finale attesa (rho)
b=1
c=2
rho.LA.fun=function(a){
  integranda.fun=function(theta){(b*(theta-a)*(theta>=a)+c*(a-theta)*(theta<a))*dbeta(theta,al.p,be.p)}
  return(integrate(integranda.fun,0,1)$value)
}
# ad esempio, in a=0.3
rho.LA.fun(0.3)
# (b) grafico di rho al variare di a
rho.LA.val=apply(a.val,1,rho.LA.fun)
plot(a.val,rho.LA.val,type="l",xlab="a",ylab = "perdita finale attesa",ylim=c(0,2))
lines(a.val,rho.ABS.val,lty=2,lwd=3)
# (c) calcolo azione ottima a.star rispetto perdita LA
# Ricorda: a.star quantile di livello b/(b+c) di pi(theta|zn)
# analiticamente
a.LA.star=qbeta(b/(b+c),al.p,be.p)
a.LA.star
# numericamente
a.val[rho.LA.val==min(rho.LA.val)]
# (d) calcolo di  rho(a.LA.star)
rho.LA.fun(a.LA.star)
# (e) ripeti con b=c=1


# Es. 3 calcolo rischio di Bayes
al=3
be=6
n=10
M=10000
th.MC=rbeta(M,al,be)
r.d.val=th.MC*(1-th.MC)/n
r.d.MC=mean(r.d.val)
r.d.MC