Taller método aceptacion y rechazo

Universidad Nacional de Colombia

Asignatura: Computación Estadística.

Fecha: 7 de septiembre de 2026

Solución

1). Para la función objetivo \(f(x)= 2x\) en el intervalo \([0,1]\) ¿cuál es el valor máximo que alcanza la función?

como \(x \in [0,1]\)

\(\Rightarrow\) el valor más grande en el intervalo es 1

\[f(x)=2(1)=2\]

\[\therefore f(x)=2\]

R/T: el valor maximo es 2

2). Si usamos una propuesta uniforme. \(g(x)=1\) ¿cuál debe ser el valor exacto de la constante \(c\) para que se cumpla \(f(x) \leq c \cdot g(x)\) en todo el dominio?.

sabemos que el valor maximo de \(2x\) en el intervalo \([0,1]\) es 2

\(\Rightarrow\) tenemos la función \(f(x) \leq c \cdot g(x)\)

y sabemos que \[f(x)=2x\] \[g(x)=1\]

\[2x \leq2\]

\[\Rightarrow 2\leq 1(2) = 2\]

\(\Rightarrow\) la constante tiene un valor exacto de 2

3). ¿Qué pasaría con la eficiencia del algoritmo si por error definimos \(c=20\) en lugar de \(c=2\)?

con \(c=2\)

la probabilidad de aceptación es de: \[\frac{1}{2}= 0.5(100) =50\%\] pero con \(c=20\) la probabilidad de aceptación es:

\[\frac{1}{20}= 0.05(100) =5\%\]

R/T: El algoritmo seguiría funcionando y los números aceptados seguirían la distribución correcta, pero sería mucho menos eficiente, Esto significa que se rechazarían muchísimos más números y el algoritmo necesitaría generar muchas más propuestas.

¿Fallarán los números generados o qué métrica se verá afectada?.

R/T: Los números generados no necesariamente estarán mal; la distribución de los números aceptados seguirá siendo correcta, lo que se verá afectado será principalmente es la tasa de aceptación y eficiencia del algoritmo.

Vamos a simular una variable aleatoria cuya función de densidad objetivo es \(f(x)=3x^{2}\) para \(x \in [0,1]\). Como distribución propuesta \(g(x)\) usaremos una Uniforme en \([0,1]\), donde \(g(x)=1\). El valor máximo de \(f(x)\) en ese intervalo es \(3\), por lo que nuestra constante será \(c=3\). La condición de aceptación se reduce a:

\[\frac{3Y^{2}}{3\cdot1}=Y^{2}\]

# Parámetros de la simulación
n_sim <- 10000  # Número de valores que queremos generar
muestras <- numeric(n_sim)
intentos <- 0
aceptados <- 0

set.seed(123) # Para reproducibilidad

while (aceptados < n_sim) {
  intentos <- intentos + 1
  
  # 1. Generar candidato Y de la propuesta Uniform(0,1)
  Y <- runif(1, min = 0, max = 1)
  
  # 2. Generar U Uniform(0,1)
  U <- runif(1, min = 0, max = 1)
  
  # 3. Prueba de aceptación (U <= f(Y) / (c * g(Y))) -> U <= Y^2 / 1
  if (U <= (Y^2)) {  # Notar que f(Y)/(c*g(Y)) = (3*Y^2)/(3*1) = Y^2
    aceptados <- aceptados + 1
    muestras[aceptados] <- Y
  }
}

# Calcular la tasa de eficiencia
eficiencia <- n_sim / intentos
cat("Total de intentos:", intentos, "\n")
## Total de intentos: 30226
cat("Eficiencia teórica vs real:", round(eficiencia * 100, 2), "%\n")
## Eficiencia teórica vs real: 33.08 %

Ejercicio Práctico

Cambia la distribución objetivo: Simula una variable aleatoria con densidad triangular \(f(x)=2x\) para \(x∈[0,1]\). Define los parámetros: Utiliza como propuesta una distribución \(g(x)=1\) (Uniforme) y determina el valor adecuado de la constante \(c\). Ajusta el código: Escribe la condición de aceptación correspondiente y genera una muestra de \(5000\) valores.Comprueba: Grafica el histograma resultante y superpone la línea de la función teórica \(2x\) usando la función curve(). Muestra la tasa de eficiencia del algoritmo en la consola.

# Parámetros de la simulación

n_sim <- 5000
muestras <- numeric(n_sim)

intentos <- 0
aceptados <- 0

set.seed(123)
# Función objetivo

f <- function(x) {
  2 * x
}

# Función propuesta
g <- function(x) {
  1
}

c <- 2


#  Generar muestras
while (aceptados < n_sim) {
  
  # Contar intento
  intentos <- intentos + 1
  
  # Generar candidato Y
  Y <- runif(1, min = 0, max = 1)
  
  # Generar número U
  U <- runif(1, min = 0, max = 1)
  
  # Condición de aceptación
  # U <= f(Y) / (c * g(Y))
  # U <= (2Y)/(2*1)
  # U <= Y
  
  if (U <= Y) {
    
    aceptados <- aceptados + 1
    muestras[aceptados] <- Y
    
  }
  
}


