# Impostazione del seme
set.seed(345857)

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