Método de transformación inversa

Alex J. Zambrano
alexzambrano@usantotomas.edu.co

Cálculo numérico y Simulación

Generación de v.a. continuas

Sea \(U\) una v.a. uniforme \((0,1)\). Para cualquier función de distribución continua F, invertible, la variable aleatoria \(X\) definida como \[ X=F^{-1}(U) \] tiene distribución F.

Para demostrar lo anterior, decimos que \(F_X\) es la función de distribución de \(X=F^{-1}(U)\), por lo cual \[ \begin{align*} F_X(x)&=P(X\leq x)\\ &=P(F^{-1}(U)\leq x)\\ &=P(F(F^{-1}(U))\leq F(x))\quad\text{pues $F$ es monotona, creciente se cumple que $a<b$ entonces $F(a)<F(b)$}\\ &=P(U\leq F(x))\quad\text{pues $F(F^{-1}(U))=U$}\\ &=F(x)\quad\text{pues $U$ es uniforme en (0,1)}\\ \end{align*} \] Un algoritmo para el método de la transformación inversa es el sguiente

  1. Generar \(U\sim U(0,1)\)
  2. Hacer \(X=F^{-1}(U)\)
  3. Mostrar X.

Ejemplo

Distribución exponencial \[ f(x)=\begin{cases} \lambda e^{-\lambda x}&\text{ si } x\geq 0\\ 0 & \text{e. o. c.} \end{cases} \] La distribución acumulada viene dada por \[ F(x)=\begin{cases} 1-e^{-\lambda x}&\text{ si } x\geq 0\\ 0 & \text{e. o. c.} \end{cases} \] Haciendo \(X=F^{-1}(U)=-\frac{1}{\lambda}\log(1-U)=-\frac{1}{\lambda}\log(U)\), ya que \(1-U\sim U(0,1)\). Por lo cual el algoritmo nos quedaría de la siguiente manera:

  1. Genera \(U\sim U(0,1)\).
  2. Hacer \(X=-\frac{1}{\lambda}\log(U)\).
  3. Mostrar X.

utilizando R, nos quedaría de la siguiente manera

exponencial <- function(n,a){
  U <- runif(n)
  X <- -1/a*log(U)
}

Ahora probemos para \(\lambda=1/5\), generemos 1000 valores.

X <-exponencial(n = 1000,a = 1/5)

Realizando un histograma y comparando los resultados tenemos lo siguiente

hist(X,freq = F)

curve(dexp(x,rate=1/5),add=T,col=2)

Una forma teórica para probar si los valores provienen de una distribución exponencial, se podría utilizar la prueba de kolmogorov-Smirnov

ks.test(X,"pexp",rate=1/5)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  X
## D = 0.017686, p-value = 0.9132
## alternative hypothesis: two-sided

Ejemplo

Distribución Weibull \[ f(x)=\begin{cases} \alpha\beta(\beta x)^{\alpha-1}e^{-(\beta x)^\alpha}&\text{ si } x\geq 0\\ 0 & \text{e. o. c.} \end{cases} \] La distribución acumulada viene dada por \[ F(x)=\begin{cases} 1-e^{-(\beta x)^\alpha}&\text{ si } x\geq 0\\ 0 & \text{e. o. c.} \end{cases} \] Haciendo \(X=F^{-1}(U)=(-\frac{1}{\beta^\alpha}\log(U))^{1/\alpha}\). Por lo cual el algoritmo nos quedaría de la siguiente manera:

  1. Genera \(U\sim U(0,1)\).
  2. Hacer \(X=(-\frac{1}{\beta^\alpha}\log(U))^{1/\alpha}\).
  3. Mostrar X.

utilizando R, nos quedaría de la siguiente manera

weibull <- function(n,a,b){
  U <- runif(n)
  X <- (-(1/b^a)*log(U))^(1/a)
}

Ahora probemos para \(\alpha=2\), \(\beta=3\) generemos 1000 valores.

X <-weibull(n = 1000,a = 2,b = 3)

Realizando un histograma y comparando los resultados tenemos lo siguiente

hist(X,freq = F)

curve(dweibull(x,shape = 2,scale = 1/3),add=T,col=2)

Realizando la prueba de bondad de ajuste se tiene que Kolmogorov-Smirnov

