\[Muestreo \ por \ Rechazo \ Adaptativo\]

Fundamento conceptual

pregunta 1 Explique con sus propias palabras qué problema resuelve el método de rejection sampling clásico y por qué puede ser ineficiente cuando la constante M es grande. Relacione su respuesta con la tasa de aceptación esperada del algoritmo.

R/T: El método de rejection sampling clásico sirve para poder generar datos que sigan una distribución que queremos, para eso se usa otra distribución más sencilla y se van generando candidatos, algunos se aceptan y otros se rechazan.

El problema es que si la constante (M) es muy grande, se van a rechazar muchos candidatos porque la distribución que estamos usando como referencia queda muy por encima de la distribución que realmente queremos. Por eso el método se vuelve menos eficiente y necesitamos hacer más intentos para conseguir un valor aceptado.

pregunta 2 Defina qué significa que una función de densidad p(x) sea log-cóncava. Dé dos ejemplos de distribuciones log-cóncavas y uno de una distribución que no lo sea (justifique brevemente por qué no lo es).

R/T: Una función de densidad \(p(x)\) es log-cóncava cuando su logaritmo, es decir,

\[\log p(x),\]

es una función cóncava. En otras palabras, al aplicar logaritmo a la densidad, la curva resultante debe tener forma cóncava.

Matemáticamente, si \(p(x)>0\), se cumple que

\[\frac{d^2}{dx^2}\log p(x)\leq 0.\]

Dos ejemplos de distribuciones log-cóncavas son la Normal y la Exponencial.

Para la Normal estándar,

\[p(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}.\]

Su logaritmo es

\[\log p(x)=-\frac{x^2}{2}-\frac{1}{2}\log(2\pi),\]

y como su segunda derivada es

\[\frac{d^2}{dx^2}\log p(x)=-1<0,\]

entonces es log-cóncava.

Para la distribución Exponencial,

\[p(x)=\lambda e^{-\lambda x},\qquad x\geq0.\]

Su logaritmo es

\[\log p(x)=\log(\lambda)-\lambda x,\]

cuya segunda derivada es \(0\), por lo que también es log-cóncava.

Un ejemplo de una distribución que no es log-cóncava es la distribución Cauchy. Su densidad es

\[p(x)=\frac{1}{\pi(1+x^2)}.\]

En este caso,

\[\log p(x)=-\log(\pi)-\log(1+x^2),\]

y su segunda derivada no es siempre negativa; de hecho,

\[\frac{d^2}{dx^2}\log p(x)=\frac{2(x^2-1)}{(1+x^2)^2}.\]

Cuando \(|x|>1\), esta expresión es positiva, por lo que \(\log p(x)\) deja de ser cóncava. Por eso la distribución Cauchy no es log-cóncava.

Pregunta 3 En ARS se construyen dos funciones a partir de los puntos de tangencia: la envolvente superior \(u(x)\) y la función de compresión (squeeze) \(l(x)\). Explique geométricamente cómo se construye cada una y demuestre por qué se cumple la desigualdad \(l(x) ≤ h(x) ≤ u(x)\) cuando \(h(x) = log p(x)\) es cóncava.

Sea \(h(x)=\log p(x)\) una función cóncava y sean \(x_1<x_2<\cdots<x_n\) los puntos de tangencia.

La envolvente superior \(u(x)\) se construye mediante las rectas tangentes:

\[t_i(x)=h(x_i)+h'(x_i)(x-x_i).\]

Como \(h(x)\) es cóncava, toda tangente queda por encima de la función:

\[h(x)\leq t_i(x),\]

por lo que

\[\boxed{h(x)\leq u(x)}.\]

La función squeeze \(l(x)\) se construye uniendo mediante segmentos de recta los puntos consecutivos \((x_i,h(x_i))\) y \((x_{i+1},h(x_{i+1}))\). Como una función cóncava queda por encima de las cuerdas que unen dos puntos de su gráfica,

\[l(x)\leq h(x).\]

Por tanto, se obtiene la desigualdad fundamental:

\[\boxed{l(x)\leq h(x)\leq u(x)}.\]

Geométricamente, \(u(x)\) actúa como un techo formado por tangentes y \(l(x)\) como un piso formado por segmentos que unen los puntos de tangencia.

pregunta 4 ¿Por qué el test de compresión (squeeze test) permite ahorrar cómputo en comparación con evaluarsiempre h(x)? Describa un escenario práctico (por ejemplo, dentro de un muestreador de Gibbs) donde este ahorro sea especialmente valioso.

El test de compresión permite ahorrar cómputo porque primero utiliza la función \(l(x)\), que es más sencilla de evaluar que \(h(x)\). Si el candidato cumple una condición que garantiza su aceptación usando \(l(x)\), no es necesario calcular \(h(x)\). Solo cuando el resultado no es concluyente se evalúa \(h(x)\).

Esto es especialmente útil en un muestreador de Gibbs cuando se generan muchos candidatos en cada iteración. Si \(h(x)\) requiere cálculos costosos, como sumas grandes, logaritmos o evaluaciones de una función de verosimilitud, el uso del squeeze evita realizar estos cálculos para muchos candidatos y reduce considerablemente el tiempo total de ejecución.

pregunta 5 El algoritmo se llama adaptativo. Explique qué se actualiza en cada iteración y por qué esto garantiza que la tasa de rechazo disminuya (o al menos no aumente) a medida que avanza el muestreo.

El algoritmo se llama adaptativo porque en cada iteración se actualizan los puntos de tangencia y, con ellos, las funciones \(u(x)\) y \(l(x)\). Cuando se obtiene un nuevo punto, se incorpora a la envolvente y esta se ajusta mejor a \(h(x)\).

Como \(u(x)\) se construye con las tangentes y cada nueva tangente toca a \(h(x)\), la envolvente se va acercando a la función objetivo. Por tanto, el área de la envolvente disminuye o se mantiene, haciendo que la probabilidad de aceptación aumente o se mantenga:

\[P(\text{aceptar})=\frac{\int e^{h(x)}\,dx}{\int e^{u(x)}\,dx}.\]

Así, al disminuir el área de \(u(x)\), la tasa de rechazo

\[P(\text{rechazar})=1-P(\text{aceptar})\]

disminuye o, como mínimo, no aumenta a medida que avanza el muestreo.

Derivaciones matemáticas

pregunta 6: Sea \(h(x)\) una función cóncava y sean \(x_{i} < x_{i+1}\) dos puntos de tangencia consecutivos. Las tangentes en esos puntos son:

\[t_{i}=h(x_{i})+h'(x_{i})(x-x_{i})\]

\[t_{i+1}=h(x_{i+1})+h'(x_{i+1})(x-x_{i+1})\]
a) Iguale ambas expresiones y despeje el punto de intersección \(z_{i+1}\). Compare su resultado con la fórmula usada en la función ars_build_envelope() del código visto en clase.

Solución:

Igualamos las tangentes, como en el punto de intersección ambas tienen el mismo valor, entonces:

\[t_i(z_{i+1})=t_{i+1}(z_{i+1})\]
Ahora igualamos y luego suatituimos:

\[h(x_{i})+h'(x_{i})(x-x_{i})=h(x_{i+1})+h'(x_{i+1})(x-x_{i+1})\]
Sustituimos.

\[h(x_{i})+h'(x_{i})(z_{i+1}-x_{i})=h(x_{i+1})+h'(x_{i+1})(z_{i+1}-x_{i+1})\]
Hacemos el producto.

\[ h(x_{i})+ h'(x_{i})z_{i+1}-h'(x_{i})x_{i}=h(x_{i+1})+h'(x_{i+1})z_{i+1}-h'(x_{i+1})x_{i+1}\]

\[h'(x_{i})z_{i+1}-h'(x_{i+1})z_{i+1}=h(x_{i+1})-h(x_{i})+h'(x_{i})x_{i}-h'(x_{i+1})x_{i+1}\]
Factorizamos

\[z_{i+1}[h'(x_{i})-h'(x_{i+1})]=h(x_{i+1})-h(x_{i})+x_{i}h'(x_{i})-x_{i+1}h'(x_{i+1})\]

Despejamos \(z_{i+1}\):

\[z_{i+1}= \frac{h(x_{i+1})-h(x_{i})+x_{i}h'(x_{i})-x_{i+1}h'(x_{i+1})}{h'(x_{i})-h'(x_{i+1})}\]

b) Explique por qué este punto \(z_{i+1}\) marca el límite entre dos tramos de la envolvente lineal a trozos \(u(x)\).

Para tramos asociandos a \(x_{i}\) tenemos:

\[u(x)=t_{i}(x)\]

y el tramos asociado a \(x_{i+1}\) tenemos:

\[u(x)=t_{i+1}(x)\]

y estas dos rectas se encuentran en \(z_{i+1}\) , entonces:

\[t_{i}(z_{i+1})=t_{i+1}(z_{i+1})\]

\(z_{i+1}\) es el limite entre dos tramos, así \(t_{i}(x)\)y \(t_{i+1}(x)\), toman el mismo valor, en otras palabras:

\[u(x)= \left\{ \begin{array}{cl} t_{i}(x) & : \ x \ \le \ z_{i+1 } \\ t_{i+1}(x) & : x \ \ge \ z_{i+1 } \end{array} \right.\]

pregunta 7: Para la distribución Normal estándar, \(h(x) =-\frac{x^{2}}{2}\) y \(h'(x) = -x\).

a) Calcule a mano la tangente a \(h(x)\) en el punto \(x_{0} = 1\).

Solución

La tangente en \(x_{i}\) es:

\[t_{i}(x)=h(x_{i})+h'(x_{i})(x-x_{i})\]

\[h(1) =-\frac{1^{2}}{2}= -\frac{1}{2}\]

\[h'(1) = -1\]

Entonces:

\[t_{i}(x)=-\frac{1}{2}+(-1)(x-1)\]

\[ = -\frac{1}{2}-x+1\]

\[t_{i}(x)= \frac{1}{2}-x\]

b) Si se usan los puntos iniciales \(x_{1} = -1.5\) y \(x_{2} = 1.0\), calcule el punto de intersección z entre ambas tangentes usando la fórmula de la pregunta 6.

\[z= \frac{h(x_{2})-h(x_{1})+x_{1}h'(x_{1})-x_{2}h'(x_{2})}{h'(x_{1})-h'(x_{2})}\]

Sabemos que:

\[h(x)=-\frac{x^{2}}{2} \ \ \ y \ \ \ h'(x) = -x \]

\[h(x_{1})=-\frac{(-1.5)^{2}}{2}=-\frac{2.25}{2}= -1.125\]

\[h(x_{2})=-\frac{(-1.0)^{2}}{2}=-0.5\]

\[h'(x_{1}) = -(-1.5)=1.5\]

\[h'(x_{2}) =-1.0\]

\[z= \frac{(-0.5)-(-1.125)+(-1.5)(1.5)-1.0(-1.0)}{(1.5)-(-1.0)}\]

\[z=\frac{-0.625}{2.5}\]

\[z=-0.25\]

c) ¿Por qué en este ejemplo se necesita al menos un punto con \(h'(x) > 0\) y otro con \(h'(x) < 0\) si el dominio es toda la recta real? (Pista: piense en el área bajo \(e^{u(x)}\) fuera de esos puntos).

Tenemos, \(h'(x)=-x\)

\[\therefore \left\{ \begin{array}{cl} si \ \ x < 0 & \Rightarrow h'(x)>0 \\ si \ \ x > 0 & \Rightarrow h'(x)<0 \end{array} \right.\]

Del punto anterior tenemos:

\[x_{1}=-1.5 \ \ y \ \ x_{2}=1.0\] 

\[\Rightarrow \left\{ \begin{array}{cl} h'(-1.5)= 1.5 > 0 \\ h'(1.0)= -1.0 < 0 \end{array} \right.\]

\(\therefore\) Se necesita al menos un punto con \(h'(x)>0\) y otro con \(h'(x)<0\), por que el dominio de la normal estandar es toda la recta real.

Pregunta 8 Demuestre que el área bajo \(e^{u_{i}(x)}\) entre dos puntos a y b, para un tramo con pendiente \(h'(x_{i})=m≠0\), es:

\[\int_{a}^{b}e^{h(x_{i})+m(x-x_{i})}dx=\frac{e^{h(x_{i})+m(b-x_{i})}-e^{h(x_{i})+m(a-x_{i})}}{m}\] Explique además qué ocurre con esta fórmula cuando \(m = 0\), y por qué el código maneja ese caso por separado (ver la función_segment_area en Python o ars_segment_area en R).

temos:

\[u_{i}(x)=e^{h(x_{i})+m(x-x_{i})}\]

con \(m=h'(x_{i}) ≠0\)

Partimos de:

\[I=\int_{a}^{b}e^{h(x_{i})+m(x-x_{i})}dx\]

Definimos:

\[t=h(x_{i})+m(x-x_{i})\]

Entonces,

\[\frac{dt}{dx}=m\]\

\[\therefore dt = mdx\]

\[ dx=\frac{dt}{m}\]

\[I=\int_{e^{h(x_{i})+m(a-x_{i})}}^{e^{h(x_{i})+m(b-x_{i})}}e^{t}\frac{dx}{m}\]

\[I=\frac{1}{m}\int_{e^{h(x_{i})+m(a-x_{i})}}^{e^{h(x_{i})+m(b-x_{i})}}e^{t}dx\]

\[I=\frac{1}{m}\left[ e^{t} \right]^{e^{h(x_{i})+m(b-x_{i})}}_{e^{h(x_{i})+m(a-x_{i})}}\]

\[I=\frac{1}{m} \ e^{h(x_{i})+m(b-x_{i})} -e^{h(x_{i})+m(a-x_{i})}\]

Así:

\[\int_{a}^{b}e^{h(x_{i})+m(x-x_{i})}dx=\frac{e^{h(x_{i})+m(b-x_{i})}-e^{h(x_{i})+m(a-x_{i})}}{m}\]

Entonces cuando \(m=0\)

\[u_{i}(x)=h(x_{i})+0(x-x_{i})\]

\[u_{i}(x)=h(x_{i})\]

\[\Rightarrow e^{u_{1}(x)}=e^{h(x_{i})}\]

\[\int_{a}^{b} e^{h(x_{i})}dx\]

\[e^{h(x_{i})}\int_{a}^{b}dx = b-a\]

\[\Rightarrow \int_{a}^{b} e^{h(x_{i})}dx\]

\[=b-a \ \ e^{h(x_{i})}\]

El codigo los maneja separados porque cuando \(m≠0\) tiene que, \(\frac{1}{m}\) y cuando \(m=0\) tiene que es, \(\frac{e^{h(x_{i})}-e^{h(x_{i})}}{m}=\frac{0}{0}\) lo cual es una indeterminación.

Aplicación y programación

Pregunta 9 Usando el código de ARS visto en clase (Python o R, a su elección), genere 5000 muestras de una distribución Gamma(shape = 5, scale = 1) usando: \[h(x) = (k - 1)·log(x) - x \ \ con \ \ k = 5\]

\[h'(x) = \frac{(k - 1)}{x - 1}\]\

# 1. Definir h(x) y su derivada

k <- 5

# h(x) = (k - 1) log(x) - x
h_gamma <- function(x) {
  (k - 1) * log(x) - x
}

# h'(x) = (k - 1)/x - 1
h_prime_gamma <- function(x) {
  (k - 1) / x - 1
}


# Crear ARS

AdaptiveRejectionSampler <- function(h, h_prime,
                                     domain = c(1e-6, Inf),
                                     init_points = c(1, 3, 6)) {

  x <- sort(init_points)

  hx <- sapply(x, h)
  hpx <- sapply(x, h_prime)


  # Construcción de la envolvente superior

  build_envelope <- function(x, hx, hpx) {

    n <- length(x)

    z <- numeric(n + 1)

    z[1] <- domain[1]
    z[n + 1] <- domain[2]

    for (i in 1:(n - 1)) {

      num <- hx[i + 1] - hx[i] -
        x[i + 1] * hpx[i + 1] +
        x[i] * hpx[i]

      den <- hpx[i] - hpx[i + 1]

      z[i + 1] <- num / den
    }

    return(z)
  }


  # Función de la envolvente superior

  u_function <- function(xv, i, x, hx, hpx) {

    hx[i] + hpx[i] * (xv - x[i])
  }


  # Área de cada segmento de la envolvente

  segment_area <- function(i, z, x, hx, hpx) {

    a <- z[i]
    b <- z[i + 1]

    slope <- hpx[i]
    h0 <- hx[i]
    x0 <- x[i]

    if (abs(slope) < 1e-12) {

      return(exp(h0) * (b - a))
    }

    upper_val <- if (is.finite(b)) {
      exp(h0 + slope * (b - x0))
    } else {
      0
    }

    lower_val <- if (is.finite(a)) {
      exp(h0 + slope * (a - x0))
    } else {
      0
    }

    return((upper_val - lower_val) / slope)
  }


  # Muestrear desde la envolvente

  sample_from_envelope <- function(z, x, hx, hpx, areas) {

    cum_areas <- cumsum(areas)

    total <- tail(cum_areas, 1)

    u1 <- runif(1, 0, total)

    i <- which(cum_areas >= u1)[1]

    previous_area <- if (i > 1) {
      cum_areas[i - 1]
    } else {
      0
    }

    target_area <- u1 - previous_area

    a <- z[i]

    slope <- hpx[i]

    h0 <- hx[i]

    x0 <- x[i]

    if (abs(slope) < 1e-12) {

      x_star <- a + target_area / exp(h0)

      return(list(x = x_star, segment = i))
    }

    lower_val <- if (is.finite(a)) {
      exp(h0 + slope * (a - x0))
    } else {
      0
    }

    val <- lower_val + slope * target_area

    x_star <- x0 + (log(val) - h0) / slope

    return(list(x = x_star, segment = i))
  }


  # Función para obtener las muestras

  sample <- function(n_samples) {

    samples <- numeric(0)

    while (length(samples) < n_samples) {

      # Construir la envolvente
      z <- build_envelope(x, hx, hpx)

      # Calcular las áreas
      areas <- sapply(
        1:length(x),
        segment_area,
        z = z,
        x = x,
        hx = hx,
        hpx = hpx
      )

      # Generar candidato
      candidate <- sample_from_envelope(
        z, x, hx, hpx, areas
      )

      x_star <- candidate$x

      i <- candidate$segment

      # Generar uniforme
      w <- runif(1)

      # Evaluar envolvente superior
      u_val <- u_function(
        x_star, i, x, hx, hpx
      )

      # Envolvente inferior

      if (x_star >= x[1] &&
          x_star <= x[length(x)]) {

        j <- findInterval(x_star, x)

        j <- max(1, min(j, length(x) - 1))

        weight <- (x_star - x[j]) /
          (x[j + 1] - x[j])

        l_val <- (1 - weight) * hx[j] +
          weight * hx[j + 1]

      } else {

        l_val <- -Inf
      }

      # Prueba de aceptación mediante squeeze

      if (is.finite(l_val) &&
          w <= exp(l_val - u_val)) {

        samples <- c(samples, x_star)

        next
      }

      # Evaluar h(x)

      h_val <- h(x_star)

      # Prueba de aceptación final

      if (w <= exp(h_val - u_val)) {

        samples <- c(samples, x_star)
      }

      # Adaptar la envolvente

      hp_val <- h_prime(x_star)

      position <- sum(x < x_star) + 1

      x <- append(x, x_star, after = position - 1)

      hx <- append(hx, h_val, after = position - 1)

      hpx <- append(hpx, hp_val, after = position - 1)
    }

    return(samples)
  }


  return(list(
    sample = sample
  ))
}


# Crear el muestreador ARS


sampler_gamma <- AdaptiveRejectionSampler(
  h_gamma,
  h_prime_gamma,
  domain = c(1e-6, Inf),
  init_points = c(1, 3, 6)
)

a) Reporte la media y varianza muestral obtenidas y compárelas con los valores teóricos \((E[X] = k, Var[X] = k)\).

set.seed(123)
# Generar 5000 muestras

muestras_gamma <- sampler_gamma$sample(5000)


media <- mean(muestras_gamma)

varianza <- var(muestras_gamma)

cat("Media muestral:", media, "\n")
## Media muestral: 4.994368
cat("Varianza muestral:", varianza, "\n")
## Varianza muestral: 4.983285

b) Grafique el histograma de las muestras junto con la densidad teórica de la Gamma(5,1). Adjunte la gráfica a su entrega.

# Histograma de las muestras ARS + densidad teorica
hist(
  muestras_gamma,
  probability = TRUE,
  breaks = 40,
  col = "lightgray",
  border = "white",
  main = "Muestras ARS vs. densidad teorica Gamma(5,1)",
  xlab = "x",
  ylab = "Densidad"
)

# densidad teorica Gamma(5,1)
curve(
  dgamma(x, shape = 5, scale = 1),
  from = 0,
  to = max(muestras_gamma),
  add = TRUE,
  col = "red",
  lwd = 2
)

# leyenda
legend(
  "topright",
  legend = "Densidad teorica Gamma(5,1)",
  col = "red",
  lwd = 2,
  bty = "n"
)

R/T: En el gráfico se observa que el histograma de las 5000 muestras generadas mediante el algoritmo ARS presenta un ajuste aproximado a la densidad teórica de la distribución Gamma(5,1). La mayor concentración de observaciones se encuentra alrededor de \(x=4\), donde la distribución alcanza su mayor densidad. A partir de este punto, la frecuencia de las observaciones disminuye y se presenta una cola hacia la derecha, característica de la distribución Gamma. Las pequeñas diferencias entre las barras del histograma y la curva teórica se deben a la variabilidad aleatoria del muestreo. En general, el gráfico muestra que el ARS reproduce adecuadamente la forma de la distribución Gamma(5,1).

c) Pruebe con distintos puntos iniciales (por ejemplo, todos muy cercanos entre sí vs. bien distribuidos). ¿Cambia el número de evaluaciones de h(x) necesarias para generar las 5000 muestras? Explique por qué.

contador_h <- 0

h_gamma_contador <- function(x) {

  contador_h <<- contador_h + 1

  4 * log(x) - x
}
contador_h <- 0

sampler_cercanos <- AdaptiveRejectionSampler(
  h_gamma_contador,
  h_prime_gamma,
  domain = c(1e-6, Inf),
  init_points = c(3.8, 4.0, 4.2)
)
muestras_cercanos <- sampler_cercanos$sample(5000)
evaluaciones_cercanos <- contador_h
cat(
  "Puntos cercanos:",
  evaluaciones_cercanos,
  "evaluaciones de h(x)\n"
)
## Puntos cercanos: 43 evaluaciones de h(x)
contador_h <- 0

sampler_distribuidos <- AdaptiveRejectionSampler(
  h_gamma_contador,
  h_prime_gamma,
  domain = c(1e-6, Inf),
  init_points = c(1, 3, 6)
)
contador_h <- 0
muestras_distribuidas <- sampler_distribuidos$sample(5000)
evaluaciones_distribuidas <- contador_h
cat(
  "Puntos bien distribuidos:",
  evaluaciones_distribuidas,
  "evaluaciones de h(x)\n"
)
## Puntos bien distribuidos: 40 evaluaciones de h(x)

R/T: Al comparar diferentes puntos iniciales se observa que el número de evaluaciones de \(h(x)\) puede cambiar. Esto ocurre porque los puntos iniciales determinan la envolvente inicial utilizada por el algoritmo ARS. Cuando los puntos están bien distribuidos, la envolvente puede representar mejor la función desde el comienzo, por lo que se necesitan menos actualizaciones y evaluaciones adicionales de \(h(x)\). En cambio, cuando los puntos están muy concentrados, la envolvente inicial puede ser menos precisa en algunas regiones y el algoritmo necesita adaptarse agregando nuevos puntos. Por esta razón, la elección de los puntos iniciales puede afectar la eficiencia del método.

Pregunta 10 Considere la densidad Beta(2,2), cuya log-densidad (sin normalizar) es \(h(x) = log(x) + log(1-x)\) en el intervalo \((0,1)\).

a) Verifique que \(h(x)\) es cóncava en (0,1) calculando \(h''(x)\).

Para verificar que \(h(x)\) es cóncava, calculamos su segunda derivada. Partimos de

\[h(x)=\log(x)+\log(1-x), \qquad 0<x<1.\]

La primera derivada es

\[h'(x)=\frac{1}{x}-\frac{1}{1-x}.\]

Derivando nuevamente,

\[h''(x)=-\frac{1}{x^2}-\frac{1}{(1-x)^2}.\]

Como \(0<x<1\), tanto \(x^2\) como \((1-x)^2\) son positivos. Por lo tanto,

\[h''(x)<0\]

para todo \(x\in(0,1)\). En consecuencia, \(h(x)\) es estrictamente cóncava en el intervalo \((0,1)\).

b) Adapte el código de ARS para este caso (dominio finito, sin necesitar las condiciones de \(\frac{h'(x1)>0} {h'(xk)<0}\) que sí se requieren en dominios infinitos). Genere 3000 muestras y compare la media muestral con la media teórica de la \(Beta(2,2)\), que es 0.5.

# ARS para Beta(2,2)
# Dominio finito: (0,1)

AdaptiveRejectionSampler_Finito <- function(h, h_prime,
                                            domain = c(0, 1),
                                            init_points = c(0.25, 0.5, 0.75)) {

  # Puntos iniciales
  x <- sort(init_points)

  # Evaluaciones iniciales
  hx <- sapply(x, h)
  hpx <- sapply(x, h_prime)


  # Construccion de la envolvente superior

  build_envelope <- function(x, hx, hpx) {

    n <- length(x)

    z <- numeric(n + 1)

    # Como el dominio es finito,
    # los extremos son 0 y 1
    z[1] <- domain[1]
    z[n + 1] <- domain[2]

    if (n > 1) {

      for (i in 1:(n - 1)) {

        num <- hx[i + 1] - hx[i] -
          x[i + 1] * hpx[i + 1] +
          x[i] * hpx[i]

        den <- hpx[i] - hpx[i + 1]

        z[i + 1] <- num / den
      }
    }

    return(z)
  }

  # Funcion de la envolvente superior u(x)

  u_function <- function(xv, i, x, hx, hpx) {

    hx[i] + hpx[i] * (xv - x[i])
  }


  # Area bajo exp(u(x))

  segment_area <- function(i, z, x, hx, hpx) {

    a <- z[i]
    b <- z[i + 1]

    slope <- hpx[i]
    h0 <- hx[i]
    x0 <- x[i]

    # Si la pendiente es aproximadamente cero
    if (abs(slope) < 1e-12) {

      return(
        exp(h0) * (b - a)
      )
    }

    upper_val <- exp(
      h0 + slope * (b - x0)
    )

    lower_val <- exp(
      h0 + slope * (a - x0)
    )

    return(
      (upper_val - lower_val) / slope
    )
  }


  # Generacion de candidato desde la envolvente

  sample_from_envelope <- function(z, x, hx, hpx, areas) {

    cum_areas <- cumsum(areas)

    total <- tail(cum_areas, 1)

    u1 <- runif(1, 0, total)

    i <- which(cum_areas >= u1)[1]

    previous_area <- if (i > 1) {
      cum_areas[i - 1]
    } else {
      0
    }

    target_area <- u1 - previous_area

    a <- z[i]

    slope <- hpx[i]
    h0 <- hx[i]
    x0 <- x[i]

    # Caso de pendiente aproximadamente cero
    if (abs(slope) < 1e-12) {

      x_star <- a +
        target_area / exp(h0)

      return(
        list(
          x = x_star,
          segment = i
        )
      )
    }

    lower_val <- exp(
      h0 + slope * (a - x0)
    )

    val <- lower_val +
      slope * target_area

    x_star <- x0 +
      (log(val) - h0) / slope

    return(
      list(
        x = x_star,
        segment = i
      )
    )
  }


  # Generacion de muestras

  sample <- function(n_samples) {

    samples <- numeric(0)

    while (length(samples) < n_samples) {

      # Construir envolvente
      z <- build_envelope(
        x,
        hx,
        hpx
      )

      # Areas de cada segmento
      areas <- sapply(
        1:length(x),
        segment_area,
        z = z,
        x = x,
        hx = hx,
        hpx = hpx
      )

      # Generar candidato
      candidate <- sample_from_envelope(
        z,
        x,
        hx,
        hpx,
        areas
      )

      x_star <- candidate$x

      i <- candidate$segment

      # Variable uniforme para aceptar/rechazar
      w <- runif(1)

      # Valor de la envolvente
      u_val <- u_function(
        x_star,
        i,
        x,
        hx,
        hpx
      )

      # Evaluacion de h(x)

      h_val <- h(x_star)

      # Probabilidad de aceptacion
      if (w <= exp(h_val - u_val)) {

        samples <- c(
          samples,
          x_star
        )

      }

      # Agregar el nuevo punto para adaptar la envolvente

      hp_val <- h_prime(x_star)

      position <- sum(x < x_star) + 1

      x <- append(
        x,
        x_star,
        after = position - 1
      )

      hx <- append(
        hx,
        h_val,
        after = position - 1
      )

      hpx <- append(
        hpx,
        hp_val,
        after = position - 1
      )
    }

    return(samples)
  }


  return(
    list(
      sample = sample
    )
  )
}

Definimos la función beta y su derivada.

h_beta <- function(x) {
  log(x) + log(1 - x)
}
h_prime_beta <- function(x) {
  1 / x - 1 / (1 - x)
}

Creamos ARS

sampler_beta <- AdaptiveRejectionSampler_Finito(
  h_beta,
  h_prime_beta,
  domain = c(0, 1),
  init_points = c(0.25, 0.5, 0.75)
)
muestras_beta <- sampler_beta$sample(3000)

Media muestral.

media_muestral <- mean(muestras_beta)

media_muestral
## [1] 0.4935299

Media teorica.

media_teorica <- 2 / (2 + 2)

media_teorica
## [1] 0.5

Comparamos las medias.

cat(
  "Media muestral:",
  media_muestral,
  "\n"
)
## Media muestral: 0.4935299
cat(
  "Media teorica:",
  media_teorica,
  "\n"
)
## Media teorica: 0.5

Calculamos el error, para mirar que tan lejos esta la media mustral de la media teorica.

error <- abs(
  media_muestral - media_teorica
)

cat(
  "Error absoluto:",
  error,
  "\n"
)
## Error absoluto: 0.00647009
hist(
  muestras_beta,
  probability = TRUE,
  breaks = 30,
  col = "lightgray",
  border = "white",
  main = "Muestras ARS de una Beta(2,2)",
  xlab = "x",
  ylab = "Densidad"
)
curve(
  dbeta(
    x,
    shape1 = 2,
    shape2 = 2
  ),
  from = 0,
  to = 1,
  add = TRUE,
  col = "red",
  lwd = 2
)

R/T:En el gráfico se observa que el histograma de las 3000 muestras generadas mediante ARS se ajusta de manera aproximada a la densidad teórica de la distribución Beta(2,2). Las observaciones presentan una mayor concentración alrededor de \(x=0.5\) y disminuyen hacia los extremos del intervalo \((0,1)\). Además, se evidencia una forma aproximadamente simétrica alrededor de \(0.5\), tal como se espera para esta distribución. Las pequeñas diferencias entre las barras del histograma y la curva teórica se deben a la variabilidad aleatoria del muestreo. En general, el gráfico muestra que el algoritmo ARS genera adecuadamente muestras de una Beta(2,2).

Pregunta 11 Reto opcional

El método Adaptive Rejection Metropolis Sampling (ARMS) es una versión ampliada del método ARS, propuesta por Gilks, Best y Tan en 1995.

El método ARS funciona bien cuando la distribución tiene una forma log-cóncava. Sin embargo, no todas las distribuciones cumplen esta condición, especialmente algunas que presentan colas pesadas o formas más complejas.

ARS utiliza las tangentes de la función

\[h(x) = \log(p(x))\]

para construir una distribución auxiliar que sirve como propuesta. El problema es que, cuando \(h(x)\) no es cóncava en todo el dominio, esta propuesta puede no ser suficiente para obtener muestras correctas.

Para solucionar esto, ARMS mantiene la idea adaptativa de ARS, pero añade un paso de Hastings-Metropolis. En términos sencillos, el procedimiento es el siguiente:

  1. Se construye una propuesta utilizando ARS.
  2. Se genera un posible nuevo valor.
  3. Se compara este valor con la distribución objetivo.
  4. El nuevo valor se acepta o se rechaza mediante el criterio de Hastings-Metropolis.
  5. Si se rechaza, la cadena conserva el valor anterior.
  6. La propuesta continúa ajustándose con la información obtenida.

Por lo tanto, la idea principal puede resumirse como:

\[\text{ARMS} = \text{ARS} + \text{Hastings-Metropolis}\]

La ventaja de ARMS es que permite trabajar con distribuciones que no son completamente log-cóncavas. Por esta razón, resulta útil para generar muestras de distribuciones más complicadas, como algunas distribuciones condicionales utilizadas en el muestreador de Gibbs.

En resumen:

Pregunta 12 Investigue y Compare

Los métodos para generar muestras tienen diferentes ventajas y limitaciones. La elección depende del tipo de distribución que se quiera simular y de si es fácil calcular su densidad o su función de distribución acumulada.

Método Ventajas Desventajas
Transformada inversa Es sencilla de comprender y programar. Si se conoce fácilmente la inversa de la función de distribución acumulada, permite generar muestras de forma directa y sin rechazos. Para muchas distribuciones no es sencillo encontrar la inversa de la función de distribución acumulada. En algunos casos, incluso puede no existir una forma fácil de calcularla.
Rejection sampling clásico Es un método general y relativamente fácil de implementar. No requiere conocer la inversa de la función de distribución acumulada. Su eficiencia depende de qué tan bien la distribución propuesta se parezca a la distribución objetivo. Si el ajuste es malo, se producen muchos rechazos y el método puede volverse lento.
ARS Puede ser muy eficiente cuando la densidad es log-cóncava. Además, la propuesta se va ajustando durante el proceso, lo que puede mejorar su rendimiento. Solo puede utilizarse directamente cuando la densidad cumple ciertas condiciones. También requiere conocer la derivada del logaritmo de la densidad, por lo que su implementación es más complicada.
Funciones nativas (rnorm(), rgamma(), numpy.random, etc.) Son rápidas, eficientes y fáciles de utilizar. Cuando la distribución ya está disponible en R o Python, normalmente representan la opción más práctica. Solo pueden utilizarse directamente para las distribuciones que ya están implementadas. Además, no permiten observar con detalle cómo funciona el proceso de generación de las muestras.

¿Cuándo sería útil implementar ARS manualmente?

En la mayoría de los casos prácticos, es más conveniente utilizar funciones ya programadas y optimizadas, como rnorm() o rgamma(). Sin embargo, implementar ARS manualmente puede ser útil cuando se trabaja con una distribución personalizada que cumple la condición de ser log-cóncava, pero que no está disponible directamente en R o Python.

También puede ser una buena opción con fines académicos, ya que permite comprender mejor cómo funcionan los métodos de rechazo, las envolventes, la función squeeze y la adaptación del algoritmo.

Además, ARS puede ser útil en problemas de inferencia bayesiana, especialmente dentro de un muestreador de Gibbs cuando alguna de las distribuciones condicionales es log-cóncava.

Conclusión

No existe un método que sea mejor para todos los casos. La elección depende de las características de la distribución y del objetivo del problema.

En general, ARS resulta especialmente interesante cuando se trabaja con densidades personalizadas log-cóncavas o cuando se desea comprender mejor el funcionamiento de los métodos de simulación.