Monte Carlo

Introducción

El método de Monte Carlo es una técnica computacional que utiliza números aleatorios para resolver problemas matemáticos y estadísticos. Se emplea cuando los métodos analíticos son complejos o imposibles de aplicar.

Concepto

La idea principal es simular muchas veces un proceso aleatorio para aproximar resultados numéricos. A medida que aumenta el número de simulaciones, la precisión mejora gracias a la Ley de los Grandes Números.

Ejemplo: Aproximación de π

Podemos estimar el valor de π simulando puntos aleatorios dentro de un cuadrado y contando cuántos caen dentro de un círculo inscrito.

set.seed(123)

n <- 10000
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)

dentro <- x^2 + y^2 <= 1

pi_estimado <- 4 * mean(dentro)
pi_estimado

## [1] 3.1576

plot(x, y, col = ifelse(dentro, "blue", "red"),
     main = "Simulación Monte Carlo para π",
     xlab = "x", ylab = "y", pch = 20)

El método de Monte Carlo es una técnica numérica que permite aproximar integrales y resolver problemas probabilísticos mediante simulación aleatoria. Su fundamento teórico se basa en la Ley de los Grandes Números.

Integración Monte Carlo

Sea el problema de calcular una integral:

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

Si generamos variables aleatorias \(X_1, X_2, ..., X_n \sim U(a,b)\), entonces:

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

Por lo tanto, podemos estimar la integral como:

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

Demostración de Convergencia

Por la Ley de los Grandes Números (LGN):

\[ \frac{1}{n} \sum_{i=1}^n f(X_i) \xrightarrow{c.s.} \mathbb{E}[f(X)] \]

Entonces:

\[ \hat{I}_n \xrightarrow{c.s.} I \]

Esto demuestra que el estimador de Monte Carlo converge casi seguramente al valor real de la integral.

Relación con la Esperanza Matemática

La clave del método es que una integral puede interpretarse como una esperanza:

\[ \mathbb{E}[f(X)] = \int f(x) p(x) dx \]

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

\[ p(x) = \frac{1}{b-a} \]

Entonces:

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

Reordenando:

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

Error y Convergencia

El error del estimador es:

\[ \text{Error} \sim \mathcal{O}\left(\frac{1}{\sqrt{n}}\right) \]

Esto proviene del Teorema Central del Límite:

\[ \sqrt{n}(\hat{I}_n - I) \xrightarrow{d} \mathcal{N}(0, \sigma^2) \]

set.seed(123)

### ==============================
### 1. Distribución Uniforme
### ==============================
x <- runif(10000, 0, 1)
mean(x)              # Esperanza teórica = 0.5

## [1] 0.4975494

var(x)               # Varianza teórica = 1/12

## [1] 0.08219329

### ==============================
### 2. Distribución Normal
### ==============================
x <- rnorm(10000, mean = 0, sd = 1)
mean(x)              # ~0

## [1] 0.00199853

var(x)               # ~1

## [1] 1.003636

### ==============================
### 3. Distribución Exponencial
### ==============================
x <- rexp(10000, rate = 2)
mean(x)              # Esperanza = 1/rate = 0.5

## [1] 0.501162

### ==============================
### 4. Aproximación de una integral
### ∫_0^1 x^2 dx = 1/3
### ==============================
n <- 10000
x <- runif(n, 0, 1)
integral <- mean(x^2)
integral

## [1] 0.3354715

### ==============================
### 5. Integral más compleja
### ∫_0^1 sin(x) dx = 1 - cos(1)
### ==============================
x <- runif(n, 0, 1)
integral <- mean(sin(x))
integral

## [1] 0.4595658

1 - cos(1)

## [1] 0.4596977

### ==============================
### 6. Método Monte Carlo general
### ∫_a^b f(x) dx
### ==============================
f <- function(x) log(x + 1)
a <- 0
b <- 1

x <- runif(n, a, b)
integral <- (b - a) * mean(f(x))
integral

## [1] 0.3840189

### ==============================
### 7. Probabilidad P(X > 1) con Normal
### ==============================
x <- rnorm(10000)
prob <- mean(x > 1)
prob