eficiencia <- n_sim / intentos
cat("Total de intentos:", intentos, "\n")
## Total de intentos: 9976
cat("Eficiencia:", 
    round(eficiencia * 100, 2), "%\n")
## Eficiencia: 50.12 %
hist(muestras,
     probability = TRUE,
     breaks = 40,
     col = "lightblue",
     main = "Muestreo por Aceptación y Rechazo",
     xlab = "X",
     ylab = "Densidad")

curve(2 * x,
      add = TRUE,
      col = "red",
      lwd = 2)

Reto de Experimentación

Modifica el parámetro de la constante en tu script:

1). Cambia temporalmente tu constante \(c\) a un valor de 15(manteniendo \(f(x)=2x\)). ¿Cuántos intentos necesitó ahora para llegar a las 5000 muestras válidas?

obs:

Total de intenos: Es la cantidad total de candidatos generados que el algoritmo tuvo que probar para obtener la muestra deseada.

Eficiencia: Es el porcentaje de esos intentos fue aceptado.(si la eficiencia es alta, significa que el algoritmo necesitó rechazar pocos candidatos y si es baja, tuvo que hacer muchos intentos para conseguir los valores aceptados.)

# Parámetros de la simulación
n_sim <- 5000
muestras <- numeric(n_sim)
intentos <- 0
aceptados <- 0

set.seed(123)
# Función objetivo

f <- function(x) {
  2 * x
}

# Función propuesta
g <- function(x) {
  1
}

c <- 15


#  Generar muestras
while (aceptados < n_sim) {
  
  # Contar intento
  intentos <- intentos + 1
  
  # Generar candidato Y
  Y <- runif(1, min = 0, max = 1)
  
  # Generar número U
  U <- runif(1, min = 0, max = 1)
  
  # Condición de aceptación
  # U <= f(Y) / (c * g(Y))
  # U <= (2Y)/(2*1)
  # U <= Y
  
  if  (U <= (2 * Y) / 15) {
    
    aceptados <- aceptados + 1
    muestras[aceptados] <- Y
    
  }
  
}

eficiencia <- n_sim / intentos
cat("Total de intentos:", intentos, "\n")
## Total de intentos: 74581

¿Qué porcentaje de eficiencia muestra la consola comparado con usar \(c=2\)?

cat("Eficiencia:", 
    round(eficiencia * 100, 2), "%\n")
## Eficiencia: 6.7 %
hist(muestras,
     probability = TRUE,
     breaks = 40,
     col = "lightblue",
     main = "Muestreo por Aceptación y Rechazo",
     xlab = "X",
     ylab = "Densidad")

curve(2 * x,
      add = TRUE,
      col = "red",
      lwd = 2)

¿Por qué ocurre esto geométricamente?

Porque el área que usamos para generar candidatos se vuelve mucho más grande, pero la función objetivo sigue siendo la misma. Entonces, una gran parte de los puntos generados queda por encima de la curva objetivo, por lo que son rechazados.

Auditoría del Modelo

1. Si el algoritmo rechaza tantos candidatos al inicio del ciclo while, ¿por qué el histograma final de los valores aceptados se ajusta tan perfectamente a la curva teórica \(2x\)?

R/T: Porque, aunque el algoritmo genere muchos valores y rechace algunos, la regla de aceptación está construida para que los valores que sí se guardan tengan la distribución \(f(x)=2x\), dicho de en otras palabras, el método genera muchos valores aleatorios entre 0 y 1, pero solo guarda aquellos que cumplen la condición de aceptación.

2. ¿Qué relación guarda el área bajo la curva roja con la probabilidad de que un candidato sea aceptado en el primer intento?

R/T: \[f(x)=2x, \ \ 0\leq x\leq 1\] entonces el area total bajo la curba es: \[ \int_{0}^{1}x^{2} dx\]

\[=[x^{2}]_{0}^{1}=1\]

Eso es una dencidad Válida por tanto probabilidad de aceptación depende del valor de la constante, entonces cuanto más grande sea el espacio que se asuma para realizar la aceptación y rechazo, menor será la proporción de valores que quedan debajo de la curva, por tanto, menor será la probabilidad de aceptación.

3. Si cambiáramos el dominio de la función a \([0,2]\), ¿cómo cambiaría la integral de la densidad y qué ajustes obligatorios tendrías que hacer en la función runif() de los candidatos?

R/T:
\[ \int_{0}^{1}x^{2} dx\]

\[=[x^{2}]_{0}^{2}=4\]

ya no sería una densidad de probabilidad válida, porque una densidad debe cumplir que la integral de la funcion sea igula a 1.