Métodos universales para la simulación de variables discretas

Ejercicio 1.

Cuando se produce una sobrecarga eléctrica en una cadena productiva el número de componentes averiados está descrito mediante una variable aleatoria, X, que puede tomar los valores 1, 2, 3, 4 y 5, con probabilidades,

\[\begin{align*} \begin{array}{rrr} p_{1} = 0.11, & p_{2}=0.22, & p_{3}=0.45, & p_{4}=0.13, & p_{5}=0.09. \end{array} \end{align*}\]

Para simular dicha variable se usa el algoritmo:

1. Generar U, con distribucion U(0,1)

2.1 Si \(U \leq 0.11\) entonces hacer \(X=1\), sino

2.2 Si \(U \leq 0.33\) entonces hacer \(X=2\), sino

2.3 Si \(U \leq 0.78\) entonces hacer \(X=3\), sino

2.4 Si \(U \leq 0.91\) entonces hacer \(X=4\), sino hacer \(X=5\)

3. Devolver \(X\).

¿Cuál es el número medio de compraciones para simular un valor de X? Construir un algoritmo basado en un árbol binario óptimo. ¿Cuál es ahora el número medio de comparaciones? Dar un algoritmo basado en el método de la tabla guía. Calcular, también en este caso, el número medio de comparaciones.

#=======================================================================================
# Algoritmo_1
#=======================================================================================
ncomp <<- 0
sim1N <- function(n=1000){
  X <- numeric(n)
  for(isim in 1:n)
  {
    U <- runif(1,0,1)
    if(U<=0.11){X[isim]<-1
                 ncomp<<-ncomp+1}
    else{if(U<=0.33){X[isim]<-2
                      ncomp<<-ncomp+2}
      else{if(U<=0.78){X[isim]<-3
                        ncomp<<-ncomp+3}
        else{if(U<=0.91){X[isim]<-4
                          ncomp<<-ncomp+4}
          else{X[isim]<-5
                ncomp<<-ncomp+4}}}}
    
  }
  return(X)
}

## Nº medio de comparaciones =  2.7749

Árbol Binario

Tenemos el siguiente árbol binario:

Donde las probabilidades asociadas a los nodos están dados por:

Nodos {1,2,3,4,5 } {1,4,5 } {2,3 } {1,5}
p 1 0.33 0.67 0.2

Asi pues, proponemos el sigueinte algoritmo.

#=======================================================================================
# Algoritmo_2
#Arbol Binario Óptimo
#Obteniendo previamente el árbol binario óptimo
#=======================================================================================

ncomp <<- 0
sim2N <- function(n=1000){
  X <- numeric(n)
  for (isim in 1:n) {
    U <- runif(1,0,1)
    if (U>0.33){if(U>0.78){X[isim]<-2
                            ncomp <<- ncomp+2}
                else{X[isim]<-3
                      ncomp <<- ncomp+2}}
    else{if(U>0.2){X[isim]<-4
                    ncomp <<- ncomp+2}
          else{if(U>0.11){X[isim]<-5
                          ncomp <<- ncomp+3}
                else{X[isim]<-1
                      ncomp <<- ncomp+3}}}
    
  }
  return(X)
  return(ncomp)
}

## Nº medio de comparaciones =  2.1998
#=======================================================================================
# Algoritmo_3
# Tabla Guía. 
#=======================================================================================
#---------------------------------------------------------------------------------------
x <- c(1:5)
p <- c(0.11,0.22,0.45,0.13,0.09)
#---------------------------------------------------------------------------------------
ncomp <<- 0
rfmp.tabla <- function(x, prob = 1/length(x), m, nsim = 1000) {
  # Simulación v.a. discreta a partir de función de masa de probabilidad
  # por tabla guia de tamaño m
  # Inicializar tabla y FD
  Fx <- cumsum(prob)
  g <- rep(1,m)
  i <- 1
  for(j in 2:m) {
    while (Fx[i] < (j-1)/m) i <- i+1
    g[j] <- i
  }
  # Generar valores
  X <- numeric(nsim)
  U <- runif(nsim)
  for(j in 1:nsim) {
    i <- g[floor(U[j]*m)+1]
    while (Fx[i] < U[j]) i <- i + 1
    X[j] <- x[i]
    ncomp <<- ncomp+1
  }
  return(X)
}

## Nº medio de comparaciones =  1

