Monte Carlo

Integración por Monte Carlo (MC) — Teoría y paso a paso

1) Idea central

Queremos aproximar una integral (a veces difícil de resolver analíticamente) usando aleatoriedad.

Monte Carlo convierte una integral en un valor esperado y luego lo estima con un promedio muestral.

2) Caso básico: integral en 1D sobre un intervalo

Sea:

\[ I=\int_a^b f(x)\,dx \]

Si tomamos una variable aleatoria:

\[ X \sim \mathrm{Uniforme}(a,b) \]

entonces:

\[ \mathbb{E}[f(X)] = \frac{1}{b-a}\int_a^b f(x)\,dx \]

Por lo tanto:

\[ I = (b-a)\,\mathbb{E}[f(X)] \]

Estimador Monte Carlo

Genera \(n\) muestras i.i.d. \(X_1,\dots,X_n \sim U(a,b)\). Define:

\[ \widehat{I}_n = (b-a)\,\frac{1}{n}\sum_{i=1}^n f(X_i) \]

3) Paso a paso (algoritmo)

Define la integral \(I=\int_a^b f(x)\,dx\).
Simula \(X_i \sim U(a,b)\) para \(i=1,\dots,n\).
Evalúa \(f(X_i)\).
Promedia:

\[ \overline{f}=\frac{1}{n}\sum_{i=1}^n f(X_i) \]
Multiplica por el largo del intervalo:

\[ \widehat{I}_n=(b-a)\overline{f} \]

4) Propiedades importantes

4.1) Insesgadez

Como \(\mathbb{E}[\overline{f}] = \mathbb{E}[f(X)]\), se cumple:

\[ \mathbb{E}[\widehat{I}_n] = I \]

4.2) Varianza del estimador

Sea \(\sigma_f^2 = \mathrm{Var}(f(X))\). Entonces:

\[ \mathrm{Var}(\widehat{I}_n) = (b-a)^2\,\frac{\sigma_f^2}{n} \]

4.3) Error típico (desviación estándar)

El error Monte Carlo decrece como \(1/\sqrt{n}\):

\[ \mathrm{SD}(\widehat{I}_n) = (b-a)\,\frac{\sigma_f}{\sqrt{n}} \]

Regla práctica: para reducir el error a la mitad, necesitas ~4 veces más simulaciones.

5) Intervalo de confianza (aprox. normal)

Por el TCL (para \(n\) grande):

\[ \widehat{I}_n \approx \mathcal{N}\left(I,\ (b-a)^2\frac{\sigma_f^2}{n}\right) \]

Como \(\sigma_f\) es desconocida, se estima con la desviación estándar muestral:

\[ s_f^2 = \frac{1}{n-1}\sum_{i=1}^n\left(f(X_i)-\overline{f}\right)^2 \]

Entonces:

\[ \widehat{\mathrm{SE}}(\widehat{I}_n)=(b-a)\,\frac{s_f}{\sqrt{n}} \]

Un IC del 95%:

\[ \widehat{I}_n \pm 1.96\,\widehat{\mathrm{SE}}(\widehat{I}_n) \]

6) Generalización: integral como esperanza con una densidad cualquiera (Importance Sampling)

Para una integral en 1D o multi-D:

\[ I=\int g(x)\,dx \]

Escoge una densidad \(p(x)\) tal que \(p(x)>0\) donde \(g(x)\neq 0\). Reescribe:

\[ I=\int \frac{g(x)}{p(x)}\,p(x)\,dx = \mathbb{E}_p\left[\frac{g(X)}{p(X)}\right] \]

Si \(X_i \sim p\), el estimador es:

\[ \widehat{I}_n = \frac{1}{n}\sum_{i=1}^n \frac{g(X_i)}{p(X_i)} \]

Idea: elegir \(p\) “parecida” a \(g\) reduce la varianza.

7) Integral en varias dimensiones (multi-D)

Sea:

\[ I=\int_{[0,1]^d} f(\mathbf{x})\,d\mathbf{x} \]

Si \(\mathbf{U}\sim U([0,1]^d)\), entonces:

\[ I=\mathbb{E}[f(\mathbf{U})] \]

Estimador:

\[ \widehat{I}_n=\frac{1}{n}\sum_{i=1}^n f(\mathbf{U}_i) \]

Si el dominio es un hipercubo \([a_1,b_1]\times\cdots\times[a_d,b_d]\), el volumen es:

\[ V=\prod_{j=1}^d (b_j-a_j) \]

\[ \widehat{I}_n = V\cdot \frac{1}{n}\sum_{i=1}^n f(\mathbf{X}_i),\quad \mathbf{X}_i\sim U(\text{dominio}) \]

8) Método geométrico (hit-or-miss) — opcional

Si conoces un cota \(0\le f(x)\le M\) en \([a,b]\), define el rectángulo de área \((b-a)M\).

Simula puntos \((X,Y)\) uniformes en \([a,b]\times[0,M]\). Define indicador:

\[ Z=\mathbf{1}\{Y\le f(X)\} \]

Entonces:

\[ I=\int_a^b f(x)\,dx = (b-a)M\,\mathbb{E}[Z] \]

Estimador:

\[ \widehat{I}_n=(b-a)M\,\frac{1}{n}\sum_{i=1}^n Z_i \]

9) Ventajas y desventajas