## [1] 0.1616

### ==============================
### 8. Simulación Ley de los Grandes Números
### ==============================
x <- rnorm(10000)
promedios <- cumsum(x) / (1:10000)
tail(promedios, 1)

## [1] 0.0113895

### ==============================
### 9. Aproximación de π
### ==============================
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)
dentro <- x^2 + y^2 <= 1
pi_est <- 4 * mean(dentro)
pi_est

## [1] 3.1512

### ==============================
### 10. Valor esperado de función
### E[X^2] para X~N(0,1) = 1
### ==============================
x <- rnorm(n)
mean(x^2)

## [1] 1.005456

### ==============================
### 11. Varianza usando simulación
### ==============================
x <- rnorm(n)
mean(x^2) - mean(x)^2

## [1] 1.002023

### ==============================
### 12. Integral en 2D
### ∫∫ x*y dx dy en [0,1]x[0,1] = 1/4
### ==============================
x <- runif(n)
y <- runif(n)
mean(x * y)

## [1] 0.2445615

### ==============================
### 13. Distribución Binomial
### ==============================
x <- rbinom(10000, size = 10, prob = 0.5)
mean(x)   # Esperanza = np = 5

## [1] 5.0035

### ==============================
### 14. Distribución Poisson
### ==============================
x <- rpois(10000, lambda = 3)
mean(x)   # Esperanza = 3

## [1] 3.0245

### ==============================
### 15. Simulación de máximo
### ==============================
x <- runif(n)
max(x)

## [1] 0.9999805

Teorema de Bayes

La distribución posterior está dada por:

\[ p(\theta \mid x) = \frac{p(x \mid \theta)p(\theta)}{p(x)} \]

donde:

\(\theta\) es el parámetro
\(x\) son los datos
\(p(\theta)\) es la distribución a priori
\(p(x \mid \theta)\) es la verosimilitud

El término marginal es:

\[ p(x) = \int p(x \mid \theta)p(\theta)\, d\theta \]

Problema Computacional

Las cantidades de interés suelen ser esperanzas posteriores:

\[ \mathbb{E}[g(\theta) \mid x] = \int g(\theta)p(\theta \mid x)\, d\theta \]

Estas integrales rara vez tienen solución cerrada.

Método de Monte Carlo

Si podemos generar muestras \(\theta_1, \dots, \theta_n \sim p(\theta \mid x)\), entonces:

\[ \mathbb{E}[g(\theta) \mid x] \approx \frac{1}{n} \sum_{i=1}^n g(\theta_i) \]

Justificación: Ley de los Grandes Números

Por la Ley de los Grandes Números:

\[ \frac{1}{n} \sum_{i=1}^n g(\theta_i) \xrightarrow{c.s.} \mathbb{E}[g(\theta) \mid x] \]

cuando \(n \to \infty\).

Error y Convergencia

Por el Teorema Central del Límite:

\[ \sqrt{n} \left( \hat{\mu}_n - \mu \right) \xrightarrow{d} \mathcal{N}(0, \sigma^2) \]

Por lo tanto, el error es del orden:

\[ \mathcal{O}\left(\frac{1}{\sqrt{n}}\right) \]

Métodos Monte Carlo Bayesianos

Cuando no se puede muestrear directamente de \(p(\theta \mid x)\), se utilizan métodos MCMC:

Metropolis-Hastings
Gibbs Sampling

Estos generan una cadena de Markov cuya distribución estacionaria es la posterior.

Ejemplo

Supongamos:

\[ X \sim \text{Binomial}(n, \theta), \quad \theta \sim \text{Beta}(a,b) \]

Entonces la posterior es:

\[ \theta \mid X \sim \text{Beta}(a + x, b + n - x) \]

############################################################
# TEORÍA: MONTE CARLO Y SU RELACIÓN CON ESTADÍSTICA BAYESIANA
############################################################

# 1. PROBLEMA BAYESIANO

# En estadística bayesiana queremos calcular la distribución posterior:

# p(θ | x) = (p(x | θ) * p(θ)) / p(x)