Notamos que a medida que vamos mejorando los algoritmos el numero medio de comparaciones disminuye, llegando a ser 1 al momento de usar la tabla guía, lo cual significaria una gran eficiencia para simular esta variable.

Ejercicio 2.

El número de procesos de usuarios activos en una estación de trabajo es una variable aleatoria, \(X\), con distribución dada por:

\[\begin{align*} \begin{array}{rrr} & P(X=0) = 0.33, & P(X=1)=0.22, & P(X=2)=0.15.\\ & P(X=3)=0.11 &, P(x=4)=0.09, & P(X=5)=0.06.\\ & P(X=6)=00.4.& &\\ \end{array} \end{align*}\]

Dar un algoritmo para simular esta variable mediante búsqueda en forma de Árbol de Huffman. Calcular el número medio de comparaciones para simular cada valor de la variable \(X\). ¿Cuál es una cota para dicho número medio de comparaciones en caso de usar el método de la tabla guía? Para este Último método, calcular el valor exacto del número medio de comparaciones.

Árbol Binario

Tenemos el siguiente árbol binario:

Donde las probabilidades asociadas a los nodos están dados por:

Nodos {0,1,2,3,4,5,6 } {2,3,4,5,6 } {0,1 } {3,4,5,6} {4,5,6 } {5,6}
p 1 0.45 0.55 0.30 0.19 0.1

Asi pues, proponemos el siguiente algoritmo.

Proponemos el siguiente algoritmo

sim21N <- function(n=10000){
  X <- numeric(n)
  
  for(i in 1:n){
    U <- runif(1)
    if(U>0.45){if(U>0.78){X[i]<-1
                          ncomp1<<-ncomp1+2}
                else{X[i]<-0
                      ncomp0<<-ncomp0+2}}
    else{if(U>0.3){X[i]<-2
                    ncomp2<<-ncomp2+3}
          else{if(U>0.19){X[i]<-3
                          ncomp3<<-ncomp3+4}
                else{if(U>0.1){X[i]<-4
                                ncomp4<<-ncomp4+5}
                  else{if(U>0.04){X[i]<-5
                                  ncomp5<<-ncomp5+6}
                        else{X[i]<-6
                              ncomp6<<-ncomp6+6}}}}}
    
  }
  return(X)
  
}

Número medio de comparaciones para cada elemento

## Número medio de comparaciones para 0  = 2.011515
##  Número medio de comparaciones para 1  = 1.945455
##  Número medio de comparaciones para 2  = 3.058
##  Número medio de comparaciones para 3  = 4.032727
##  Número medio de comparaciones para 4  = 5.072222
##  Número medio de comparaciones para 5  = 5.73
##  Número medio de comparaciones para 6  = 6.255

Usando el algoritmo de tabla guía se espera que el número de comparaciones para cada elemento sea 1. Asi pues, tenemos:

ncomp <<- 0
rfmp.tabla <- function(x, prob = 1/length(x), m, nsim = 1000) {
  # Simulación v.a. discreta a partir de función de masa de probabilidad
  # por tabla guía de tamaño m
  # Inicializar tabla y FD
  Fx <- cumsum(prob)
  g <- rep(1,m)
  i <- 1
  for(j in 2:m) {
    while (Fx[i] < (j-1)/m) i <- i+1
    g[j] <- i
  }
  # Generar valores
  X <- numeric(nsim)
  U <- runif(nsim)
  for(j in 1:nsim) {
    i <- g[floor(U[j]*m)+1]
    while (Fx[i] < U[j]) i <- i + 1
    X[j] <- x[i]
    ncomp <<- ncomp+1
  }
  return(X)
}

## Nº medio de comparaciones =  1

Ejercicio 3.

El número de veces que se produce una caida de un sistema informatico durante un mes es una variable aleatoria, X, con distribucion dada por \(P(X=0)=0.09, P(X=1)=0.21, P(X=2)=0.39, P(X=3)=0.19, P(X=4)=0.08, P(X=5)=0.04\). Dar un algoritmo para simular esta variable mediante busqueda en forma de Árbol de Huffman. Calcular el número medio de comparaciones para simular cada valor de la variable \(X\). ¿Qué ventajas se obtienen para este método de simulación en comparación con la implementación mediante el método secuencial (tanto en su versión directa como en la del etiquetado mas eficiente)?

Árbol Binario

Tenemos el siguiente árbol binario:

Donde las probabilidades asociadas a los nodos están dados por:

Nodos {0,1,2,3,4,5 } {0,3,4,5 } {1,2 } {0,4,5} {4,5}
p 1 0.40 0.60 0.21 0.12

Asi pues, proponemos el siguiente algoritmo.

#Datos Iniciales
s <- c(0:5)
p <- c(0.09,0.21,0.39,0.19,0.08,0.04)
sim31N <- function(n=1000){
  X <- numeric(n)
  for (i in 1:n) {
    U <- runif(1)
    if(U>0.4){if(U>0.79){X[i]<-1
                          ncomp1 <<- ncomp1+2}
              else{X[i]<-2
                    ncomp2 <<- ncomp2+2}}
    else{if(U>0.21){X[i]<-3
                    ncomp3 <<- ncomp3+3}
          else{if(U>0.12){X[i]<-0
                          ncomp0 <<- ncomp0+4}
                else{if(U>0.08){X[i]<-5
                                ncomp5 <<- ncomp5+5}
                      else{X[i]<-4
                            ncomp4 <<- ncomp4+5}}}
    
  }
  
  
    }
  return(X)
}
system.time(u <- sim31N(nsim))
##    user  system elapsed 
##    0.08    0.00    0.08

## Numero de comparaciones para 0 = 4.08
##  Numero de comparaciones para 1 = 1.948571
##  Numero de comparaciones para 2 = 2.041026
##  Numero de comparaciones para 3 = 2.97
##  Numero de comparaciones para 4 = 5.10625
##  Numero de comparaciones para 5 = 4.475
#=======================================================================
# Método de busqueda secuencial
#=======================================================================
rfmp <- function(x, prob = 1/length(x), nsim = 1000) {
  # Simulación nsim v.a. discreta a partir de fmp
  # por inversión generalizada (transformación cuantil)
  # Inicializar FD
  Fx <- cumsum(prob)
  # Simular
  X <- numeric(nsim)
  U <- runif(nsim)
  for(j in 1:nsim) {
    i <- 1
    while (Fx[i] < U[j]) i <- i + 1
    X[j] <- x[i]
    ncomp <<- ncomp + i
  }
  return(X)
}
system.time(v <- rfmp(s,p,nsim))
##    user  system elapsed 
##    0.02    0.00    0.01

## Nº medio de comparaciones =  3.0689
#=======================================================================
# Método secuencial, ordenando previamente
#=======================================================================

so <- c(2,1,3,0,4,5)
po <- c(0.39,0.21,0.19,0.09,0.08,0.04)
system.time(v <- rfmp(so,po,nsim))
##    user  system elapsed 
##       0       0       0

## Nº medio de comparaciones =  2.3845

El método de árbol binario tiene ventaja al ser mas sencillo de implementar, aun asi, no es tan eficiente como los demas métodos, al menos para este ejercicio, como se ve, el tiempo de computo del método de árbol binario es el mas grande.

Ejercicio 4.

El código de error que puede suministrar una rutina implementada en una lenguaje de alto nivel puede describirse mediante una variable aleatoria, \(C\), que puede tomar los valores 0, 1, 3, 4 y 7, con probabilidades.

\[\begin{align*} \begin{array}{rrr} &p_{0}=0.28, &p_{1}=0.35, & p_{3}=0.16,& p_{4}=0.12, &p_{7}=0.09. \end{array} \end{align*}\]

Dentro de un complejo modelo de simulación se quiere generar dicha variable usando para cada valor generado, un unico numero pseudoaleatorio U(0,1). Encontrar algoritmos, para tal fin, basados en la construcción de un árbol de Huffman y en el método de las tablas guía. Estudiar la eficiencia de ambos algoritmos y compararla.

Árbol Binario

Tenemos el siguiente árbol binario:

Donde las probabilidades asociadas a los nodos están dados por:

Nodos {0,1,3,4,7 } {3,4,7 } {0,1 } {4,7}
p 1 0.37 0.63 0.21

Asi pues, proponemos el siguiente algoritmo.

s <- c(0,1,3,4,7)
p <- c(0.28,0.35,0.16,0.12,0.09)
sim41N <- function(n=10000){
  X<-numeric(n)
  for (i in 1:n) {
    U <- runif(1)
    if(U>0.37){if(U>0.72){X[i]<-0}
                else{X[i]<-1}}
    else{if(U>0.21){X[i]<-3}
          else{if(U>0.12){X[i]<-7}
                else{X[i]<-4}}}
    
  }
  return(X)
}
system.time(u <- sim41N(nsim))
##    user  system elapsed 
##    0.06    0.00    0.06

