D.3 Capítulo 5 Simulación de variables discretas

D.3.1 Ejercicio 5.1

Enunciado 5.1 (Simulación de una distribución mixta mediante el método de inversión generalizado):

Consideramos la variable aleatoria con función de distribución dada por: \[F(x)=\left\{ \begin{array} [c]{cl}0 & \mbox{si $x<0$}\\ \frac{x}{2}+\frac{1}{10} & \mbox{si $x\in[0,\frac{1}{5})$}\\ x+\frac{1}{10} & \mbox{si $x\in[\frac{1}{5},\frac{9}{10}]$}\\ 1 & \mbox{en otro caso} \end{array} \right.\]

Esta función está implementada en el siguiente código:

fdistr <- function(x) {
ifelse(x < 0, 0,
    ifelse(x < 1/5, x/2 + 1/10,
        ifelse(x <= 9/10, x + 1/10, 1) ) )
}
# Empleando ifelse la función es vectorial (y podemos emplear curve...)
curve(fdistr, from = -0.1, to = 1.1, type = 's', 
      main = 'Función de distribución')
# Discontinuidades en 0 y 1/5
abline(h = c(1/10, 2/10, 3/10), lty = 2) 

Nota: Esta variable toma los valores 0 y 1/5 con probabilidad 1/10.

1. Diseñar un algoritmo basado en el método de inversión generalizado para generar observaciones de esta variable.

2. Implementar el algoritmo en una función que permita generar \(nsim\) valores de esta variable.

---

1. El algoritmo general es siempre el mismo. Empleando la función cuantil: \[Q\left( u\right) = \inf \left\{ x\in \mathbb{R}:F\left( x\right) \geq u\right\},\] el algoritmo sería:

   1. Generar \(U\sim \mathcal{U}\left( 0,1\right)\)

   2. Devolver \(X=Q\left( U\right)\)

   En este caso concreto:

   1. Generar \(U\sim \mathcal{U}\left( 0,1\right)\)

   2. Si \(U < \frac{1}{10}\) devolver \(X = 0\)

   3. Si \(U < \frac{2}{10}\) devolver \(X = 2(U - \frac{1}{10})\)

   4. Si \(U < \frac{3}{10}\) devolver \(X = \frac{2}{10}\)

   5. En caso contrario devolver \(X = U - \frac{1}{10}\)

2. El algoritmo de simulación se puede implementar a partir de la función cuantil (vectorial):

   # Función cuantil:
   fquant <- function(u) 
     ifelse(u < 1/10, 0,
            ifelse(u < 2/10, 2*(u - 1/10),
                   ifelse(u < 3/10, 1/5, u - 1/10) ) )
   # Función para generar nsim valores:
   rx <- function(nsim) fquant(runif(nsim))

   Ejemplo:

   set.seed(1)
   nsim <- 10^4
   system.time(simx <- rx(nsim))

   ##    user  system elapsed 
   ##       0       0       0

   hist(simx, breaks = "FD", freq = FALSE)

   

   En este caso como no es una variable absolutamente continua mejor emplear la función de distribución para compararla con la teórica:

   curve(ecdf(simx)(x), from= -0.1, to = 1.1, type = "s")
   curve(fdistr(x), type = "s", lty = 2, add = TRUE)