# LDS - Prova finale 1 - dicembre 2024 - Compito B

# Es. 1
rm(list = ls())
delta = 3
lambda = 1
M = 50000
mi = 4
sd = sqrt(3)
omega.MC = rnorm(M, mi, sd)
loss.MC = abs(delta - omega.MC)
mean(loss.MC > lambda) 
# 0.625


# Es. 2
rm(list = ls())
M = 50000 
n = 20
mi = 5
sd = 2
x.MC = rnorm(M*n, mi, sd)^3
x.MC = matrix(x.MC, M, n)
sample.means = apply(x.MC, 1, mean)
mean(sample.means) # Simulazione del valore atteso sullo spazio campionario
# 185


# Es. 3
rm(list = ls())
M = 50000
alp.prior = 4
beta.prior = 6
n = 18
sn = 8
alp.post = alp.prior + n
beta.post = beta.prior + sn
theta.MC = rgamma(M, alp.post, rate = beta.post )
theta.MC = theta.MC^(-3)
mean(theta.MC) # 0.3433449
# Verifica con momento teorico:
# igamma(alp.post-3,0)/igamma(alp.post,0)/beta.post^(-3)
#  0.344


# Es. 4
rm(list = ls())
theta.vero = 5
n = 10
M = 50000
x.MC = rexp(n*M, theta.vero)
x.samples = matrix(x.MC, M, n)
bounds = matrix(NA, M, 2)
copertura = 0 
# Se intervallo contiene theta, copertura = copertura + 1
for(camp in 1:M){
  sample.mean = mean(x.samples[camp,])
  L = 1/sample.mean - 1/ (sqrt(n)*sample.mean)
  U = 1/sample.mean + 1/ (sqrt(n)*sample.mean)
  bounds[camp, 1:2] = c(L,U)
  if(L<theta.vero & U > theta.vero){
    copertura = copertura + 1
  }
}
copertura/M 
# Su M campioni, quanti generano intevalli che contengono theta?
# 0.69
# In alternativa
rm(list = ls())
theta.vero = 5
n = 10
M = 50000
x.matr=rgamma(M*n,1,rate=theta.vero)
x.matr=matrix(x.matr,M,n)
sn.mc=apply(x.matr,1,sum)
L=n/sn.mc - n/(sn.mc*sqrt(n))
U=n/sn.mc + n/(sn.mc*sqrt(n))
# copertura
mean(L<theta.vero&U>theta.vero)
# 0.69


# Es. 5
rm(list = ls())
n = 15
sn = 6
alp.prior = 4
beta.prior = 5
alp.post = alp.prior + sn
beta.post = beta.prior + n
gg=1-0.8
L = qgamma(gg/2, shape = alp.post, rate = beta.post)
U = qgamma(1-gg/2, shape = alp.post, rate = beta.post)
c(exp(-U), exp(-L))
# 0.49, 0.73


# Es. 6
rm(list = ls())
n = 8
M = 50000
theta.MC = rbeta(M, 7, 5)
var.MC = (theta.MC)*(1-theta.MC)/n
r = mean(var.MC)
# 0.028


# Es. 7
rm(list = ls())
k = 2
mi = 1
theta.vero = 5
sd = sqrt(theta.vero)
M = 50000
n = 15
x.MC = rnorm(M*n, mi, sd )
x.MC = matrix(x.MC, M, n)
scarti = matrix(NA, M, n)
for(camp in 1:M){
  media = mean( x.MC[camp,] )
  scarti[camp,] = abs( x.MC[camp,] - media )
}
scarti = apply(scarti, 1, mean)
mean(scarti < k)
# 0.794


# Es. 8
rm(list = ls())
M = 50000
theta0 = 3
theta1 = 2
theta.vero = 2.5
x.sim = rexp(M, theta.vero) 
# la somma è, per n=1, il valore della singola Xi
lambda = (theta0/theta1)*exp(-x.sim)
mean(lambda)
# 1.07