system.time(v <- rfmp.tabla(s,p,4,nsim))
##    user  system elapsed 
##    0.03    0.00    0.03

## Nº medio de comparaciones =  1

La eficiencia entre uno y otro algoritmo es muy simular, estudiando los tiempos de computo vemos que la diferencia entre el método usando arboles binarios y el método usando Tablas Guía, es apenas de un 0.01s, aun asi es claro que el método de Tablas Guía es más eficiente.

Ejercicio 5.

Considerese la variable aleatoria discreta, X, número de fallos mensuales en el suministro electrico de un equipo informatico, que tiene la siguiente masa de probabilidad:

x 0 1 2 3 4 5 6 7
P(X=x) 0.35 0.20 0.16 0.12 0.08 0.05 0.03 0.01

Construir una tabla guía y, basandose en ella, dar un algoritmo para simular esta variable. ¿Cuál es el número medio exacto de compraciones para generar un valor mediante ese algoritmo?

s <- c(0:7)
p <- c(0.35,0.20,0.16,0.12,0.08,0.05,0.03,0.01) 

ncomp <<- 0
rfmp.tabla <- function(x, prob = 1/length(x), m, nsim = 1000) {
  # Simulación v.a. discreta a partir de función de masa de probabilidad
  # por tabla guia de tamaño m
  # Inicializar tabla y FD
  Fx <- cumsum(prob)
  g <- rep(1,m)
  i <- 1
  for(j in 2:m) {
    while (Fx[i] < (j-1)/m) i <- i+1
    g[j] <- i
  }
  # Generar valores
  X <- numeric(nsim)
  U <- runif(nsim)
  for(j in 1:nsim) {
    i <- g[floor(U[j]*m)+1]
    while (Fx[i] < U[j]) i <- i + 1
    X[j] <- x[i]
    ncomp <<- ncomp+1
  }
  return(X)
}
system.time(v <- rfmp.tabla(s,p,7,nsim))
##    user  system elapsed 
##    0.03    0.00    0.03

## Nº medio de comparaciones =  1

Ejercicio 6.

Dado por supuesto un algoritmo que permita simular la variable aleatoria del ejercicio anterior, deducir un método para aproximar por simulación la probabilidad de que en un año haya mas de 30 fallos en el suministro electrico. Plantear el calculo analitico de la cantidad anterior.¿Sería sencillo obtener teoricamente la masa de probabilidad de la variable aleatoria número de falllos anuales?.

  1. Generar 12 muestras (1 año, doce meses) \(X_{i}\) de tamaño n.

  2. Hacer \(X=\sum_{1}^{12} X_{i}\).

  3. Obtener el número m, \(m \leq n\) de posibles resultados \(X[j] \geq 30\), con \(j \in [1,n]\)

  4. Devolver \(p=\frac{m}{n}\).

Proponemos el siguiente algoritmo:

set.seed(123)
n <- 12
nsim <- 1000
X<<-numeric(nsim)
for (i in 1:n) {
  u <-  rfmp.tabla(s,p,7,nsim)
  X <<- X+u
}
hist(X,freq = FALSE)

#Como vemos, la variable aleatoria numero de fallos por año para seguir una distribucion 
# N(nsigma, nmu)
probabilidad <- function(X,nsim){
  aux1 <-0
  for (i in 1:nsim) {
    if(X[i]>=30){aux1 <- aux1+1}
    
  }
  return(aux1/nsim)
}

probabilidad(X,nsim)
## [1] 0.07

Teoricamente podemos encontrar la masa de probabilidad de la variable aleatoria numero de fallos anuales. Modificando el algoritmo inicial dado, podremos obtener la probabilidad de que hayan de 0 a 84 fallos.

Ejercicio 7.

Se pretende simular \(nsim=10^{4}\) observaciones de una variable hipergeométrica (dhyper(x, m, n, k)) de parámetros m=el nº de grupo multiplicado por 10, \(n=100*m\) y \(k=20\)

