Un modelo estadístico \(F\) es un conjuto de distribuciones (o densidades o funciones de regresión). Un modelo paramétrico es un conjunto \(F\) que puede ser parametrizado por un número finito de parámetros. Por ejemplo, si suponemos que los datos provienen de una distribución Normal, el modelo es:

\[F=\bigg\{p(x;\mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}}exp\bigg(-\frac{1}{2\sigma^2}(x-\mu)^2\bigg), \mu \in \mathbb{R}, \sigma>0\bigg\}\]

En general, un modelo paramétrico tiene la forma

\[F=\bigg\{p(x;\theta):\theta \in \Theta \bigg\}\]

donde \(\theta\) es un parámetro desconocido (o un vector de parámetros) que puede tomar valores en el espacio paramétrico \(\Theta\).

Un modelo no paramétrico es un conjunto \(F\) que no puede ser parametrizado por un número finito de parámetros.

Variables aleatorias

Recordemos brevemente algunos conceptos de variables aleatorias.

Una variable aleatoria es un mapeo \[X:\Omega \to \mathbb{R}\] que asigna un número real \(X(\omega)\) a cada elemento \(\omega\) en el espacio de resultados.


Ejemplo. Lanzamos una moneda justa dos veces, definimos \(X\) como en el número de soles, entonces la variable aleatoria se pueden resumir como:

\(\omega\) \(P(\{\omega\})\) \(X(\omega)\)
AA 1/4 0
AS 1/4 1
SA 1/4 1
SS 1/4 2

y su distribución:

x \(P(X=x)\)
0 1/4
1 1/2
2 1/4
La función de distribución acumulada es la función \(P_X:\mathbb{R}\to[0,1]\) definida como: \[P_X(x)=P(X\leq x)\]


En el ejemplo: \[ p_X(x) = \left\{ \begin{array}{lr} 0 & x < 0\\ 1/4 & 0 \leq x < 1 \\ 3/4 & 1 \leq x < 2 \\ 1 & x \ge 2 \end{array} \right. \]

Una variable aleatoria \(X\) es discreta si toma un número contable de valores \(\{x_1,x_2,...\}\). En este caso definimos la función de probabilidad o la función masa de probabilidad de X como \(p_X(x)=P(X=x)\).


Notemos que \(p_X(x)\geq 0\) para toda \(x \in \mathbb{R}\) y \(\sum_i p_X(x)=1\). Más aún, la función de distribución acumulada esta relacionada con \(p_X\) por \[P_X(x)=P(X \leq x)= \sum_{x_i\leq x} = \sum_{x_i\leq x}p_{X}(x_i)\]

Sea \(X\) una variable aleatoria con FDA \(P_X\). La función de distribución acumulada inversa o función de cuantiles se define como: \[P_X^{-1}(q) = inf\{x:P_X(x)>q\}\] para \(q \in [0,1]\).


Llamamos a \(P^{-1}(1/4)\) el primer cuartil, a \(P^{-1}(1/2)\) la mediana y \(P^{-1}(3/4)\) el tercer cuartil.

Familias discretas importantes

Distribución Uniforme

Sea \(k>1\) un entero dado. Decimos que \(X\) tiene una distribución uniforme en \(\{1,...,k\}\) si tiene una función de probailidaddada por:

\[ p(x) = \left\{ \begin{array}{lr} 1/k & x \in \{1,...,k\}\\ 0 & e.o.c. \\ \end{array} \right. \]

\(E(X) = (a+b)/2, Var(X)=(n^2-1)/12\)

Distribución Bernoulli

Sea \(X\) la variable aleatoria que representa un lanzamiento de moneda, entonces \(P(X=1)=p\) y \(P(X=0)=1-p\) para alguna \(p\in[0,1]\). Decimos que \(X\) tiene una distribución Bernoulli (\(X \sim Bernoulli(p)\)), y su función de distribución es:

\[ p(x) = \left\{ \begin{array}{lr} p^x(1-p)^{1-x} & x \in \{0,1\}\\ 0 & e.o.c. \\ \end{array} \right. \]

\(E(X) = p, Var(X)=p(1-p)\)

Distribución Binomial

Supongamos que tenemos una moneda que cae en sol con probabilidad \(p\), para alguna \(0 \leq p \leq 1\). Lanzamos la moneda \(n\) veces y sea \(X\) el número de soles. Suponemos que los lanzamientos son independientes, entonces la función de distribución es:

\[ p(x) = \left\{ \begin{array}{lr} {n \choose x}p^x(1-p)^{n-x} & x \in \{0,1,...,n\}\\ 0 & e.o.c. \\ \end{array} \right. \]

\(E(X) = np, Var(X)=np(1-p)\)
Si \(X_1 \sim Binomial(n_1, p)\) y \(X_2 \sim Binomial(n_2,p)\) entonces \(X_1 + X_2 \sim Binomial(n_1+n_2, p)\). En general la distribución binomial describe el comportamiento de una variable \(X\) que cuenta número de éxitos tal que: 1) el número de observaciones \(n\) esta fijo, 2) cada observación es independiente, 3) cada observación representa uno de dos posibles eventos (éxito o fracaso) y 3) la probabilidad de éxito \(p\) es la misma en cada observación.