# donde:
# θ = parámetro desconocido
# x = datos observados
# p(θ) = prior
# p(x | θ) = verosimilitud
# p(θ | x) = posterior

# El problema es que:

# p(x) = ∫ p(x | θ)p(θ) dθ

# suele ser difícil de calcular analíticamente.

############################################################
# 2. IDEA DE MONTE CARLO EN BAYES

# En lugar de resolver la integral, usamos simulación:

# Si θ_1, θ_2, ..., θ_n ~ p(θ | x)

# entonces cualquier esperanza posterior se aproxima por:

# E[g(θ) | x] ≈ (1/n) Σ g(θ_i)

############################################################
# 3. RELACIÓN CON INTEGRALES

# En Bayes queremos calcular:

# E[g(θ) | x] = ∫ g(θ) p(θ | x) dθ

# Monte Carlo lo aproxima como:

# (1/n) Σ g(θ_i)

# Esto es una aplicación directa de la Ley de los Grandes Números.

############################################################
# 4. DEMOSTRACIÓN DE CONVERGENCIA

# Si θ_i ~ p(θ | x), entonces:

# (1/n) Σ g(θ_i) → E[g(θ) | x]

# cuando n → ∞

# por la Ley de los Grandes Números.

############################################################
# 5. MÉTODOS IMPORTANTES

# (a) Monte Carlo directo:
# cuando podemos muestrear directamente de la posterior.

# (b) MCMC (Markov Chain Monte Carlo):
# cuando no podemos muestrear directamente.

# Ejemplos:
# - Metropolis-Hastings
# - Gibbs Sampling

############################################################
# 6. EJEMPLO CONJUGADO (ANALÍTICO + MONTE CARLO)

# Supongamos:
# X ~ Binomial(n, θ)
# θ ~ Beta(a, b)

# Entonces:
# θ | X ~ Beta(a + x, b + n - x)

# Simulamos de la posterior:

set.seed(123)

n <- 10
x_obs <- 7
a <- 2
b <- 2

theta_post <- rbeta(10000, a + x_obs, b + n - x_obs)

# Estimación de la esperanza posterior:
mean(theta_post)

## [1] 0.6429457

############################################################
# 7. INTERPRETACIÓN

# Monte Carlo permite:
# - Aproximar integrales bayesianas
# - Estimar parámetros posteriores
# - Resolver modelos donde no hay solución cerrada

############################################################
# 8. ERROR

# Error Monte Carlo:

# Error ≈ O(1/sqrt(n))

# Más simulaciones → mayor precisión

Queremos aproximar:

\[ \mathbb{E}[g(\theta) \mid x] = \int g(\theta)p(\theta \mid x)\, d\theta \]

usando simulación Monte Carlo y MCMC.

1. Convergencia Monte Carlo (Ley de los Grandes Números)

\[ \hat{\mu}_n = \frac{1}{n} \sum_{i=1}^n g(\theta_i) \]

\[ \hat{\mu}_n \xrightarrow{c.s.} \mathbb{E}[g(\theta)] \]

Código + gráfica

set.seed(123)

n <- 5000
x <- rnorm(n)

promedios <- cumsum(x) / (1:n)

plot(promedios, type = "l", col = "blue",
     main = "Convergencia (Ley de los Grandes Números)",
     ylab = "Media acumulada", xlab = "Iteraciones")

abline(h = 0, col = "red", lwd = 2)

n_vals <- seq(10, 5000, by = 10)
errores <- sapply(n_vals, function(n) {
  x <- rnorm(n)
  abs(mean(x))
})

plot(n_vals, errores, type = "l", col = "darkgreen",
     main = "Error Monte Carlo",
     xlab = "n", ylab = "Error")

set.seed(123)

n <- 5000
theta <- numeric(n)
theta[1] <- 0

for (i in 2:n) {
  propuesta <- rnorm(1, mean = theta[i-1], sd = 1)
  
  alpha <- min(1, exp(-(propuesta^2 - theta[i-1]^2)/2))
  
  if (runif(1) < alpha) {
    theta[i] <- propuesta
  } else {
    theta[i] <- theta[i-1]
  }
}
plot(theta, type = "l", col = "blue",
     main = "Trace plot MCMC",
     xlab = "Iteraciones", ylab = expression(theta))

