Untitled

Metodos de simulación para la distribución posterior

Información

Prior Beta(1,1)

Verosimiltud Binomial (20,14)

Posterior verdadera Beta(15,7)

library(ggplot2)

# Parámetros de la distribución beta
a <- 1 + 14
b <- 1 + 6

# Generar datos desde una distribución beta
x <- seq(0, 1, length = 1000)
y <- dbeta(x, shape1 = a, shape2 = b)

# Crear un marco de datos con los datos
df <- data.frame(x, y)

# Graficar la distribución beta con ggplot2
ggplot(df, aes(x, y)) +
  geom_line(color = "blue", size = 1) +
  labs(title = "Distribución posterior verdadera Beta(15,7)", x = "Valores", y = "Densidad") 

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

media = 0.68182

sigma = 0.00943

Rejection Sampling

set.seed(123)
n <- 20
h <- 14 # numero de exitos
it <- 1000  # numero de interaciones
samples <- numeric()

for (i in 1:it) {
  prior_sample <- rbeta(1, 1, 1) # genera un numero aleatorio de una distribucion beta(1,1)
  heads_sample <- rbinom(1, n, prior_sample) # genera el numero de exitos de una distribucion bernoulli 
  if (heads_sample == h) { # verifica condición si el numero de exitos de la simulación coincide con el número de exitos observados
    samples <- c(samples, prior_sample) # guarda el valor del parámetro si cumple la condición
  }
}

# Crear un marco de datos con los datos
df <- data.frame(samples)

# Graficar el histograma y la densidad
ggplot(df, aes(x = samples)) +
  geom_density(color = "red", size = 1) +  # Agregar densidad
  labs(title = "Rejection Sampling - Distribución de la Muestra 1000 de iteraciones",
       x = "Valor de la Muestra",
       y = "Densidad")+
  xlim(0, 1)

mean(samples)

## [1] 0.6302586

var(samples)

## [1] 0.01042802

set.seed(123)
n <- 20
h <- 14 # numero de exitos
it <- 10000  # numero de interaciones
samples <- numeric()

for (i in 1:it) {
  prior_sample <- rbeta(1, 1, 1) # genera un numero aleatorio de una distribucion beta(1,1)
  heads_sample <- rbinom(1, n, prior_sample) # genera el numero de exitos de una distribucion bernoulli 
  if (heads_sample == h) { # verifica condición si el numero de exitos de la simulación coincide con el número de exitos observados
    samples <- c(samples, prior_sample) # guarda el valor del parámetro si cumple la condición
  }
}

# Crear un marco de datos con los datos
df <- data.frame(samples)

# Graficar el histograma y la densidad
ggplot(df, aes(x = samples)) +
  geom_density(color = "red", size = 1) +  # Agregar densidad
  labs(title = "Rejection Sampling - Distribución de la Muestra 10000 de iteraciones",
       x = "Valor de la Muestra",
       y = "Densidad")+
  xlim(0, 1)

mean(samples)

## [1] 0.6809181

var(samples)

## [1] 0.009268504

set.seed(123)
n <- 20
h <- 14 # numero de exitos
it <- 100000  # numero de interaciones
samples <- numeric()

for (i in 1:it) {
  prior_sample <- rbeta(1, 1, 1) # genera un numero aleatorio de una distribucion beta(1,1)
  heads_sample <- rbinom(1, n, prior_sample) # genera el numero de exitos de una distribucion bernoulli 
  if (heads_sample == h) { # verifica condición si el numero de exitos de la simulación coincide con el número de exitos observados
    samples <- c(samples, prior_sample) # guarda el valor del parámetro si cumple la condición
  }
}

# Crear un marco de datos con los datos
df <- data.frame(samples)

# Graficar el histograma y la densidad
ggplot(df, aes(x = samples)) +
  geom_density(color = "red", size = 1) +  # Agregar densidad
  labs(title = "Rejection Sampling - Distribución de la Muestra 100000 de iteraciones",
       x = "Valor de la Muestra",
       y = "Densidad")+
  xlim(0, 1)

mean(samples)

## [1] 0.6827453

var(samples)

## [1] 0.009610484

set.seed(123)
n <- 20
h <- 14
it <- 100000
samples <- numeric()
mean_evolution <- numeric()  # Vector para almacenar la evolución de la media

for (i in 1:it) {
  prior_sample <- rbeta(1, 1, 1)
  heads_sample <- rbinom(1, n, prior_sample)
  if (heads_sample == h) {
    samples <- c(samples, prior_sample)
  }
  
  # Calcular la media en cada iteración y almacenarla
  mean_evolution <- c(mean_evolution, mean(samples))
}

# Crear un marco de datos con los datos de la evolución de la media
df_evolution <- data.frame(iteration = 1:it, mean = mean_evolution)

# Graficar la evolución de la media
ggplot(df_evolution, aes(x = iteration, y = mean)) +
  geom_line(color = "red", size = 1) +
  labs(title = "Sampling Rejection - Evolución de la Media a lo largo de las Iteraciones",
       x = "Iteración",
       y = "Media")