1. Comprobar que el rango de posibles valores de esta variable es \(\max(0,k-n): \min(m,k)\). Generar los valores empleando el método de la transformación cuantil usando búsqueda secuencial. Obtener el tiempo de CPU empleado. Aproximar por simulación la función de masa de probabilidad, representarla gráficamente y compararla con la teórica. Calcular también la media muestral (compararla con la teórica \(km/(m+n)\)) y el número medio de comparaciones para generar cada observación.

Si los valores de la variable estan fuera de ese intervalo, tenemos que \(phyper = 0\)

#===============================================================================
# Distribucion Hipergeometrica
#===============================================================================
system.time(u<-rfmp(s,p,nsim))
##    user  system elapsed 
##    0.03    0.00    0.03

## Simulacion de la media muestral =  11.3939
##  Media Teorica =  11.42857

Esta bien aproximado.

2. Repetir el apartado anterior ordenando previamente las probabilidades en orden decreciente y también: empleando la función sample de R, mediante una tabla guía (con \(k-1\) subintervalos) y usando el método de Alias.

#===============================================================================
# Ordenando las probabilidades
#===============================================================================

ind <- order(vectorp(m,n,k), decreasing=TRUE)
system.time(rx <- rfmp(s[ind],p[ind],nsim))
##    user  system elapsed 
##       0       0       0
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("rx", "psim")
plot(as.matrix(res), col='green',lwd=3,type="h")
points(c(0:k),vectorp(m,n,k),pch=18)

#===============================================================================
# Usando la funcion Sample
#===============================================================================

system.time( rx <- sample(s, nsim, replace = TRUE, prob = p) )
##    user  system elapsed 
##       0       0       0
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("rx", "psim")
plot(as.matrix(res), col='brown',lwd=3,type="h")
points(c(0:k),vectorp(m,n,k),pch=16)

#===============================================================================
# Usando Tablas Guia
#===============================================================================
system.time( rx <- rfmp.tabla(s, p, k-1, nsim) )
##    user  system elapsed 
##    0.03    0.00    0.03
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("rx", "psim")
plot(as.matrix(res), col='violet',lwd=3,type="h")
points(c(0:k),vectorp(m,n,k),pch=13)

#===============================================================================
#Usando el metodo de Alias
#===============================================================================
system.time( rx <- rfmp.alias(s,p,nsim) )
##    user  system elapsed 
##    0.04    0.00    0.05
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("rx", "psim")
plot(as.matrix(res), col='orange',lwd=3,type="h")
points(c(0:k),vectorp(m,n,k),pch=9)

Métodos especificos para la simulación de variables discretas.

** Proponemos los siguientes algoritmos **

#=========================================================================
# Distribucion de Beornulli
#=========================================================================

rbeornulli <- function(p){
  U <- runif(1)
  if(U<=p){return(X<-1)}
  else{return(X<-0)}
}

rbeornulliN <- function(p,n=10000){
  X<-numeric(n)
  for (i in 1:n) {
    X[i]<-rbeornulli(p)
  }
  return(X)
}
system.time( rx <- rbeornulliN(0.5) )
##    user  system elapsed 
##    0.06    0.00    0.09
hist(rx)

#res <- as.data.frame(table(rx)/nsim)
#names(res) <- c("rx", "psim")
#plot(as.matrix(res), col='orange',lwd=3,type="h")
#points(c(0,1),c(0.5,0.5),pch=9)
#=========================================================================
# Distribucion Hipergeometrica
#=========================================================================
#=======================================================================
# Metodo de busqueda secuencial
#=======================================================================
ncomp <<- 0
rfmp <- function(x, prob = 1/length(x), nsim = 1000) {
  # Simulación nsim v.a. discreta a partir de fmp
  # por inversión generalizada (transformación cuantil)
  # Inicializar FD
  Fx <- cumsum(prob)
  # Simular
  X <- numeric(nsim)
  U <- runif(nsim)
  for(j in 1:nsim) {
    i <- 1
    while (Fx[i] < U[j]) i <- i + 1
    X[j] <- x[i]
    ncomp <<- ncomp + i
  }
  return(X)
}
#======================================================================
#Vector de probabilidades
vectorp <-function(m,n,k){
  s <- c(0:k)
  p <- (choose(m,s)*choose(n,k-s))/choose(m+n,k)
  return(p)
}


Hgeom <- function(m,n,k,nsim=10^4){
  s <- c(0:k)
  p <- (choose(m,s)*choose(n,k-s))/choose(m+n,k)
  X <- rfmp(s,p,nsim)
  return(X)
}