hist(theta, probability = TRUE, col = "lightblue",
     main = "Distribución MCMC vs Normal")

curve(dnorm(x, 0, 1), col = "red", lwd = 2, add = TRUE)

media_acum <- cumsum(theta) / (1:n)

plot(media_acum, type = "l", col = "purple",
     main = "Convergencia de la media (MCMC)",
     xlab = "Iteraciones", ylab = "Media")

abline(h = 0, col = "red", lwd = 2)

acf(theta, main = "Autocorrelación MCMC")

Demostración de Mínimos Cuadrados en Regresión Lineal Simple y Estadístico t de Student

1. Modelo de Regresión Lineal Simple

El modelo de regresión lineal simple se define como:

[ Y_i = _0 + _1 X_i + _i, i = 1,2,,n]

donde:

(Y_i): variable respuesta
(X_i): variable explicativa
(_0, _1): parámetros desconocidos
(_i): término de error, con (E(_i)=0), (Var(_i)=^2)

2. Método de Mínimos Cuadrados

El objetivo es estimar (_0) y (_1) minimizando:

[ S(_0,1) = {i=1}^{n} (Y_i - _0 - _1 X_i)^2]

Derivación

Derivamos respecto a (_0) y (_1) e igualamos a cero:

[ = -2 (Y_i - _0 - _1 X_i) = 0]

[ = -2 X_i (Y_i - _0 - _1 X_i) = 0]

Simplificando:

[ Y_i = n_0 + _1 X_i]

[ X_i Y_i = _0 X_i + _1 X_i^2]

3. Solución de los Estimadores

Definimos medias:

[ {X} = X_i, {Y} = Y_i]

Pendiente:

[ _1 = ]

Intercepto:

[ _0 = {Y} - _1 {X}]

4. Varianza del Error

El estimador de la varianza es:

[ ^2 = (Y_i - _i)^2]

5. Varianza del Estimador de la Pendiente

[ Var(_1) = ]

Estimación:

[ (_1) = ]

6. Estadístico t de Student para la Pendiente

Queremos probar:

[ H_0: _1 = 0 H_1: _1 ]

Estadístico:

[ t = ]

Bajo (H_0):

[ t = t*{n-2}]

7. Demostración de la Distribución t

Sabemos que:

[ _1 N(_1, )]

Entonces:

[ Z = N(0,1)]

Además:

[ ^2_{n-2}]

y es independiente de (_1).

Por definición:

[ t = t_{n-2}]

# Instalar paquetes si no los tienes
# install.packages("ggplot2")

library(ggplot2)

# Crear datos de ejemplo
set.seed(123)
x <- 1:50
y <- 2 + 0.5*x + rnorm(50, mean = 0, sd = 3)

data <- data.frame(x, y)

# Modelo de regresión lineal simple
modelo <- lm(y ~ x, data = data)

# Resumen del modelo
summary(modelo)

## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0336 -1.8667 -0.2458  1.9977  6.4790 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.2071     0.8056   2.739  0.00861 ** 
## x             0.4959     0.0275  18.036  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.806 on 48 degrees of freedom
## Multiple R-squared:  0.8714, Adjusted R-squared:  0.8687 
## F-statistic: 325.3 on 1 and 48 DF,  p-value: < 2.2e-16

# Gráfico con ggplot2
ggplot(data, aes(x = x, y = y)) +
  geom_point(color = "blue") +                 # puntos
  geom_smooth(method = "lm", color = "red") +  # recta de regresión
  labs(title = "Regresión Lineal Simple",
       x = "Variable X",
       y = "Variable Y") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Regresión Lineal Simple: Derivación por Mínimos Cuadrados

1. Modelo

El modelo de regresión lineal simple es:

\[ Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i \]

donde:

\(Y_i\): variable dependiente
\(X_i\): variable independiente
\(\beta_0, \beta_1\): parámetros
\(\varepsilon_i\): error aleatorio

2. Función de Mínimos Cuadrados

Se desea minimizar:

\[ S(\beta_0, \beta_1) = \sum_{i=1}^{n} (Y_i - \beta_0 - \beta_1 X_i)^2 \]

3. Derivadas Parciales

Derivamos respecto a \(\beta_0\):

\[ \frac{\partial S}{\partial \beta_0} = -2 \sum (Y_i - \beta_0 - \beta_1 X_i) \]

Derivamos respecto a \(\beta_1\):

\[ \frac{\partial S}{\partial \beta_1} = -2 \sum X_i (Y_i - \beta_0 - \beta_1 X_i) \]

4. Ecuaciones Normales

Igualando a cero:

\[ \sum Y_i = n\beta_0 + \beta_1 \sum X_i \]

\[ \sum X_i Y_i = \beta_0 \sum X_i + \beta_1 \sum X_i^2 \]

5. Solución de los Estimadores

Definimos:

\[ \bar{X} = \frac{1}{n}\sum X_i, \quad \bar{Y} = \frac{1}{n}\sum Y_i \]

Pendiente

\[ \hat{\beta}_1 = \frac{\sum (X_i - \bar{X})(Y_i - \bar{Y})}{\sum (X_i - \bar{X})^2} \]

Intercepto

\[ \hat{\beta}_0 = \bar{Y} - \hat{\beta}_1 \bar{X} \]

6. Forma Alternativa de la Pendiente

\[ \hat{\beta}_1 = \frac{\sum X_i Y_i - n\bar{X}\bar{Y}}{\sum X_i^2 - n\bar{X}^2} \]

7. Interpretación

\(\hat{\beta}_1\): cambio promedio en \(Y\) por unidad de \(X\)
\(\hat{\beta}_0\): valor esperado de \(Y\) cuando \(X=0\)

# Instalar paquetes si es necesario
# install.packages("ggplot2")

library(ggplot2)

# ---------------------------
# 1. Datos simulados
# ---------------------------
set.seed(123)
x <- 1:50
y <- 5 + 2*x + rnorm(50, mean = 0, sd = 5)

data <- data.frame(x, y)

# ---------------------------
# 2. Modelo de regresión
# ---------------------------
modelo <- lm(y ~ x, data = data)

# ---------------------------
# 3. Inferencia
# ---------------------------
summary(modelo)   # incluye t-test para la pendiente

## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0560  -3.1111  -0.4097   3.3295  10.7983 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.34508    1.34274   3.981 0.000232 ***
## x            1.99321    0.04583  43.494  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.676 on 48 degrees of freedom
## Multiple R-squared:  0.9753, Adjusted R-squared:  0.9747 
## F-statistic:  1892 on 1 and 48 DF,  p-value: < 2.2e-16

# Intervalos de confianza
confint(modelo)

##                2.5 %   97.5 %
## (Intercept) 2.645329 8.044835
## x           1.901072 2.085354

# ---------------------------
# 4. Gráfico con ggplot2
# ---------------------------
ggplot(data, aes(x = x, y = y)) +
  geom_point(color = "blue") +
  geom_smooth(method = "lm", 
              se = TRUE,        # banda de confianza
              color = "red") +
  labs(title = "Regresión Lineal Simple con Inferencia",
       subtitle = "Incluye intervalo de confianza",
       x = "Variable X",
       y = "Variable Y") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Regresión Lineal Múltiple: Formulación Matricial y Derivación por Mínimos Cuadrados

1. Modelo de Regresión Lineal Múltiple

El modelo se define como:

\[ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_k X_{ik} + \varepsilon_i \]

para \(i = 1, \dots, n\), donde:

\(Y_i\): variable respuesta
\(X_{ij}\): variables explicativas
\(\beta_j\): parámetros desconocidos
\(\varepsilon_i \sim N(0, \sigma^2)\), independientes

2. Forma Matricial

Se puede escribir como:

\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

donde:

\[ \mathbf{Y} = \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{pmatrix}, \quad \mathbf{X} = \begin{pmatrix} 1 & X_{11} & \cdots & X_{1k} \\ 1 & X_{21} & \cdots & X_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & X_{n1} & \cdots & X_{nk} \end{pmatrix} \]