## Warning: Removed 14 rows containing missing values (`geom_line()`).

Metropolis Hastings

# Establece la semilla para la reproducibilidad
set.seed(123)

# Parámetros
n <- 20
h <- 14
it <- 1000
samples <- numeric()
current_sample <- 0.5
rejection_count <- 0

# Bucle de simulación
for (i in 1:it) {
  prop_sample <- pmax(0, pmin(1, rnorm(1, current_sample, 0.1)))
  current_prob <- dbinom(h, n, current_sample)
  prop_prob <- dbinom(h, n, prop_sample)
  accept_ratio <- prop_prob / current_prob
  rand <- runif(1)
  
  if (rand <= accept_ratio) {
    current_sample <- prop_sample
  } else {
    rejection_count <- rejection_count + 1
  }
  
  samples <- c(samples, current_sample)
}

# Crear un marco de datos con los datos
df <- data.frame(samples)

# Graficar la densidad kernel
ggplot(df, aes(x = samples)) +
  geom_density(color = "blue", size = 1) +
  labs(title = "Metropolis Hastings - Distribución de la Muestra 1000 iteraciones",
       x = "Valor de la Muestra",
       y = "Densidad") +
  xlim(0, 1)

mean(samples)

## [1] 0.6844295

var(samples)

## [1] 0.0101563

# Establece la semilla para la reproducibilidad
set.seed(123)

# Parámetros
n <- 20
h <- 14
it <- 10000
samples <- numeric()
current_sample <- 0.5
rejection_count <- 0

# Bucle de simulación
for (i in 1:it) {
  prop_sample <- pmax(0, pmin(1, rnorm(1, current_sample, 0.1)))
  current_prob <- dbinom(h, n, current_sample)
  prop_prob <- dbinom(h, n, prop_sample)
  accept_ratio <- prop_prob / current_prob
  rand <- runif(1)
  
  if (rand <= accept_ratio) {
    current_sample <- prop_sample
  } else {
    rejection_count <- rejection_count + 1
  }
  
  samples <- c(samples, current_sample)
}

# Crear un marco de datos con los datos
df <- data.frame(samples)

# Graficar la densidad kernel
ggplot(df, aes(x = samples)) +
  geom_density(color = "blue", size = 1) +
  labs(title = "Metropolis Hastings - Distribución de la Muestra 10000 iteraciones",
       x = "Valor de la Muestra",
       y = "Densidad") +
  xlim(0, 1)

mean(samples)

## [1] 0.6804315

var(samples)

## [1] 0.00990304

# Establece la semilla para la reproducibilidad
set.seed(123)

# Parámetros
n <- 20
h <- 14
it <- 100000
samples <- numeric()
current_sample <- 0.5
rejection_count <- 0

# Bucle de simulación
for (i in 1:it) {
  prop_sample <- pmax(0, pmin(1, rnorm(1, current_sample, 0.1)))
  current_prob <- dbinom(h, n, current_sample)
  prop_prob <- dbinom(h, n, prop_sample)
  accept_ratio <- prop_prob / current_prob
  rand <- runif(1)
  
  if (rand <= accept_ratio) {
    current_sample <- prop_sample
  } else {
    rejection_count <- rejection_count + 1
  }
  
  samples <- c(samples, current_sample)
}

# Crear un marco de datos con los datos
df <- data.frame(samples)

# Graficar la densidad kernel
ggplot(df, aes(x = samples)) +
  geom_density(color = "blue", size = 1) +
  labs(title = "Metropolis Hastings - Distribución de la Muestra 100000 iteraciones",
       x = "Valor de la Muestra",
       y = "Densidad") +
  xlim(0, 1)

mean(samples)

## [1] 0.6809516

var(samples)

## [1] 0.009296026

# Establece la semilla para la reproducibilidad
set.seed(123)

# Parámetros
n <- 20
h <- 14
it <- 100000
samples <- numeric()
current_sample <- 0.5
rejection_count <- 0
mean_evolution <- numeric()  # Vector para almacenar la evolución de la media

# Bucle de simulación
for (i in 1:it) {
  prop_sample <- pmax(0, pmin(1, rnorm(1, current_sample, 0.1)))
  current_prob <- dbinom(h, n, current_sample)
  prop_prob <- dbinom(h, n, prop_sample)
  accept_ratio <- prop_prob / current_prob
  rand <- runif(1)
  
  if (rand <= accept_ratio) {
    current_sample <- prop_sample
  } else {
    rejection_count <- rejection_count + 1
  }
  
  samples <- c(samples, current_sample)
  
  # Calcular la media en cada iteración y almacenarla
  mean_evolution <- c(mean_evolution, mean(samples))
}

# Crear un marco de datos con los datos de la evolución de la media
df_evolution <- data.frame(iteration = 1:it, mean = mean_evolution)

# Graficar la evolución de la media
ggplot(df_evolution, aes(x = iteration, y = mean)) +
  geom_line(color = "blue", size = 1) +
  labs(title = "Metropolis Hastings - Evolución de la Media a lo largo de las Iteraciones",
       x = "Iteración",
       y = "Media")