Ejercicio 1

#2 mezcla de gammas.
# Y = .2*Gamma((1,10,1/30)) + .8*Gamma((1,50,1/20))
# Parámetros
m = 1000
x = numeric(m)  # Inicializa el vector x
b = numeric(m)  # Inicializa el vector b
a1 = 0.20
a2 = 0.80
b = rbinom(m, 1, a2)  # Asignación de mezcla

# Simulación de la mezcla de distribuciones gamma
for (i in 1:m) {
  if (b[i] == 0) {
    x[i] = rgamma(1, shape = 20, rate = 7)  # Gamma(20, 7)
  } else {
    x[i] = rgamma(1, shape = 50, rate = 7)  # Gamma(50, 7)
  }
}

# Histograma de la mezcla de gammas
x_hist = hist(x, probability = TRUE, breaks = 40)

# Graficar histograma
hist(x, probability = TRUE, breaks = 40, 
     main = expression(bold("Gamma(20,7) y Gamma(50,7)")), 
     xlab = "Valores", 
     col = "lightblue", 
     border = "white",  # Color del histograma y borde
     xlim = c(0, 12.5), # Ajuste para el eje x
     ylim = c(0, 0.7))  # Ajuste para el eje y

# Añadir las curvas teóricas
curve(dgamma(x, shape = 20, rate = 7), from = 0, to = 12.5, 
      add = TRUE, col = "red", lwd = 2)  # Curva Gamma(20, 7)
curve(dgamma(x, shape = 50, rate = 7), from = 0, to = 12.5, 
      add = TRUE, col = "blue", lwd = 2)  # Curva Gamma(50, 7)

# Etiquetas para las curvas
text(5, 0.5, "Gamma(20, 7)", col = "red", cex = 1.2)  # Etiqueta para Gamma(20, 7)
text(5, 0.4, paste("Media:", round(20/7, 3)), col = "red", cex = 0.8)

text(10, 0.3, "Gamma(50, 7)", col = "blue", cex = 1.2)  # Etiqueta para Gamma(50, 7)
text(10, 0.2,paste("Media:", round(50/7, 3)), col = "blue", cex = .8)

Imprimimos los valores de nuestra función de densidad.

# Mostrar las probabilidades del histograma
breaks = x_hist$breaks
density = x_hist$density

# Mostrar los breaks y densidad calculados
print(breaks)

 [1]  1.2  1.4  1.6  1.8  2.0  2.2  2.4  2.6  2.8  3.0  3.2  3.4  3.6  3.8  4.0
[16]  4.2  4.4  4.6  4.8  5.0  5.2  5.4  5.6  5.8  6.0  6.2  6.4  6.6  6.8  7.0
[31]  7.2  7.4  7.6  7.8  8.0  8.2  8.4  8.6  8.8  9.0  9.2  9.4  9.6  9.8 10.0
[46] 10.2 10.4 10.6

print(density)

 [1] 0.005 0.000 0.040 0.055 0.095 0.095 0.090 0.140 0.150 0.045 0.065 0.045
[13] 0.070 0.035 0.025 0.005 0.010 0.015 0.010 0.035 0.030 0.110 0.120 0.115
[25] 0.255 0.270 0.330 0.280 0.285 0.355 0.325 0.270 0.275 0.205 0.185 0.100
[37] 0.125 0.145 0.065 0.035 0.040 0.015 0.020 0.005 0.000 0.000 0.010

Ahora calcularemos las probabilidades acumuladas. Usamos la siguiente formula.

# Calcular P(Y < 200)
p1 = pgamma(200, shape = 20, rate = 7)  # P(Y < 200 | Gamma(20, 7))
p2 = pgamma(200, shape = 50, rate = 7)  # P(Y < 200 | Gamma(50, 7))

# Probabilidad total P(Y < 200)
probabilidad = a1 * p1 + a2 * p2

# Imprimir resultado
print(paste("P(Y < 200) =", round(probabilidad, 4)))

[1] "P(Y < 200) = 1"

De modo que tenemos para P<400 y P<500 los siguientes datos.

# Calcular P(Y < 400)
p400_1 = pgamma(400, shape = 20, rate = 7)  
p400_2 = pgamma(400, shape = 50, rate = 7)  

# Probabilidad total P(Y < 400)
probabilidad2 = a1 * p400_1 + a2 * p400_2

# Calcular P(Y < 500)
p500_1 = pgamma(500, shape = 20, rate = 7)  
p500_2 = pgamma(500, shape = 50, rate = 7)  

# Probabilidad total P(Y < 500)
probabilidad3 = a1 * p500_1 + a2 * p500_2

# Imprimir resultado
print(paste("P(Y < 400) =", round(probabilidad2, 4)))

[1] "P(Y < 400) = 1"

print(paste("P(Y < 500) =", round(probabilidad3, 4)))

[1] "P(Y < 500) = 1"

Ahora para la esperanza tenemos:

Por lo tanto

# Esperanzas de cada distribución
E_gamma_20_7 = 20 / 7  # Esperanza de Gamma(20, 7)
E_gamma_50_7 = 50 / 7  # Esperanza de Gamma(50, 7)

# Esperanza de la mezcla
E_Y = a1 * E_gamma_20_7 + a2 * E_gamma_50_7

# Imprimir resultado
print(paste("E[Y] =", round(E_Y, 4)))

[1] "E[Y] = 6.2857"

Ahora para el último caso tenemos lo siguiente

# Calcular P(Y < E[Y])
P_Y_less_than_EY = a1 * pgamma(E_gamma_20_7, shape = 20, rate = 7) +
                    a2 * pgamma(E_gamma_50_7, shape = 50, rate = 7)

# Imprimir resultado
print(paste("P(Y < E[Y]) =", round(P_Y_less_than_EY, 4)))

[1] "P(Y < E[Y]) = 0.521"