\[ \boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{pmatrix}, \quad \boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{pmatrix} \]

3. Supuestos del Modelo

\(E(\boldsymbol{\varepsilon}) = 0\)
\(Var(\boldsymbol{\varepsilon}) = \sigma^2 I_n\)
Independencia de los errores
Normalidad: \(\boldsymbol{\varepsilon} \sim N(0, \sigma^2 I)\)

4. Método de Mínimos Cuadrados

Se busca minimizar:

\[ S(\boldsymbol{\beta}) = (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})'(\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}) \]

5. Derivación del Estimador

Expandimos:

\[ S(\boldsymbol{\beta}) = \mathbf{Y}'\mathbf{Y} - 2\boldsymbol{\beta}'\mathbf{X}'\mathbf{Y} + \boldsymbol{\beta}'\mathbf{X}'\mathbf{X}\boldsymbol{\beta} \]

Derivamos respecto a \(\boldsymbol{\beta}\):

\[ \frac{\partial S}{\partial \boldsymbol{\beta}} = -2\mathbf{X}'\mathbf{Y} + 2\mathbf{X}'\mathbf{X}\boldsymbol{\beta} \]

Igualando a cero:

\[ \mathbf{X}'\mathbf{X}\boldsymbol{\beta} = \mathbf{X}'\mathbf{Y} \]

6. Solución (Ecuaciones Normales)

Si \(\mathbf{X}'\mathbf{X}\) es invertible:

\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1} \mathbf{X}'\mathbf{Y} \]

7. Valores Ajustados y Residuos

\[ \hat{\mathbf{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{H}\mathbf{Y} \]

donde:

\[ \mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' \]

Residuos:

\[ \mathbf{e} = \mathbf{Y} - \hat{\mathbf{Y}} \]

8. Propiedades del Estimador

Esperanza

\[ E(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta} \]

Varianza

\[ Var(\hat{\boldsymbol{\beta}}) = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1} \]

9. Estimación de la Varianza del Error

\[ \hat{\sigma}^2 = \frac{\mathbf{e}'\mathbf{e}}{n - k - 1} \]

10. Inferencia: Estadístico t

Para cada coeficiente \(\beta_j\):

\[ t_j = \frac{\hat{\beta}_j}{\sqrt{\hat{\sigma}^2 c_{jj}}} \]

donde \(c_{jj}\) es el elemento \(j\)-ésimo de \((\mathbf{X}'\mathbf{X})^{-1}\).

\[ t_j \sim t_{n-k-1} \]

11. Inferencia Global: Estadístico F

\[ F = \frac{(SSR / k)}{(SSE / (n-k-1))} \]

donde:

SSR: suma de cuadrados de regresión
SSE: suma de cuadrados del error

\[ F \sim F_{k, n-k-1} \]

12. Interpretación

Cada \(\beta_j\) mide el efecto marginal de \(X_j\) sobre \(Y\), manteniendo constantes las demás variables.
El modelo generaliza la regresión simple.
La forma matricial simplifica cálculos.

13. Condición de Existencia

\[ \text{rank}(\mathbf{X}) = k+1 \]

No debe haber multicolinealidad perfecta.

# Instalar paquetes si es necesario
# install.packages("ggplot2")

library(ggplot2)

# ---------------------------
# 1. Datos simulados
# ---------------------------
set.seed(123)

n <- 100
x1 <- rnorm(n, 10, 2)
x2 <- rnorm(n, 5, 1)

# Modelo real: Y = 3 + 2*x1 - 1*x2 + error
y <- 3 + 2*x1 - 1*x2 + rnorm(n, 0, 3)

data <- data.frame(x1, x2, y)

# ---------------------------
# 2. Modelo de regresión múltiple
# ---------------------------
modelo <- lm(y ~ x1 + x2, data = data)

# ---------------------------
# 3. Inferencia
# ---------------------------
summary(modelo)     # t-test para cada coeficiente

## 
## Call:
## lm(formula = y ~ x1 + x2, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6190 -1.9822 -0.3733  1.8641  6.2395 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0456     2.2334   2.259  0.02611 *  
## x1            1.8002     0.1573  11.444  < 2e-16 ***
## x2           -0.9286     0.2970  -3.127  0.00233 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.854 on 97 degrees of freedom
## Multiple R-squared:  0.5986, Adjusted R-squared:  0.5903 
## F-statistic: 72.32 on 2 and 97 DF,  p-value: < 2.2e-16

confint(modelo)     # intervalos de confianza

##                 2.5 %    97.5 %
## (Intercept)  0.613011  9.478189
## x1           1.488037  2.112448
## x2          -1.517997 -0.339135

# ---------------------------
# 4. Predicciones
# ---------------------------
data$y_pred <- predict(modelo)

# ---------------------------
# 5. Gráfico (relación parcial con x1)
# ---------------------------
ggplot(data, aes(x = x1, y = y)) +
  geom_point(color = "blue") +
  geom_line(aes(y = y_pred), color = "red") +
  labs(title = "Regresión Lineal Múltiple",
       subtitle = "Relación de Y con X1 (controlando X2)",
       x = "X1",
       y = "Y") +
  theme_minimal()

```rmarkdown id=“supreglm” — title: “Supuestos de la Regresión Lineal Múltiple” output: html_document —

Supuestos de la Regresión Lineal Múltiple

Para que los estimadores de mínimos cuadrados sean válidos (insesgados, eficientes y permitan inferencia), el modelo debe cumplir ciertos supuestos.

1. Linealidad

El modelo es lineal en los parámetros:

\[ Y = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k + \varepsilon \]

Esto implica que la relación entre las variables explicativas y la respuesta es lineal en los coeficientes.

2. Independencia de los errores

\[ Cov(\varepsilon_i, \varepsilon_j) = 0 \quad \text{para } i \neq j \]

Los errores no deben estar correlacionados entre sí.

Esto es importante especialmente en: - series de tiempo
- datos longitudinales

3. Homocedasticidad

La varianza de los errores es constante:

\[ Var(\varepsilon_i) = \sigma^2 \]

Esto significa que la dispersión de los errores no depende de los valores de \(X\).

4. Normalidad de los errores

\[ \varepsilon_i \sim N(0, \sigma^2) \]

Este supuesto es necesario para realizar inferencia: - pruebas t
- intervalos de confianza
- pruebas F

5. No multicolinealidad perfecta

Las variables explicativas no deben ser combinaciones lineales exactas:

\[ \text{rank}(X) = k+1 \]

Si existe multicolinealidad perfecta: - no se puede invertir \(X'X\)
- no existe solución única

6. Exogeneidad (media cero)

\[ E(\varepsilon_i \mid X) = 0 \]

Implica que: - no hay variables omitidas relevantes
- no hay correlación entre errores y regresores

7. Forma matricial de los supuestos

\[ E(\boldsymbol{\varepsilon}) = 0 \]

\[ Var(\boldsymbol{\varepsilon}) = \sigma^2 I \]

8. Consecuencias de violar los supuestos

No linealidad → modelo mal especificado
Heterocedasticidad → errores estándar incorrectos
Autocorrelación → inferencia inválida
No normalidad → afecta pruebas en muestras pequeñas
Multicolinealidad → estimaciones inestables
Endogeneidad → estimadores sesgados

# =========================
# 1. Paquetes
# =========================
# install.packages(c("ggplot2","lmtest","car"))

library(ggplot2)
library(lmtest)

## Cargando paquete requerido: zoo

## 
## Adjuntando el paquete: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(car)

## Cargando paquete requerido: carData

# =========================
# 2. Datos simulados
# =========================
set.seed(123)

n <- 100
x1 <- rnorm(n, 10, 2)
x2 <- rnorm(n, 5, 1)

y <- 3 + 2*x1 - 1*x2 + rnorm(n, 0, 3)

data <- data.frame(x1, x2, y)

# =========================
# 3. Modelo
# =========================
modelo <- lm(y ~ x1 + x2, data = data)

summary(modelo)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6190 -1.9822 -0.3733  1.8641  6.2395 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0456     2.2334   2.259  0.02611 *  
## x1            1.8002     0.1573  11.444  < 2e-16 ***
## x2           -0.9286     0.2970  -3.127  0.00233 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.854 on 97 degrees of freedom
## Multiple R-squared:  0.5986, Adjusted R-squared:  0.5903 
## F-statistic: 72.32 on 2 and 97 DF,  p-value: < 2.2e-16

# =========================
# 4. TEST DE SUPUESTOS
# =========================

# 4.1 Normalidad (Shapiro-Wilk)
shapiro.test(residuals(modelo))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo)
## W = 0.97815, p-value = 0.09534

# 4.2 Heterocedasticidad (Breusch-Pagan)
bptest(modelo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 0.039949, df = 2, p-value = 0.9802

# 4.3 Autocorrelación (Durbin-Watson)
dwtest(modelo)

## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.8296, p-value = 0.1989
## alternative hypothesis: true autocorrelation is greater than 0

# 4.4 Multicolinealidad (VIF)
vif(modelo)

##       x1       x2 
## 1.002459 1.002459

# =========================
# 5. GRÁFICOS DE DIAGNÓSTICO
# =========================

# Base R (muy usados en econometría)
par(mfrow = c(2,2))
plot(modelo)

# =========================
# 6. GRÁFICOS CON ggplot2
# =========================

# Residuos vs ajustados
data$residuos <- residuals(modelo)
data$ajustados <- fitted(modelo)

ggplot(data, aes(x = ajustados, y = residuos)) +
  geom_point(color = "blue") +
  geom_hline(yintercept = 0, color = "red") +
  labs(title = "Residuos vs Ajustados",
       x = "Valores ajustados",
       y = "Residuos") +
  theme_minimal()

# QQ-plot normalidad
ggplot(data, aes(sample = residuos)) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "QQ-Plot de Residuos") +
  theme_minimal()

Monte Carlo y Mínimos cuadrados

Jhonier Rangel

2026-03-19

Monte Carlo

Introducción

Concepto

Ejemplo: Aproximación de π

Integración Monte Carlo

Demostración de Convergencia

Relación con la Esperanza Matemática

Error y Convergencia

Teorema de Bayes

Problema Computacional

Método de Monte Carlo

Justificación: Ley de los Grandes Números

Error y Convergencia

Métodos Monte Carlo Bayesianos

Ejemplo

1. Convergencia Monte Carlo (Ley de los Grandes Números)

Código + gráfica

Demostración de Mínimos Cuadrados en Regresión Lineal Simple y Estadístico t de Student

1. Modelo de Regresión Lineal Simple

2. Método de Mínimos Cuadrados

Derivación

3. Solución de los Estimadores

Pendiente:

Intercepto:

4. Varianza del Error

5. Varianza del Estimador de la Pendiente

6. Estadístico t de Student para la Pendiente

Estadístico:

7. Demostración de la Distribución t

Regresión Lineal Simple: Derivación por Mínimos Cuadrados

1. Modelo

2. Función de Mínimos Cuadrados

3. Derivadas Parciales

4. Ecuaciones Normales

5. Solución de los Estimadores

Pendiente

Intercepto

6. Forma Alternativa de la Pendiente

7. Interpretación

Regresión Lineal Múltiple: Formulación Matricial y Derivación por Mínimos Cuadrados

1. Modelo de Regresión Lineal Múltiple

2. Forma Matricial

3. Supuestos del Modelo

4. Método de Mínimos Cuadrados

5. Derivación del Estimador

6. Solución (Ecuaciones Normales)

7. Valores Ajustados y Residuos

8. Propiedades del Estimador

Esperanza

Varianza

9. Estimación de la Varianza del Error

10. Inferencia: Estadístico t

11. Inferencia Global: Estadístico F

12. Interpretación

13. Condición de Existencia

Supuestos de la Regresión Lineal Múltiple

1. Linealidad

2. Independencia de los errores

3. Homocedasticidad

4. Normalidad de los errores

5. No multicolinealidad perfecta

6. Exogeneidad (media cero)

7. Forma matricial de los supuestos

8. Consecuencias de violar los supuestos