# Impostazione del seme
set.seed(345857)


# ESERCIZIO 1

# Parametri 
a <- runif(1, min = 2.5, max = 5.5)  # parametro 'a' estratto da U(2.5, 5.5)
b <- 1                                # parametro 'b'
p <- 0.5                              # probabilità per Bernoulli
lambda <- runif(1, min = 2, max = 4)  # parametro 'lambda' estratto da U(2, 4)
xseq = seq(0, 15, by=0.001)
n_samples <- 10000  # numero di campioni Monte Carlo

# 1. Generazione campioni marginali di Y e stima Monte Carlo della CDF

# Generazione campioni per Z da una distribuzione Bernoulli con probabilità p
Z <- rbinom(n_samples, 1, p)  # Z ~ Bern(p)

# Generazione di Y condizionato al valore di Z
Y <- ifelse(Z == 0, rpois(n_samples,lambda), rgamma(n_samples, a, b))

cum_Y = c()
for (i in 1:length(xseq)) {
  cum_Y[i] = sum(Y<=xseq[i])/n_samples
}

plot(xseq, cum_Y, type='l', xlab = "Y", ylab = "F_Y(y)")






# 2. Stima con Monte Carlo di P(Y ∈ [3.5, 4.5]) e P(Y ∈ [3.5, 4.5] | Z = 0)
# e stima della varianza degli stimatori

# Calcoliamo la probabilità P(Y ∈ [3.5, 4.5]) come proporzione di campioni
p_Y_interval <- mean(Y >= 3.5 & Y <= 4.5)

# Condizionatamente a Z = 0, stima di P(Y ∈ [3.5, 4.5] | Z = 0)
Y_Z0 <- Y[Z == 0]
p_Y_interval_Z0 <- mean(Y_Z0 >= 3.5 & Y_Z0 <= 4.5)

# Varianza degli stimatori
var_p_Y_interval <- p_Y_interval + (a-p_Y_interval) / n_samples
var_p_Y_interval_Z0 <- p_Y_interval_Z0 * (p_Y_interval_Z0) / length(Y_Z0)

# Output dei risultati
cat("P(Y ∈ [3.5, 4.5]):", p_Y_interval, "\n")
cat("Varianza stimatore di P(Y ∈ [3.5, 4.5]):", var_p_Y_interval, "\n")
cat("P(Y ∈ [3.5, 4.5] | Z = 0):", p_Y_interval_Z0, "\n")
cat("Varianza stimatore di P(Y ∈ [3.5, 4.5] | Z = 0):", var_p_Y_interval_Z0, "\n")

# 3. Stima con Monte Carlo di f(Z = 0 | Y ∈ [1.5, 3.5]) e f(Z = 1 | Y ∈ [1.5, 3.5])
#Y_interval <- Y >= 1.5 & Y <= 3.5
#f_Z0_given_Y_interval <- mean(Z[Y_interval] == 0)
#f_Z1_given_Y_interval <- mean(Z[Y_interval] == 1)
Y_Z0 <- Y[Z == 0]
Y_Z1 <- Y[Z == 1]
p_Y_interval_Z0 <- mean(Y_Z0 >= 1.5 & Y_Z0 <= 3.5)
p_Y_interval_Z1 <- mean(Y_Z1 >= 1.5 & Y_Z1 <= 3.5)
p_Y_interval <- mean(Y >= 1.5 & Y <= 3.5)
f_Z0_given_Y_interval <- p_Y_interval_Z0*(1-p)/p_Y_interval
f_Z1_given_Y_interval <- p_Y_interval_Z1*p/p_Y_interval


# Output dei risultati
cat("f(Z = 0 | Y ∈ [1.5, 3.5]):", f_Z0_given_Y_interval, "\n")
cat("f(Z = 1 | Y ∈ [1.5, 3.5]):", f_Z1_given_Y_interval, "\n")

# 4. Stima dei quantili empirici per le probabilità date (0.1, 0.2, 0.5, 0.75)
quantiles_Y <- quantile(Y, probs = c(0.1, 0.2, 0.5, 0.75))
quantiles_Z <- quantile(Z, probs = c(0.1, 0.2, 0.5, 0.75))