xy <- data.frame(x = 0:20)
ggplot(xy) +
  geom_bar(aes(x = x, y = dbinom(x, size = 20, prob = 0.5)), stat = "identity")

xy <- data.frame(x = 0:20)
ggplot(xy) +
  geom_bar(aes(x = x, y = dbinom(x, size = 20, prob = 0.1)), stat = "identity")



Distribución Poisson

\(X\) tienen una distribución Poisson con prámetro \(\lambda\) si \[ p(x) = \left\{ \begin{array}{lr} e^{-\lambda} \frac{\lambda^x}{x!} & x \in \{0,1,...\}\\ 0 & e.o.c. \\ \end{array} \right. \]

\(E(X) = \lambda, Var(X)=\lambda\)


La distribución Poisson se utiliza con frecuencia para modelar conteos de eventos raros, por ejemplo número de accidentes de tráfico. La distribución Poisson es un caso límite de la distribución binomial cuando el número de casos es muy grande y la probabilidad de éxito \(p\) es chica. Si \(X_1 \sim Poisson(\lambda_1)\) y \(X_2 \sim Poisson(\lambda_2)\) entonces \(X_1 + X_2 \sim Poisson(\lambda_1 + \lambda_2)\).

xy <- data.frame(x = 0:20)
ggplot(xy) +
  geom_bar(aes(x = x, y = dpois(x, lambda = 2)), stat = "identity")


#### Distribución geométrica \(X\) tiene distribución geométrica con parámetro \(p \in (0,1)\), \(X \sim Geom(p)\) si, \[ p(x) = \left\{ \begin{array}{lr} p(1-p)^{k-1} & x \in \{1,2,...\}\\ 0 & e.o.c. \\ \end{array} \right. \]

\(E(X)=1/p, Var(X)=(1-p)/p^2\)
con \(k \geq 1\). Podemos pensar en \(X\) como el número de lanzamientos necesarios hasta que obtenemos el primer sol en los lanzamientos de una moneda.

xy <- data.frame(x = 0:12)
ggplot(xy) +
  geom_bar(aes(x = x, y = dgeom(x, prob = .5)), stat = "identity")



Variables aleatorias continuas

Una variable aleatoria \(X\) es continua si existe una función \(p_x\) tal que \(p_X(x) \geq 0\) para toda \(x\), \(\int_{-\infty}^{\infty}p_X(x)dx=1\) y para toda \(a\leq b\),

\[P(a < X < b) = \int_{a}^b p_X(x)dx\]

La función \(p_X(x)\) se llama la función de densidad de probabilidad (fdp). Tenemos que \[P_X(x)=\int_{-\infty}^x p_X(t)dt\] y \(p_X(x)=P_X^{\'}(x)\) en todos los puntos \(x\) en los que la FDA \(P_X\) es diferenciable.


Ejemplo. Supongamos que elegimos un número al azar entre cero y uno, entonces

\[ p(x) = \left\{ \begin{array}{lr} \frac{1}{b-a} & x \in [0, 1]\\ 0 & e.o.c. \\ \end{array} \right. \] es claro que \(p_X(x) \geq 0\) para toda \(x\) y \(\int_{-\infty}^{\infty}p_X(x)dx=1\), la FDA esta dada por

\[ P_X(x) = \left\{ \begin{array}{lr} 0 & x < 0 \\ x & x \in [0,1]\\ 1 & x>b \\ \end{array} \right. \]

Vale la pena notar que en el caso de variables aleatorias continuas \(P(X=x)=0\) para toda \(x\) y pensar en \(p_X(x)\) como \(P(X=x)\) solo tiene sentido en el caso discreto.

Familias Continuas importantes

Distribución Uniforme

\(X\) tiene una distribución \(Uniforme(a,b)\) si

\[ p(x) = \left\{ \begin{array}{lr} \frac{1}{b-a} & x \in [a,b]\\ 0 & e.o.c. \\ \end{array} \right. \]

donde \(a < b\). La función de distribución es

\[ P_X(x) = \left\{ \begin{array}{lr} 0 & x < a \\ \frac{x-a}{b-a} & x \in [a,b]\\ 1 & x>b \\ \end{array} \right. \]

\(E(X) = (a+b)/2, Var(X)= (b-a)^2/12\)


Distribución Normal

\(X\) tiene una distribución normal con parámetros \(\mu\) y \(\sigma\), denotado \(X\sim N(\mu, \sigma^2)\) si

\[p(x) = \frac{1}{\sigma\sqrt{2\pi}}exp\bigg(-\frac{1}{2\sigma^2}(x-\mu)^2\bigg)\]

\(E(X)=\mu, Var(X)=\sigma^2\)


donde \(\mu \in \mathbb{R}\) y \(\sigma>0\).

Decimos que \(X\) tiene una distribución Normal estándar si \(\mu=0\) y \(\sigma=1\). Una variable aleatoria Normal estándar se denota tradicionalmente por \(Z\), su función de densidad de probabilidad por \(\phi(z)\) y la función de probabilidad acumulada por \(\Phi(z)\).

Algunas porpiedades importantes son:

  1. Si \(X \sim N(\mu, \sigma^2)\), entonces \(Z=(X-\mu)/\sigma \sim N(0,1)\).

  2. Si \(Z \sim N(0, 1)\) entonces \(X = \mu + \sigma Z \sim N(\mu, sigma^2)\).

  3. Si \(X_i \sim N(\mu_i, \sigma_i^2)\), \(i=1,...,n\) independientes, entonces: \[\sum_{i=1}^n X_i \sim N(\sum_{i=1}^n \mu_i, \sum_{i=1}^n \sigma_i^2)\]

  4. Se sigue de 1 que si \(X\sim N(\mu, \sigma^2)\), entonces \[P(a<X<b) = P\big(\frac{a-\mu}{\sigma} < Z < \frac{b-\mu}{\sigma}\big)= \Phi\big(\frac{b-\mu}{\sigma}\big) - \Phi\big(\frac{a-\mu}{\sigma}\big)\]

xy <- data.frame(x = -5:5, y = seq(0, 1, length.out = 11))
ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = dnorm, color = "orange") + 
  stat_function(fun = dnorm, args = list(mean = 1), color = "steelblue") +
  stat_function(fun = dnorm, args = list(sd = 2), color = "darkgray") +
  labs(y = "", title = "fdp's normales")

ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = pnorm, color = "orange") + 
  stat_function(fun = pnorm, args = list(mean = 1), color = "steelblue") +
  stat_function(fun = pnorm, args = list(sd = 2), color = "darkgray") +
  labs(y = "", title = "FDA's normales")

xy <- data.frame(x = seq(0.001, 0.999, length.out = 11), y = -5:5)
ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = qnorm, color = "orange") + 
  stat_function(fun = qnorm, args = list(mean = 1), color = "steelblue") +
  stat_function(fun = qnorm, args = list(sd = 2), color = "darkgray") +
  labs(y = "", title = "Funciones de cuantiles normales")

Distribución Exponencial

Una variable aleatoria \(X\) tienen distribución Exponencial con parámetro \(\beta\), \(X \sim Exp(\beta)\) si,

\[ p(x) = \left\{ \begin{array}{lr} \frac{1}{\beta}e^{-x/\beta} & x >0\\ 0 & e.o.c. \\ \end{array} \right. \] \(E(X)=\beta, Var(X)=\beta^2\)

donde \(\beta > 0\). La distribución exponencial se utiliza para modelar tiempos de espera hasta un evento, por ejemplo modelar el tiempo de vida de un componente electrónico o el tiempo de espera entre llamadas telefónicas.

  • Realiza una gráfica de las funciones de densidad y de distribución acumulada correspondientes a distribuciones exponenciales con parámetros: \(\beta=1, \beta = 2\) y \(\beta=2/3\).


Distribución Gamma

Comencemos definiendo la función Gamma: para \(\alpha>0\), \(\Gamma(\alpha)=\int_0^{\infty}y^{\alpha-1}e^{-y}dy\), esta función es una extensión de la función factorial, tenemos que si \(n\) es un entero positivo, \(\Gamma(n)=(n-1)!\).

Ahora, \(X\) tienen una distribución Gamma con parámetros \(\alpha\), \(\beta\), denotado como \(X \sim Gamma(\alpha, \beta)\) si \[ p(x) = \left\{ \begin{array}{lr} \frac{1}{\beta^\alpha \Gamma(\alpha)}x^{\alpha-1}e^{-x/\beta} & x >0\\ 0 & e.o.c. \\ \end{array} \right. \] \(E(X)=\alpha \beta, Var(X)=\alpha \beta^2\)

Vale la pena notar que una distribución exponencial es una \(Gamma(1, \beta)\). Si \(X_i \sim Gamma(\alpha_i, \beta)\) independientes, entonces \(\sum_{i=1}^n X_i \sim Gamma(\sum_{i=1}^n \alpha_i, \beta)\).

En la práctica la distribución Gamma se ha usado para modelar el tamaño de las reclamaciones de asegurados, en neurociencia se ha usado para describir la distribución de los intervalos entre los que ocurren picos. Finalmente, la distribución Gamma es muy usada en estadística bayesiana como a priori conjugada para el parámetro de precisión de una distribución Normal.

xy <- data.frame(x = 0:12, y = seq(0, 1, length.out = 13))
ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = dgamma, args = list(shape = 1), color = "orange") + 
  stat_function(fun = dgamma, args = list(scale = 0.5, shape = 2), color = "steelblue") +
  stat_function(fun = dgamma, args = list(scale = 3, shape = 4), color = "darkgray") +
  labs(y = "", title = "fdp's gamma")


Distribución Beta

\(X\) tiene una distrinución Beta con parámetros \(\alpha > 0\) y \(\beta >0\), \(X \sim Beta(\alpha, \beta)\) si

\[ p(x) = \left\{ \begin{array}{lr} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1} & 0 < x < 1\\ 0 & e.o.c. \\ \end{array} \right. \] \(E(X)=\alpha/(\alpha+\beta), Var(X)=\alpha \beta /[(\alpha+\beta)^2(\alpha + \beta + 1)]\)


La distribución Beta se ha utilizado para describir variables aleatorias limitadas a intervalos de longitud finita, por ejemplo, distribución del tiempo en sistemas de control o administración de proyectos, proporción de minerales en rocas, etc.

xy <- data.frame(x = 0:1, y = 0:1)
ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = dbeta, args = list(shape1 = 1, shape2 = 3), color = "orange") + 
  stat_function(fun = dbeta, args = list(shape1 = 5, shape2 = 1), color = "steelblue") +
  stat_function(fun = dbeta, args = list(shape1 = 0.5, shape2 = 0.5), color = "darkgray") +
  labs(y = "", title = "fdp's Beta")


Teoría básica de simulación

La simulación de un modelo consiste en la construcción de un programa computacional que permite obtener los valores de las varibles de salida para distintos valores de las variables de entrada con el objetivo de obtener conclusiones del sistema que apoyen la toma de decisiones (explica y/o predecir el comportamiento del sistema.

Requerimientos prácticos de un proyecto de simulación:

Números pseudoaleatorios

El objetivo es generar sucesiones \(\{u_i\}_{i=1}^N\) de números independientes que se puedan considerar como observaciones de una distribución uniforme en el intervalo \((0, 1)\).

  1. Verdaderos números aleatorios Los números completamente aleatorios (no determinísticos) son fáciles de imaginar conceptualmente, por ejemplo podemos imaginar lanzar una moneda, lanzar un dado, una lotería.

En general los números aleatorios se basan en alguna fuente de aleatoreidad física que puede ser teóricamente impredecible (cuántica) o prácticamente impredecible(caótica). Un ejemplo es random.org que genera aleatoreidad a través de ruido atmosférico.

La desventaja de éstos métodos es que son costosos, tardados y no son reproducibles.

  1. Números pseudoaleatorios Los números pseudoaleatorios se generan de manera secuencial con un algoritmo determinístico, formalmente se define por una función de inicialización, un estado, una función de transición y una función que genera salidas:

  • Función de inicialización. Recibe un número (la semilla) y pone al generador en su estado inicial.

  • Función de transición. Transforma el estado del generador.

  • Función de salidas. Transforma el estado para producir un número fijo de bits (0 ó 1).

Una sucesión de bits pseudoaleatorios se obtiene definiendo la semilla y llamando repetidamente la función de transición y la función de salidas. Esto implica, entre otras cosas, que la sucesión de números generados esta completamente determinada por la semilla.

Ejemplo. Por muchos años (antes de 1995) el generador de la función rand en Matlab fue el generador congruencial:

\[X_{n+1} = (7^5)X_n mod(2^{31}-1)\]

library(plyr) # raply
library(reshape2) # melt
library(dplyr)
library(ggplot2) # ggplot

## construyamos sucesiones de longitud n usando el algoritmo de rand
sucesion <- function(n = 1500, semilla = runif(1, 0, 2 ^ 31 - 1)){
  x <- rep(NA, n)
  u <- rep(NA, n)
  x[1] <- semilla # semilla
  u[1] <- x[1] / (2 ^ 31 - 1) # transformamos al (0, 1)
  for(i in 2:n){
    x[i] <- (7 ^ 5 * x[i - 1]) %% (2 ^ 31 - 1)
    u[i] <- x[i] / (2 ^ 31 - 1)
  }
  u
}
sucesiones <- t(raply(12, sucesion())) # generamos 12 sucesiones con 1500 entradas cada una
head(sucesiones)
##         [,1]      [,2]      [,3]       [,4]       [,5]      [,6]
## 1 0.09426519 0.9829031 0.5703200 0.69339678 0.92521191 0.5936606
## 2 0.31504363 0.6521957 0.3684073 0.91965308 0.03651533 0.6542054
## 3 0.93827899 0.4534935 0.8209898 0.60936146 0.71322507 0.2299235
## 4 0.65497656 0.8651914 0.3755790 0.53807845 0.17381546 0.3250295
## 5 0.19105466 0.2722057 0.3558644 0.48450641 0.31641338 0.7709532
## 6 0.05567415 0.9618297 0.0134776 0.09924817 0.95966455 0.4102017
##           [,7]       [,8]       [,9]     [,10]     [,11]      [,12]
## 1 0.4076398057 0.90468862 0.59294541 0.8210418 0.1248874 0.49570726
## 2 0.2022142853 0.10160126 0.63358702 0.2494361 0.9821182 0.35185586
## 3 0.6154924527 0.61245535 0.69698617 0.2732725 0.4613434 0.64140143
## 4 0.5816529433 0.53704316 0.24662618 0.8911644 0.7985221 0.03387472
## 5 0.8410186747 0.08439161 0.04628959 0.8001057 0.7613733 0.33241925
## 6 0.0008663737 0.36976967 0.98921069 0.3767827 0.4014595 0.97032177

Una propiedad deseable es que la sucesión de \(u_i\) parezca una sucesión de observaciones independientes de una \(Uniforme(0,1)\).

sucesiones_m <- melt(sucesiones) # cambiamos los datos a forma "larga"
head(sucesiones_m)
##   Var1 Var2      value
## 1    1    1 0.09426519
## 2    2    1 0.31504363
## 3    3    1 0.93827899
## 4    4    1 0.65497656
## 5    5    1 0.19105466
## 6    6    1 0.05567415
# veamos una gráfica del índice de simulación contra el valor obtenido
ggplot(sucesiones_m, aes(x = Var1, y = value)) + 
  geom_point(alpha = 0.5, size = 1.5) +     # alpha controla la transparencia
  facet_wrap(~ Var2)

# veamos una gráfica del índice de simulación contra el valor obtenido
ggplot(sucesiones_m, aes(x = value)) + 
  geom_histogram(aes(y = ..density..), binwidth = 0.1) +
  facet_wrap(~ Var2)

Una buena secuencia de números pseudoaleatorios no muestra ningún patrón o regularidad aparente desde un punto de vista estadístico. Otra caraceterística importante es que dada una semilla inicial, el número de variables que se pueden generar antes de repetir el ciclo es grande. Construir un buen algoritmo de números pseudo aleatorios es complicado (pseudoaleatorio vs aleatorio).

Ejemplo. Un generador de números aleatorios ampliamente utilizado en los 60´s y 70´s fue RANDU, que se define como: \[X_{n + 1}= (2 ^ {16} + 3)X_n mod(2^{31})\]

A primera vista las sucesiones se asemejan a una uniforme; sin embargo, cuando se grafican ternas emergen patrones no deseados.

# correr este código como shiny app para que funcione la siguiente gráfica
# interactiva
library(tourr)
## 
## Attaching package: 'tourr'
## 
## The following object is masked from 'package:plyr':
## 
##     ozone
library(ggvis)
## 
## Attaching package: 'ggvis'
## 
## The following object is masked from 'package:ggplot2':
## 
##     resolution
library(shiny)


n <- 1500 # longitud de la sucesión
x <- rep(NA, n)
u <- rep(NA, n)
 
x[1] <- 4798374 # semilla
u[1] <- x[1] / (2 ^ 31 - 1) # transformamos al (0, 1)
for(i in 2:n){
  x[i] <- ((2 ^ 16 + 3) * x[i - 1]) %% (2 ^ 31)
  u[i] <- x[i] / (2 ^ 31)
}

aps <- 2
fps <- 30

mat <- scale(matrix(x, ncol = 3, byrow= TRUE))
tour <- new_tour(mat, grand_tour(), NULL)
start <- tour(0)

proj_data <- reactive({
  invalidateLater(1000 / fps, NULL);
  step <- tour(aps / fps)
  data.frame(center(mat %*% step$proj))
})

proj_data %>% ggvis(~X1, ~X2) %>%
  layer_points() %>%
  scale_numeric("x", domain = c(-1, 1)) %>%
  scale_numeric("y", domain = c(-1, 1)) %>%
  set_options(duration = 20)
## Warning: Can't output dynamic/interactive ggvis plots in a knitr document.
## Generating a static (non-dynamic, non-interactive) version of the plot.

El generador actual de números aleatorios de R se conoce como “Mersenne-Twister” pero hay varios generadores disponbiles (ver ?Random).

En el ejemplo de arriba vimos que el generador randu produce ternas que no se asemejan a una distribución uniforme en \(\mathbb{R}^3\), hacer gráficas no siempre es la mejor herramienta para evaluar un generador, sino que existen pruebas estadísticas para determinar la calidad de un generador de números pseudoaleatorios, por ejemplo la prueba de bondad de ajuste \(\chi^2\), prueba de permutación de índices y prueba de espera (Gap-test).

Generación de variables aleatorias discretas

Método de Inversión

Supongamos que deseamos generar el valor de una variable aleatoria discreta \(X\) con función de probabilidad: \[P(X=x_j) = p_j\] con \(j=1,2,..\).

Para lograr esto generamos un número aleatorio \(U\), esto es \(U\sim Uniforme(0,1)\) y definimos

\[ X = \left\{ \begin{array}{lr} x_0 & U < p_0\\ x_1 & p_0 \leq U < p_0 + p_1\\ \vdots &\\ x_j & \sum_{i=0}^{j-1}p_i \leq U < \sum_{i=0}^j p_i \\ \vdots & \\ \end{array} \right. \]

Como para \(0<a<b<1\) tenemos que \(P(a\leq U < b)=b-a\), tenemos que \[P(X=x_j)=P\bigg\{\sum_{i=0}^{j-1}p_i \leq U < \sum_{i=0}^{j}p_i \bigg \}=p_j\] y por tanto \(X\) tiene la distribución deseada.

Método de inversión

  1. Genera un número aleatorio \(U\), tal que \(U \in (0,1)\).
    Si \(U<p_0\) define \(X=x_0\) y para.
    Si \(U< p_0+p_1\) define \(X = x_1\) y para.
    Si \(U < p_0 + p_1 + p_2\) define \(X=x_2\) y para.
    \(\vdots\)

  2. Si las \(x_i\), están ordenadas de tal manera que \(x_0<x_1<x_2<\cdots\) y si denotamos por \(P\) la función de distribución acumulada de \(X\), entonces \(P(x_k)=\sum_{i=0}^kp_i\) y por tanto, \(X\) será igual a \(x_j\) si \[P(x_{j-1}) \leq U \leq P(x_j)\]. En otras palabras, tras generar un número aleatorio \(U\) determinamos el valor de \(X\) encontrando el intervalo \([P(x_{j-1}),P(x_j))\) en el que cae \(U\), esto es equivalente a encontrar la inversa de \(P(U)\).

El tiempo que uno tarda en generar una variable aleatoria discreta usando el método de arriba es proporcional al número de intervalos que uno debe buscar, es por esto que en ocasiones vale la pena considerar los posibles valores \(x_j\) en orden decreciente de \(p_j\).

  • Utiliza la función runif de R y el método de inversión para generar 1000 simulaciones de una variable aleatoria \(X\) tal que \(p_1=0.20, p_2= 0.15, p_3=0.25, p_4=0.40\) donde \(p_j=P(X=j)\).
u <- runif(1)

u <- runif(1000)
x <- 1 * (u < 0.20) + 2 * (0.20 <= u & u < 0.20 + 0.15) + 
     3 * (0.20 + 0.15 <= u & u < 0.20 + 0.15 + .25) + 
     4 * (0.20 + 0.15 + .25 <= u)

hist(x)

Ejemplo Uniforme discreta: Supongamos que deseamos simular de una variable aleatoria uniforme discreta que toma valores \(1,...,k\), usando los resultados anteriores tenemos que:

\(X=j\) si \(\frac{j-1}{n} \leq U < \frac{j}{n}\)

Entonces \(X=[kU] + 1\), donde \([x]\) representa la parte entera de x.

# uniforme discreta: donde n es el número de simulaciones y k el número de elementos
runifD <- function(n = 1, k) floor(k*runif(n)) + 1
# veamos un histograma de 1000 simulaciones de una distribución Uniforme
# discreta con parámetro k = 20
hist(runifD(n = 1000, k = 20))

# También podmeos usar la función sample de R
# hist(sample(1:20, size = 1000, replace= TRUE))

Ejemplo Poisson: la clave para usar el método de la transformación inversa para generar esta variable aleatoria es notar que: \[p_{i+1}=\frac{\lambda}{i+1}p_i\]

donde \(p_i=P(X=i) = e^-{\lambda} \lambda^i/i!\), con \(i=0,1,...\). Ahora, la cantidad \(i\) se refiere al valor que estamos considerando, \(p=p_i\) es la probabilidad de \(X = i\) y \(P=P(i)\) es la probabilidad de \(X\leq i\). Entonces, para generar una observación sequimos los siguientes pasos:

  1. Generar un número aleatorio \(U\).
  2. Inicializar: \(i=0\), \(p=e^{-\lambda}\), \(F=p\).
  3. Si \(U<F\), definir \(X=i\) y parar.
  4. \(p=\lambda p/(i+1)\), \(F=F+p\), \(i=i+1\).
  5. Volver a 3.
# Poisson usando Inversión
rpoisI <- function(lambda = 1){
  U <- runif(1)
  i <- 0
  p <- exp(-lambda)
  P <- p
  while(U >= P){
    p <- lambda * p / (i + 1)
    P <- P + p
    i <- i + 1
  }
  i
}
hist(rdply(1000, rpoisI)$V1)

El algoritmo que propusimos verifica de manera sucesiva si el valor es 0, 1, etc. por lo que el número de comparaciones necesarias será uno más que el valor de la variable. Ahora, el valor esperado de una variable aleatoria Poisson es \(\lambda\) por lo que en promedio se harían \(1+\lambda\) busquedas. Cuando \(\lambda\) es grande se puede mejorar el algotitmo buscando primero en valores cercanos a \(\lambda\).

  • Escribe una función en R que genere simulaciones de una variable aleatoria Poisson de la siguiente manera: define \(I=[\lambda]\), y usa que \(p_{i+1}=\lambda p_i /(i+1)\) para determinar \(F\) de manera recursiva. Genera un número aleatorio \(U\), determina si \(X \leq I\) comparando si \(U \leq F(I)\). Si \(X \leq I\) busca hacia abajo comenzando en \(I\), de lo contrario busca hacia arriba comenzando por \(I+1\). Compara el tiempo que tardan los dos algoritmos en 5000 simulaciones de una variable aleatoria Poisson con parámetro \(\lambda=10, 200, 500\).

Aceptación y rechazo

Supongamos que tenemos un método eficiente para generar simulaciones de una variable aleatoria con función de probabilidad masa \(\{q_j, j\geq 0\}\), podemos usarla como la base para simular de una distribución que tiene función de probabilidad masa \(\{p_j, j \geq 0\}\), para hacer esto comenzamos simulando una variable aleatoria \(Y\) con función \(\{q_j\}\) y después aceptamos o rechazamos el valor simulado con una probabilidad proporcional a \(p_Y/q_Y\). En particular, sea \(c\) una constante tal que \[\frac{p_j}{q_j}\leq c\] para toda \(j\) con \(p_j > 0\). Entonces el método de aceptación y rechazo para simular una variable aleatoria \(X\) con función masa de probabilidad \(p_j=P(x=j)\) es como sigue:

Método de aceptación y rechazo

  1. Simula el valor de \(Y\), con función de probabilidad masa \(q_j\).
  2. Genera un número aleatorio \(U\).
  3. Si \(U < p_Y/(cq_y)\) definimos \(X=Y\) y paramos, en otro caso regresamos a 1.

Ejemplo: Supongamos que queremos simular el valor de una variable aleatoria \(X\) que toma uno de los valores \(1, 2,..., 10\) con probabilidades \(0.15,0.12, 0.09, 0.08, 0.12, 0.10,0.06,0.09,0.10,0.10\). Usemos el método de aceptación y rechazo con \(q\) la densidad uniforme en \(1,...,10\). Podemos escoger \(c\) como \[c=max\frac{p_j}{q_j}=1.2\]

simAR <- function(){
  p <- c(0.12, 0.12, 0.09, 0.05, 0.15, 0.06, 0.09, 0.09, 0.10, 0.10)
  # Generamos un valor de la uniforme discreta
  u_d <- floor(10 * runif(1)) + 1
  # Generamos un segundo número aleatorio
  u <- runif(1)
  # Calculamos p_j
  p_j <- p[u_d]
  # Si u <=p/0.12 x = y
  while(u >= p_j / (1/10 * 1.5)){  
    u_d <- floor(10 * runif(1)) + 1
    u <- runif(1)
    p_j <- p[u_d]
  }
  x <- u_d 
}

sims <- rdply(10000, simAR())

round( prop.table(table(sims$V1)), 2)
## 
##    1    2    3    4    5    6    7    8    9   10 
## 0.12 0.12 0.10 0.05 0.16 0.06 0.09 0.09 0.11 0.11

En promedio este algoritmo requiere \(1.2\) iteraciones para obtener un valor generado para \(X\).

Referencias