# 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, 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))