# Output dei quantili
cat("Quantili empirici di Y:", quantiles_Y, "\n")
cat("Quantili empirici di Z:", quantiles_Z, "\n")

##########################################################################


# ESERCIZIO 2

# Parametri dati dall'esercizio
set.seed(345857)  # Imposta un seed per riproducibilità
mu <- runif(1, -1.5, 1.5)      # Simula mu da U(-1.5, 1.5)
sigma2 <- runif(1, 0.5, 1.5)   # Simula sigma^2 da U(0.5, 1.5)
sigma <- sqrt(sigma2)           # Deviazione standard

# Funzione di densità logistic-normal
logistic_normal_density <- function(x, mu, sigma2) {
  (1 / (x * (1 - x) * sqrt(2 * pi * sigma2))) * 
    exp(-0.5 * ((log(x / (1 - x)) - mu)^2 / sigma2))
}


# 1. campioni U e Y

n = 10 # numero di campioni
m = 1.75 # plottando la funzione, si vede che il valore massimo è circa 1.75
Y = runif (n ,0 ,1)
U = runif (n ,0 , m)
X = Y [U < logistic_normal_density(Y, mu, sigma2)]
U_X = U[U < logistic_normal_density(Y, mu, sigma2)]

# plot densità
xseq = seq (0 ,1 , by =0.01)
plot (xseq , logistic_normal_density(xseq, mu, sigma2) , ylim =c (0 , m) ,type ="l" , lwd =2 ,)
points (Y ,U , pch =20 , cex = 0.1)
points (X ,U_X , pch =20 , cex = 0.1 , col =2)
lines ( density (X , from =0 , to = 1) , col =2 , lwd =2)












# 2. prior per mu: N(0, 100)
mean_mu <- 0
var_mu <- sigma2
prior_mu <- dnorm(xseq, mean_mu, sqrt(var_mu))

# Calcola media e varianza posteriori
posterior_sigma2 <- 1 / (1 / var_mu + n / sigma2)
posterior_mu <- posterior_sigma2 * (prior_mu * n * mean(log(X / (1 - X))) / sigma2  )

# Simula campioni dalla distribuzione a posteriori di mu
posterior_samples_mu <- rnorm(1000, posterior_mu, sqrt(posterior_sigma2))



## Visualizza la distribuzione a posteriori di mu
#hist(posterior_samples_mu, breaks = 30, main = "Distribuzione a posteriori di mu", xlab = "mu")


# 3.

n2 <- 100
Y2 = runif (n2 ,0 ,1)
U2 = runif (n2 ,0 , m)
X2 = Y2 [U2 < logistic_normal_density(Y2, mu, sigma2)]
U_X2 = U2[U2 < logistic_normal_density(Y2, mu, sigma2)]

# Calcola media e varianza posteriori
posterior_sigma2_2 <- 1 / (1 / var_mu + n2 / sigma2)
posterior_mu_2 <- posterior_sigma2 * (prior_mu * n2 * mean(log(X2 / (1 - X2))) / sigma2  )

# Simula campioni dalla distribuzione a posteriori di mu
posterior_samples_mu_2 <- rnorm(1000, posterior_mu_2, sqrt(posterior_sigma2_2))

# confronto grafico
pri_dens <- density(rnorm(xseq, mean_mu, sqrt(var_mu)))
post_dens <- density(posterior_samples_mu)
post_dens_2 <- density(posterior_samples_mu_2)

plot(pri_dens, type = 'l', main = '', xlim= c(-2,2), ylim = c(0, 2.5), col=4, lwd = 2)
lines(post_dens, col=2, lwd = 2)
lines(post_dens_2, col=3, lwd = 2)

legend("topright", legend = c('Prior', 'Posterior_n_10', 'Posterior_n_100', col = 'blue', 'red', 'green'))





# 4. Importance sampling per stimare E(X)
alpha <- 1.2
beta <- 1.2
g <- function(x) dbeta(x, alpha, beta)  # Funzione di densità della Beta(1.2, 1.2)

# Genera campioni dalla distribuzione di importanza
importance_samples <- rbeta(1000, alpha, beta)

# Stima di E(X) usando importance sampling
E_X <- mean(importance_samples * logistic_normal_density(importance_samples, mu, sigma) / g(importance_samples))