# Ex. 1

M = 10^4

delta = 0.3

lambda = 0.25

omega = rbeta(M, 2, 3)

W = 2 * abs(delta-omega)

mean(W>lambda)

# 0.5762

# Ex. 2

rm(list = ls())

n = 20

M = 10^4

x = log( rgamma(M*n, shape = 3, rate = 4) ) 

x = matrix(x, M, n)

x = rowSums(x)/n

mean(x)

# -0.4645

# Possibile confrontare con momento teorico; E[ln X] = digamma(3) - log(4) = -0.4635

# Ex. 3

rm(list = ls())

M = 10^4

n = 20

sn = 10

alpha.post = 4 + sn

beta.post = 6 + n - sn

theta = rbeta(M, alpha.post, beta.post)

Y = -theta*log(theta) - (1-theta)*log(1-theta)

mean(Y)

# 0.67

# Ex. 4

# Per risolvere l'esercizio, le strade sono molteplici. Vediamone due.

# Primo approccio.

rm(list = ls())

n = 10

M = 10^4

theta = 5 # rate = 1/theta o, equivalentemente, scale = theta 

x = rexp(M*n, rate=1/theta)

x = matrix(x, M, n)

x = rowSums(x)/n

L = x/1.42

U = x/0.62

mean(L<theta & U>theta)

# 0.8 circa

# Secondo approccio, basato su distribuzione della media.

rm(list = ls())

n = 10

M = 10^4

theta = 5 

medie = rgamma(M, shape = n, rate = n/theta)

L = medie/1.42

U = medie/0.62

mean(L<theta & U>theta)


# 0. 7965


# Ex. 5 (si sfruttano le medie precedentemente generate per confronto)

L1 = medie - 1.281552 * ( medie/sqrt(n) )

U1 = medie + 1.281552 * ( medie/sqrt(n) )

len = U-L

len1 = U1-L1

mean(len) # 4.53

mean(len1) # 4.05


# Ex. 6

rm(list = ls())

M = 10^4

n = 8

# Il rischio è pari alla varianza theta/n perché L = (theta-d)^2 e E[d] = theta. 

# Generiamo i theta a priori per simulare rischio di Bayes.

theta = rgamma(M, shape = 7, rate = 5)

r = theta/n

mean(r)

# 0.1751


# Ex. 7

rm(list = ls())

theta = 4

mi = 4

n = 25

k = 40

M = 10^4

x = rnorm(M*n, mi, sqrt(theta))

x = matrix(x, M, n)

medie = rowSums(x)/n

scarti = matrix(M*n, M, n)

for(camp in 1:M){
  for(i in 1:n){
    
    scarti[camp, i] = x[camp,i] - medie[camp]
    
  }
}

scarti = scarti^4

mom.camp = rowSums(scarti)/n

mean(mom.camp<k)

# 0.5438


# Ex. 8 

# a)

# Anche in questo caso si può sfruttare
# l'uguaglianza in distribuzione tra somma di Esponenziali e Gamma e procedere in due modi.

# Primo approccio.

rm(list = ls())

M = 10^4

n = 30

theta0 = 3

k = 7.2

x = rexp(M*n, theta0)

x = matrix(x, M, n)

somme = rowSums(x)

mean(somme<k)

# circa 0.05


# Secondo approccio:

rm(list = ls())

M = 10^4 # M = numero di campioni poiché simuliamo lo spazio campionario (approccio frequentista)

n = 30

theta0 = 3

k = 7.2

somme = rgamma(M, shape = n, rate = theta0) 

mean(somme<k)

# circa 0.05


# b), sfruttando secondo approccio

theta1 = 4 # sotto hp. che theta = theta1 = 4

somme = rgamma(M, shape = n, rate = theta1) 

mean(somme<k) # P(Sn < k | theta = 4) cioè la f.ne di potenza in theta = 4

# 0.4363


# Ex.9 (facoltativo)


rm(list = ls())

# a)

optimal = function(n){
  
  set.seed(123)
  
  M = 10^4
  
  gamma = 0.05
  
  theta = 5
  
  z = qnorm(1-gamma/2)
  
  x = rexp(M*n, rate = 1/theta )
  
  x = matrix(x, M, n)
  
  medie = rowSums(x)/n
  
  L = medie - z*( medie/sqrt(n) )
  
  U = medie + z*( medie/sqrt(n) )
  
  copertura = mean(L<theta & U>theta)
  
  return(copertura)
  
}

# b)

optimal = Vectorize(optimal)

coperture = optimal(1:100) # n varia da 1 a 100

min ( which( coperture > 0.85 ) ) 

# n* = 4