ks.test(X,"pweibull",shape=2,scale=1/3)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  X
## D = 0.024788, p-value = 0.5706
## alternative hypothesis: two-sided

Ejercicios

A continuación se describen algunas distribuciones las cuales tienen expersiones analiticas simples.

  1. Distribución Burr \[ f(x)=ckx^{c-1}(1+x^c)^{-k-1} \]
  2. Distribución Cauchy \[ f(x)=\frac{\beta}{\pi[\beta^2+(x-\alpha)^2]} \]
  3. Distribución Laplace \[ f(x)=\frac{1}{2}e^{|-x|} \]
  4. Distribución Logística \[ f(x)=\frac{e^{-x}}{[1+e^{-x}]^2} \]
  5. Distribución Pareto \[ f(x)=\frac{a}{x^a+1} \]
  6. Distribución Fisher \[ f(\theta)=\frac{ke^{k\cos(\theta)}\sin(\theta)}{2\sinh(k)}\quad 0<\theta<\pi,k>0 \]

Generación de v.a. discretas

Supongamos que deseamos generar el valor de una v.a. discreta teniendo función de masa de probabildad \[ P\left\{X=x_i\right\}=p_i\quad i=0,1,\ldots,\quad \sum_{i}p_i=1 \] Teniendo esto, generamos un número aleatorio \(U\sim U(0,1)\) y el conjunto \[ X=\begin{cases} x_0 & \text{ si } U<p_0\\ x_1 & \text{ si } p_0\leq U<p_0+p_1\\ x_2 & \text{ si } p_0+p_1 \leq U<p_0 + p_1 + p_2\\ \vdots & \\ x_i & \text{ si } \sum_{j=1}^{i-1}p_j\leq U<\sum_{j=1}^{i}p_j\\ \vdots & \\ \end{cases} \] Puesto que \(F_i=\sum_{j=1}^{i}p_j\) para \(i=0,1,2,\dots\), lo anterior se puede reescribir de la siguiente manera \[ X=\begin{cases} x_0 & \text{ si } U< F_0\\ x_1 & \text{ si } F_0\leq U< F_1\\ x_2 & \text{ si } F_1 \leq U< F_2\\ \vdots & \\ x_i & \text{ si } F_{i-1}\leq U< F_i\\ \vdots & \\ \end{cases} \] Escribiendo lo anterior en forma algoritmica tenemos que

  1. Generar \(U\sim U(0,1)\)

  2. Si \(U < F_0\) entonces \(X=x_0\)

  3. Si \(U < F_1\) entonces \(X=x_1\)

  4. \(\vdots\)

Ejemplo

Generar una v.a. discreta finita, dada por

\(i\) 1 2 3 4
\(p_i\) 0.20 0.15 0.25 0.4
\(F_i\) 0.20 0.35 0.60 1

Los pasos para el algoritmo serían los siguientes:

  1. Generar \(U\sim U(0,1)\).
  2. Si \(U < 0.20\) entonces \(X=1\)
  3. Si \(U < 0.35\) entonces \(X=2\)
  4. Si \(U < 0.60\) entonces \(X=3\)
  5. En otro caso \(X=4\).

Utilizando R se tiene que

U <- runif(1000)
X<-numeric(1000)
for(i in 1:1000){
  if(U[i]< 0.20){
    X[i] <- 1
  } else if (U[i]< 0.35){
    X[i] <- 2
  } else if (U[i]< 0.60){
    X[i] <- 3
  } else{
    X[i] <- 4
  }
}
tabla <- table(X)/1000
tabla
## X
##     1     2     3     4 
## 0.219 0.137 0.247 0.397
barplot(tabla)

Es más eficiente el tiempo de busqueda si se ordenan los \(p_i\) en orden decreciente. Por lo cual realizando en R tenemos que

U <- runif(1000)
X<-numeric(1000)
for(i in 1:1000){
  if(U[i]< 0.40){
    X[i] <- 4
  } else if (U[i]< 0.65){
    X[i] <- 3
  } else if (U[i]< 0.85){
    X[i] <- 1
  } else{
    X[i] <- 2
  }
}
table(X)/1000
## X
##     1     2     3     4 
## 0.195 0.125 0.258 0.422
tabla <- table(X)/1000
tabla
## X
##     1     2     3     4 
## 0.195 0.125 0.258 0.422
barplot(tabla)