Ventajas - Funciona bien en alta dimensión (comparado con mallas tipo cuadratura). - Fácil de implementar. - Estimación directa del error (SE e intervalos).

Desventajas - Converge lento: \(O(n^{-1/2})\). - Puede tener varianza alta si \(f\) es muy dispersa. - Requiere técnicas (importance sampling, estratificación, antitéticas) para mejorar eficiencia.

10) Resumen en una sola línea

Para \(X\sim U(a,b)\):

\[ \int_a^b f(x)\,dx \approx (b-a)\frac{1}{n}\sum_{i=1}^n f(X_i) \]

n<-10
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

n<-50
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

n<-100
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

n<-500
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

n<-50
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

4*mean(ux^2+uy^2<=1)

## [1] 2.96

n<-100
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

4*mean(ux^2+uy^2<=1)

## [1] 3.32

n<-500
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

4*mean(ux^2+uy^2<=1)

## [1] 3.088

n<-1000
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

4*mean(ux^2+uy^2<=1)

## [1] 3.064

n<-10000
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
plot(ux,uy)

4*mean(ux^2+uy^2<=1)

## [1] 3.1576

n<-20000
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
4*mean(ux^2+uy^2<=1)

## [1] 3.1588

n<-1000000
ux<-runif(n,-1,1)
uy<-runif(n,-1,1)
4*mean(ux^2+uy^2<=1)

## [1] 3.1402

n<-seq(10,10000, by=1)
piest<-NULL
for(i in 1:length(n)){
  ux<-runif(n[i],-1,1)
uy<-runif(n[i],-1,1)
p<-4*mean(ux^2+uy^2<=1)
  piest[i]<-p
}
plot(n,piest,type="l")

plot(n, piest, type = "l",
     xlab = "Número de simulaciones (n)",
     ylab = "Estimación de π",
     main = "Convergencia Monte Carlo de π")

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

ggplot(data = data.frame(n = n, piest = piest),
       aes(x = n, y = piest)) +
  geom_line() +
  geom_hline(yintercept = pi, linetype = "dashed", color = "red") +
  labs(
    title = "Convergencia Monte Carlo de π",
    x = "Número de simulaciones (n)",
    y = "Estimación de π"
  )

x <- 35   # éxitos
n <- 50   # ensayos

# Modelo 1: prior más informativo
a1 <- 10
b1 <- 10

# Modelo 2: prior más difuso
a2 <- 1
b2 <- 1

post1 <- c(a1 + x, b1 + n - x)
post2 <- c(a2 + x, b2 + n - x)

mean1 <- post1[1] / sum(post1)
mean2 <- post2[1] / sum(post2)

mean1

## [1] 0.6428571

mean2

## [1] 0.6923077

set.seed(123)

M <- 1e6
p1 <- rbeta(M, post1[1], post1[2])
p2 <- rbeta(M, post2[1], post2[2])

mean(p1 > p2)

## [1] 0.279258

# Dos encuestas (dos grupos)
xA <- 330; nA <- 600   # Grupo A
xB <- 235; nB <- 500   # Grupo B

pA_hat <- xA/nA
pB_hat <- xB/nB

pA_hat

## [1] 0.55

pB_hat

## [1] 0.47

Previas no informativas

a0 <- 1
b0 <- 1

# Posteriores Beta-Binomial (conjugadas)
aA <- a0 + xA
bA <- b0 + (nA - xA)

aB <- a0 + xB
bB <- b0 + (nB - xB)

# ==========================
# 1) Media posterior (estimación bayesiana)
# ==========================
meanA <- aA / (aA + bA)
meanB <- aB / (aB + bB)

meanA

## [1] 0.5498339

meanB

## [1] 0.4701195

# ==========================
# 2) Intervalo creíble 95% para cada grupo
# ==========================
CI_A <- qbeta(c(0.025, 0.975), aA, bA)
CI_B <- qbeta(c(0.025, 0.975), aB, bB)

CI_A

## [1] 0.5099820 0.5893721

CI_B

## [1] 0.4266336 0.5138310

# ==========================
# 3) Probabilidad de que A sea mayor que B: P(pA > pB)
# ==========================
set.seed(123)
M <- 200000

pA <- rbeta(M, aA, bA)
pB <- rbeta(M, aB, bB)

prob_A_mayor <- mean(pA > pB)
prob_A_mayor

## [1] 0.99559

# (Opcional) Diferencia posterior y su intervalo creíble
d <- pA - pB
mean(d)

## [1] 0.07975193

quantile(d, c(0.025, 0.975))

##       2.5%      97.5% 
## 0.02036266 0.13839746

Monte Carlo

Jhonier Rangel

2026-01-20

Monte Carlo

Integración por Monte Carlo (MC) — Teoría y paso a paso

1) Idea central

2) Caso básico: integral en 1D sobre un intervalo

Estimador Monte Carlo

3) Paso a paso (algoritmo)

4) Propiedades importantes

4.1) Insesgadez

4.2) Varianza del estimador

4.3) Error típico (desviación estándar)

5) Intervalo de confianza (aprox. normal)

6) Generalización: integral como esperanza con una densidad cualquiera (Importance Sampling)

7) Integral en varias dimensiones (multi-D)

8) Método geométrico (hit-or-miss) — opcional

9) Ventajas y desventajas

10) Resumen en una sola línea

Previas no informativas