Ejemplo Aceptacion rechazo

Metodo aceptación rechazo

Ejemplo

Generar una muestra de \(X \sim Beta(2, 4)\)

\[ \begin{aligned*} f(x) &= \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1}\\ f(x) &= \frac{\Gamma(2+4)}{\Gamma(2)\Gamma(4)}x^1(1-x)^3\\ f(x) &= 20x(1-x)^3 \end{aligned*} \]

Utilizando la distribución instrumental $ YU(0,1)$

\[ \begin{aligned*} \text{Note que:}\\ \frac{f(x)}{g(x)} &=20x(1-x)^3 \\ \text{Derivando}:\\ \frac{d}{dx}\left(\frac{f(x)}{g(x)}\right)&=20[(1-x)^3-3(1-x)^2] \\ \text{igualando a cero}: \\ (1-x)^3-3(1-x)^2&=0 \\ \text{Luego:}: \\ (1-x)^3 &=3(1-x)^2 \\ \text{simplificando:}\\ (1-x)&=3x \\ 1&=4x\\ x& = 1/4 \end{aligned*} \] Encontrado el x evaluar en la funcion:

\[ \begin{aligned*} \text{Note que:}\\ \frac{f(x)}{g(x)}&=20x(1-x)^3 \\ \frac{f(x)}{g(x)}&=20(0.25)(1-(0.25))^3 \\ \frac{f(x)}{g(x)}&= 2.1094 = c\\ \\ c &= 2.1094 \\ \text{Finalmente:}\\ \frac{1}{c}\frac{f(x)}{g(x)} &= \frac{1}{2.1094}20x(1-x)^3 \\ \frac{1}{c}\frac{f(x)}{g(x)} &= 9.4814 x(1-x)^3\\ \end{aligned*} \]

Si se quiere un tamaño de muestra de n = 10000, entonces se ocupan alrededor de 21094.0 simulaciones esto por que

\[ nc = 2.1094(10000) =21094.0 \]

feg <- function (x) 20*x*(1-x)^3 # objetio/ instruental
valC<-optimize(feg, interval=c(0,1), maximum=T)$objective
fegC <- function (x) (20*x*(1-x)^3)*(1/valC) # Beta (2,4)

n <- 10000 # tamaño de la muestra
k <- 0 # contador de la aceptación
i <- 0 # iterations
x <- numeric(n)

while(k <= n)
  {
  y <- runif(1) # Paso 1: se genera un valor de la instrumental
  u <- runif(1) # Paso 2: se genera un valor de una uniforme
  i <- i + 1 # iteraciones
  
  if(u< fegC(y) ) # paso 3
    {
    # se acepta x
    k <- k + 1 #se actualiza el contador
    x[k] <- y
  }
}

# contador
k

## [1] 10001

ks.test(x,"pbeta",2,4)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.0075956, p-value = 0.611
## alternative hypothesis: two-sided

hist(x,breaks=40,freq=F,main="Muestra generada",col="blue")
sop.x<-seq(min(x),max(x),0.01) # soporte de f
lines(sop.x, dbeta(sop.x,2,4),col="red",lwd=4)
legend("topright",
       legend = c("Histograma", "Beta(2,4)"),
       col = c("blue", "red"),
       lwd = c(10, 3),
       bty = "n")

Ejemplo de clase 25/03/26

repeat{
  y <- runif(1, min=-1, max=1)
  u <- runif(1)
  
  if (u <= 1 - y^2){
    x<-y
    break
  }
}

x

## [1] 0.5210733

n <- 10000
x <- numeric(n)

for (i in 1:n){
  repeat{
    y <- runif(1, min=-1, max=1)
    u <- runif(1)
    
    if (u<=1 -y^2){
      x[i] <- y
      break
    }
  }
}


hist(x, col="blue", probability = TRUE, breaks = 50,
     border = "white", main="Aceptación - Rechazo")

curve((3/4)*(1-x^2),from =-1, to =1,
      add=TRUE, lwd =5,col="red")

legend("topright",
       legend = c("Histograma", expression(f(x) == frac(3,4)*(1 - x^2))),
       col = c("blue", "red"),
       lwd = c(10, 3),
       bty = "n")

Ejercicio 2.7 Casella

El Ejemplo 2.2 no proporcionó un algoritmo general para simular variables aleatorias beta \(\mathrm{Be}(\alpha, \beta)\). Sin embargo, podemos construir un algoritmo sencillo (de prueba) basado en el método de Aceptación–Rechazo, utilizando como distribución instrumental la distribución uniforme \(U[0,1]\), cuando tanto \(\alpha\) como \(\beta\) son mayores que 1. (La función genérica no impone esta restricción).

El límite superior \(M\) es entonces el máximo de la densidad beta, obtenido por ejemplo mediante (o su alias ).

# Encontrar el valor Máximo
optimize (f=function(x){dbeta(x,2.7,6.3)},interval=c(0,1),maximum=TRUE)$objective

## [1] 2.669744

alfa = 2.7
beta = 6.3
Moda = (alfa-1)/(alfa+beta-2)
Moda

## [1] 0.2428571