Método general de la transformación inversa para el caso discreto

Rubinstein (1981), propone un algoritmo por el método de transformación inversa para utilizar con funciones de probabilidad ya sea discretas finitas o infinitas. Este algortimo se basa en compara el valor U con cada valor de \(F_i\) y se llama algoritmo IT-2

  1. \(p=P_0\).
  2. \(F=p\).
  3. \(i=0\).
  4. Generar \(U\sim U(0,1)\).
  5. Si \(U\leq F\), hacer \(X=i\).
  6. \(p=A_{k+1}p\).
  7. \(F=F+p\).
  8. \(i=i+1\).
  9. Retorne al paso 5.

Donde \(P_0\) y \(A_{k+1}=P_{k+1}/P_k\) se distribuyen dependientemente.

Ejemplo

Generar valores de una v.a. poisson, cuya función de densidad viene dada por: \[ p_i=P\{X=i\}=e^{-\lambda}\frac{\lambda^i}{i!}\quad i=0,1,\ldots \] Recordemos que \(p_0=e^{-\lambda}\) y además \(A_{i+1}=p_{i+1}/p_{i}=\lambda/(i+1)\), por lo cual el algoritmo queda de la siguiente manera

Poisson.gen <- function(n,lambda){
  X <- numeric(n)
  for(k in 1:n){
    p <- exp(-lambda)
    F <- p
    i <- 0
    U <- runif(1)
    j <- 0
    while(j==0){
      if(U < F){j <- 1}
      X[k] <- i
      p <- lambda/(i+1)*p
      F <- F + p
      i <- i+1
    }
  }
  X
}

Ahora generamos 10000 valores para la distribución poisson con \(\lambda=5\)

X <- Poisson.gen(n = 10000,lambda = 5)
tabla <- table(X)/10000
tabla
## X
##      0      1      2      3      4      5      6      7      8      9 
## 0.0079 0.0362 0.0839 0.1436 0.1759 0.1727 0.1487 0.0979 0.0688 0.0350 
##     10     11     12     13     14 
## 0.0175 0.0077 0.0029 0.0009 0.0004

Comparando con los valores teóricos tenemos que

sort(unique(X))
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14
round(dpois(x=sort(unique(X)),lambda = 5),4)
##  [1] 0.0067 0.0337 0.0842 0.1404 0.1755 0.1755 0.1462 0.1044 0.0653 0.0363
## [11] 0.0181 0.0082 0.0034 0.0013 0.0005

Comparanda los diagrama de barras para ambos resultados tenemos que

par(mfrow=c(1,2))
barplot(tabla,main="Valores simulados")
barplot(height = dpois(x=sort(unique(X)),lambda = 5),names.arg = sort(unique(X)),col="red",main="Valores teóricos")

par(mfrow=c(1,1))

Otra forma comprobar que los valores simulados siguen la distribución propuesta es calculando los parámetros de interes tales como la media y la varianza. Para esté caso la media y la varianza deben ser iguales a 5.

mean(X)
## [1] 4.9566
var(X)
## [1] 4.947011

Esté mismo método funciona para distribuciones que provienen de una tabla de distribución discreta.

Ejercicios

Mediante el método general de la transformación inversa para el caso discreto, estudie los siguientes casos:

  1. Distribución binomial
  2. Distribución geométrica
  3. Distribución binomial-negativa.
  4. Escriba un programa de computadora tal que, dada una función de masa de probabilidad \(\{p_j, j = 1,\ldots,n\}\) como entrada, proporcione como salida el valor de una variable aleatoria con esta función de masa.
  5. Suponga que en un banco hay un cierto número de personas esperando a que lo habran para ser atendidas. Al momento de abrir las puertas del banco una persona queda en el puesto 50 en la única fila que hay habilitada para la atención. Si minutos más tarde dos cajeros del banco comienzan a atender con tiempos de atención exponencial con medias 2 y 3 minutos respectivamente. Cuál es el tiempo esperado que esa persona durará en la fila hasta ser atendido teniendo en cuenta solo el tiempo a partir del momento en que comienzan a atender.

Bibliografía

Ross, Sheldon M. 1999. Simulación. Segunda. Prentice Hall.

Rubinstein, Reuven Y. 1981. Simulation and Monte Carlo Method. Jhon Wiley & Sons.