Simulacion de Variables Continuas y Discretas

Simulacion Variables Aleatorias Continuas (Seccion 3.3)

Ejercicio 1

Dar un algoritmo para simular \[f(x)= (\frac{2x}{e^x})^2 ,\, si \,\,\, x\geq 0 \] utilizando , de forma auxiliar , el generador de una exponencial de parametro 1 , comentar su eficiencia.

Solucion

Simularemos por aceptacion rechazo , donde la densidad auxiliar sera :\(g(x)= e^{-x}\)

c_opt = optimize(f=function(x){(2*x/exp(x))^2/dexp(x,1)}, maximum=TRUE, interval=c(0,4))
c.opt <- c_opt$objective
c.opt

## [1] 2.165365

lambda.opt <- 1
ngen <- 0

rnormAR <- function() {
  # Simulacion por aceptacion-rechazo
  while (TRUE) {
    U <- runif(1)
    X <- rexp(1)
    ngen <<- ngen+1 
    if (c.opt * U * dexp(X, lambda.opt) <= (2*X/exp(X))^2) return(X)
  }
}
rnormARn <- function(n=1000) {
# Simulacion n valores N(0,1)
    x <- numeric(n)
    for(i in 1:n) x[i]<-rnormAR()
    return(x)
}


# Grafico
curve(c.opt * dexp(x), xlim = c(0, 4), lty = 2, col="green")
curve((2*x/exp(x))^2, add = TRUE)

system.time(x <- rnormARn(10^4))

##    user  system elapsed 
##    0.14    0.03    0.17

hist(x, breaks="FD", freq=FALSE)
curve((2*x/exp(x))^2, add=TRUE, col="red")

Conclusiones:

Se puede decir que el algoritmo es eficiente pues tan solo demora tan solo 0.18 segundos en la ejecucion.
En la primera grafica observamos el comportamiento de nuestra funcion auxiliar\(g(x)\) multiplicada por el valor optimo 2,12653565 en el el intervalo [0,4], frente a la densidad teorica \(f(x)\) ,las graficas se encuentan de color verde y negro respectivamente, con esto, podemosnotar que la funcion auxiliar fue una buena eleccion.
Con el metodo de aceptacion rechazo y la ayuda de la funcion auxiliar \(g(x)= e^{-x}\) vemos que se tiene una buena simulacion de \(f(x)\) para el intervalo [0,4],generando 1000 valores aleatorios, se puede ver en la segunda grafica que la densidad teorica coloreada de rojo es bien aproximada.

Ejercicio 2

Dada una variable aleatoria continua con funcion de densidad:

\[ f=\begin{cases} 6x(1-x)& \text{ si } x\in[0,1]\\ 0& \text{ en otro caso } \end{cases} \] Dar un algoritmo para simular valores de la misma. Analizar la eficiencia de dicho algoritmo.

Solucion

Para la implementacion de este ejercicio utilizaremos el metodo de Aceptacion - Rechazo

ddensidad<- function(x){
  6*x*(1-x)
}
#Dado que en el intervalo [0,1] la funcion de densidad se asemeja a una parabola centrada en 0.5; utilizaremos como densidad auxiliar una funcion normal.
dauxiliar<- function(x){
  dnorm(x,0.5,0.265)
}
#Buscamos el mejor valor para c tq f(x)<=c*g(x)
c_opt = optimize(f=function(x){ddensidad(x)/dauxiliar(x)}, maximum=TRUE, interval=c(0,1))
c.opt <- c_opt$objective
c.opt

## [1] 1.221112

ngen<-0
ejercicio2<-function(){ 
 while (TRUE) {
    U <- runif(1)
    X <- rnorm(1,0.5,0.265) # generaremos segun la densidad auxiliar
    ngen <<- ngen+1 
    if (c.opt * U * dauxiliar(X) <= ddensidad(X)) return(X)
  }
}
sim2<- function(n=1000) {
# Simulacion n valores para la función
    x <- numeric(n)
    for(i in 1:n) x[i]<-ejercicio2()
    return(x)
}

curve(c.opt * dauxiliar(x), xlim = c(0, 1), lty = 2, col="green")
curve(ddensidad(x), add = TRUE)

nsim <- 10^4 #numero de simulaciones
system.time(x <- sim2(nsim))

##    user  system elapsed 
##    0.14    0.00    0.14

hist(x, breaks="FD", freq=FALSE)
curve(ddensidad(x), add=TRUE, col="red")

Conclusiones:

Dado que en las graficas se observa una aproximacion muy cercana entre las densidades, y se obtiene los histogramas de lo simulado con c= 1.2211, entonces la probabilidad de aceptacion en el paso es p=1/c = 0.8189, concluyendo que el algoritmo es bastante eficiente.

Ejercicio 3

Al bombardear una lamina circular de radio 1cm, echa de plata con particulas , la distancia de cada impacto al centro del circulo resulta ser una variable aleatoria continua con funcion de densidad dada por: \[ f(x)= 3x^2 \,\,\, , si \,\, 0\leq x\leq 1\] Detallar un algoritmo para simular la distancia al centro de la lamina en sucesivos impactos.

Solucion

Resolvemos por el metodo de Inversión, considerando como funcion de distribucion asociada a nuestra funcion de densidad: \[F(x)=x^{3}\]

ddensidad<- function(x){
  3*(x^2)
}

rddistrib <- function(){
  U <- runif(1)
    return(U^(1/3))
}

rddistribn <- function(n = 1000) {
  x <- numeric(n)
  for(i in 1:n) x[i]<-rddistrib()
  return(x)
}

set.seed(54321)
system.time(x <- rddistribn(10^3))

##    user  system elapsed 
##    0.02    0.00    0.01

hist(x, breaks = "FD", freq = FALSE)
 curve(ddensidad(x), add = TRUE,col="red")

Conclusiones

Para la densidad \(f(x)\) propuesta , resulta sencillo hacer la simulacion por el metodo de inversion, la inversa de esta funcion es facil de hallar, tambien vemos que la implementacion es eficiente, pues su timepo de ejecucion es de 0.2 segundos.
Para comparar la densidad teorica con la de la muestra tomada ; realizamos 1000 simulaciones, en el intervalo [0,1] y podemos apreciar que la simulacion de la densidad \(f(x)\) llega a tener un comportamiento muy simular ala densidad teorica , graficada en color verde.

Ejercicio 4

El tiempo de respuesta de un servidor informático es una variable con funcion de distribucion dada por: \[F(x)=1-(x+1)\exp^{-x}, \ \ si \ \ x\geq0\] y cero en otro caso. Dar un algoritmo, lo mas detallado posible, para simular valores de dicho tiempo de respuesta.

Solucion

La resolucion se realizara utilizando el metodo de Aceptacion Rechazo, considerando como densidad auxiliar a la distribucion exponencial de parametro 1.

ddensidad<- function(x){
  exp(-x)*x
}
#Como se busca acotar mi funcion, tomaremos como densidad auxiliar a la asociada a una distribución exponencial con parametro 1
dauxiliar<- function(x){
  dexp(x,1)
}
c_opt = optimize(f=function(x){ddensidad(x)/dauxiliar(x)}, maximum=TRUE, interval=c(0,10))
c.opt <- c_opt$objective
c.opt

## [1] 9.999944

#Notar que dependera del intervalo, por lo cual se suguiere tomar otra distribucion auxiliar, sin embargo por motivos de aprendizaje, continuaremos con la distribucion exponencial y analizaremos los resultados de la simulacion. 
ejercicio4<-function(){ 
  while(TRUE){ 
    u<- runif(1) 
    y<- rexp(1,1)# generaremos segun la densidad auxiliar
    if(u<= y) return(y)
  }
}
sim4 <- function(n){
  replicate(n,ejercicio4())
} 

curve(c.opt * dauxiliar(x), xlim = c(0, 10), lty = 2, col="green")
curve(ddensidad(x), add = TRUE)

nsim <- 10^4 #numero de simulaciones

system.time(x <- sim4(nsim))

##    user  system elapsed 
##    0.08    0.00    0.07

hist(x, breaks="FD", freq=FALSE)
curve(ddensidad(x), add=TRUE, col="red")

Conclusiones:

Como se menciono anteriormente, la simulacion de esta distribucion no es tan certera, en efecto, se puede observar en el histograma que existen muchos valores sobre la densidad teorica. Asi, se propone buscar otra densidad auxiliar tal que f(x)/g(x) este acotada por un valor c.
Recordemos que nuestro c dependia del intervalo, así la probabilidad de aceptacion en el paso es p=1/c = 1/max(x)=0 con x suficientemente grande.

Ejercicio 5

La densidad de probabilidad de una variable aleatoria viene dada por

\[ f(x) = \left \{ \begin{matrix} (2x^2+1)e^{-2x} & \mbox{si }x \geq 0 \\ 0 & \mbox{caso contrario }\end{matrix}\right. \] Encontrar un algoritmo para simular valores precedentes de dicha distribucion y comentar su eficiencia.

Solucion

Para simular esta distribucion , lo haremos por aceptacion rechazo tomando como funcion de densidad auxiliar \(g(x)= 2e^{-2x}\) la correspondiente a una variable aletoria \(X\sim exp(2)\) .

c_opt <- optimize(f=function(x){(2*(x^2)+1)*exp(-2*x)/dexp(x,2)}, maximum =TRUE, interval=c(0,6))

c.opt <- c_opt$objective
c.opt

## [1] 36.49937

ngen <- 0

rnormAR <- function() {
  # Simulacion por aceptacion-rechazo
  while (TRUE) {
    
    U <- runif(1)
    X <- rexp(1,2)
    #ngen <<- ngen+1 
    if((c.opt* U * dexp(X,2)) <= (2*(X^2)+1)*exp(-2*X)) return(X)
  }
}

system.time(X<- replicate(1000,rnormAR()))

##    user  system elapsed 
##    0.21    0.00    0.22

# Grafico
curve(c.opt * dexp(x,2), xlim = c(0,6), lty = 2)
curve((2*(x^2)+1)*exp(-2*x), add = TRUE)

hist(X, breaks="FD", freq=FALSE, )
curve((2*(x^2)+1)*exp(-2*x), add=TRUE, col="red")

Conclusiones

Se puede notar que el algoritmo implementado por el metodo de aceptacion rechazo para esta funcion es eficiente pues tiene un tiempo de ejecucion de 0.25 segundos.
La funcion auxiliar \(g(x)\) tomada fue adecuada, pues para una generacion de 1000 variables aleatorias, se tiene un comportamiento similar a la densidad teorica graficada en color verde, en el segundo grafico.

Ejercicio 6

Dada la funcion de densidad \[f(x)=\frac{x^{3}-12x+20}{48}, \ \ si \ \ x\in [0,4]\] y cero en el resto, detallar un algoritmo que permita simular valores de la misma.

Solucion

ddensidad<- function(x){
  (x^3-12*x+20)/48
}
# Dado que el intervalo esta entre 0 y 4, se propone tomar como desidad auxiliar, la densidad asociada a una distribucion uniforme en [0,4]
dauxiliar<- function(x){
  dunif(x,0,4)
}
c_opt = optimize(f=function(x){ddensidad(x)/dauxiliar(x)}, maximum=TRUE, interval=c(0,4))
c.opt <- c_opt$objective
# valor teorico c.opt <- 2
c.opt

## [1] 2.999819

ejercicio6<-function(){ 
  while(TRUE){ 
    u<- runif(1) 
    y<- runif(1,0,4)# generaremos segun la densidad auxiliar 
    if(u<= (y^3-12*y+20)/24) return(y)
  }
}
sim6 <- function(n){
  replicate(n,ejercicio6())
} 


curve(c.opt * dauxiliar(x), xlim = c(0,4), lty = 2, col="green")
curve(ddensidad(x), add = TRUE)

#Notar que las curvas se encuentran distanciadas, esto se debe a que el valor c calculado no se acerca al teórico y por lo tanto se suguiere buscar una mejor densidad aproximada.

nsim <- 10^4 #numero de simulaciones

system.time(x <- sim6(nsim))

##    user  system elapsed 
##    0.09    0.00    0.10

hist(x, breaks="FD", freq=FALSE)
curve(ddensidad(x), add=TRUE, col="red")

Conclusiones:

Dado que en la gráfica del histograma de los valores simulados, se acercan a la densidad f(x) se puede afirmar que el algoritmo da una buena simulacion, sin embargo no es la optima en efecto la probabilidad de aceptacion en el paso es: 1/c = 0.333353.

Ejercicio 7

El tiempo (en anios) entre dos inspecciones fiscales a una emrpesa es una variable aleatoria con densidad dada por

\[ f(x) = \left \{ \begin{matrix} 0.11 -0.0018x -0.00003x^2 & \mbox{si }x \in[0,10] \\ 0 & \mbox{caso contrario }\end{matrix}\right. \]

Utilizando un generador congruencial con parametros \(a=5\) ,\(c=33\) y \(m_0=1024\) y semilla \(x_0\) =27, simular tres tiempos entre inspecciones, por medio del algoritmo encontrado. Describir un metodo para simular el numero de inspecciones realizadas desde una dada hasta dentro de veinte anios.

Solucion

Dar un algoritmo lo mas sencillo posible que permita simular valores de esa variable. Comentar el grado de eficiencia de dicho algoritmo.

# densidad objetivo: f(x)
# densidad auxiliar: densidad e una uniforme en [0.10]

c_opt = optimize(f=function(x){(0.11-0.0018*x-0.00003*x^2)/dunif(x,0,10)}, maximum=TRUE, interval=c(0,10))

c.opt

## [1] 2.999819

c.opt <- c_opt$objective

ngen <- 0

rnormAR <- function() {
  # Simulacion por aceptacion-rechazo
  while (TRUE) {
    U <- runif(1)
    X <- runif(1,0,10)
    ngen <<- ngen+1
  
    if (c.opt * U * 1/10 <= (0.11-0.0018*X-0.00003*X^2)) return(X)
  }
}


 system.time( x<-replicate(1000,rnormAR()))

##    user  system elapsed 
##    0.01    0.00    0.01

# Grafico

curve(c.opt *dunif(x,0,10) , xlim = c(0, 10), lty = 2)
curve(0.11-0.0018*x-0.00003*x^2, add=TRUE,xlim = c(0, 10),col="green")

hist(x,  freq=FALSE,xlim = c(0, 10))
curve(0.11-0.0018*x-0.00003*x^2, add=TRUE,xlim = c(0, 10), col="green")

Conclusiones

Podemos ver que por el metodo de aceptacion rechazo , la aproximacion no es precisa, pero es bastante cercana, en la primera grafica notamos que la densidad auxiliar tomada para simular nuestra densidad objetivo\(f(x)\) es lo suficientem mente cercana, además en el segundo gráfico , podemos notar la densidad teorica con color verde, vemos que el comportamiendo es bastante cercano.
El timepo de ejecucion es de 0.04 segundos, por lo que podriamos decir que es eficiente.

Utilizando un generador congruencial con parametros \(a=5\) ,\(c=33\) y \(m_0=1024\) y semilla \(x_0\) =27, simular tres tiempos entre inspecciones, por medio del algoritmo encontrado. Describir un m?todo para simular el numero de inspecciones realizadas desde una dada hasta dentro de veinte anios.

initRANDC <- function(semilla=as.numeric(Sys.time()), a=5, c=33, m=1024) {
  .semilla <<- as.double(semilla) %% m  #Calculos en doble precision
  .a <<- a
  .c <<- c
  .m <<- m
  return(invisible(list(semilla=.semilla,a=.a,c=.c,m=.m))) #print(initRANDC())
}

# --------------------------------------------------
# RANDC()
#   Genera un valor pseudoaleatorio con el generador congruencial
#   Actualiza la semilla (si no existe llama a initRANDC)
RANDC <- function() {
    if (!exists(".semilla", envir=globalenv())) initRANDC()
    .semilla <<- (.a * .semilla + .c) %% .m
    return(.semilla/.m)
}

# --------------------------------------------------
# RANDCN(n)
#   Genera un vector de valores pseudoaleatorios con el generador congruencial
#   (por defecto de dimension 1000)
#   Actualiza la semilla (si no existe llama a initRANDC)
RANDCN <- function(n=1000) {
    x <- numeric(n)
    for(i in 1:n) x[i]<-RANDC()
    return(x)
    # return(replicate(n,RANDC()))  # Alternativa mas rapida    
}

Para simular 3 inspecciones tenemos:

x<-RANDCN(3)
x

## [1] 0.8702465 0.3834593 0.9495228

barplot(x)

Para simular las inspecciones dentro de 20 anios tenemos:

y<-RANDCN(20)
y

##  [1] 0.779840683 0.931429976 0.689376441 0.479108767 0.427770396
##  [6] 0.171078541 0.887619269 0.470322908 0.383841102 0.951432073
## [11] 0.789386925 0.979161189 0.928032507 0.672389098 0.394172055
## [16] 0.003086836 0.047660743 0.270530277 0.384877947 0.956616298

barplot(y)

Ejercicio 8

El tiempo, en milisegundos de anticipación (o retraso) de una señal respecto de otra en un sistema de comunicaciones puede considerarse una variable aleotoria con función de distribucion dada por \[F(x)=\frac{exp^{x}}{1+exp^{x}}\] Dar un algoritmo, lo mas detallado y simple posible, para simular valores de dicho tiempo de respuesta. Comparar la eficiencia computacional del metodo con el algoritmo clásico para simular la distribución exponencial.

Solucion

Implementaremos un algoritmo basado en el metodo de Inversion, para lo cual necesitaremos la funcion de densidad para verificar el correcto funcionamiento de la simulacion.

ddensidad <- function(x){
  exp(x)/(1+exp(x))^2
}

rdejer8 <- function(){
  U <- runif(1)
  
    return(log(-U/(U-1)))
}


rdejer8n <- function(n = 1000) {
    x <- numeric(n)
    for(i in 1:n) x[i]<-rdejer8()
    return(x)
}
nsim <- 10^4 #numero de simulaciones

system.time(x <- rdejer8n(nsim))

##    user  system elapsed 
##    0.03    0.00    0.03

hist(x, breaks = "FD", freq = FALSE)
curve(ddensidad(x), add = TRUE, col="red")

Concluciones:

El algoritmo se ajusta de manera adecuada y esto se puede evidenciar en el grafico del histograma comparando con la funcion de densidad calculada.

Ejercicio 9

El numero de horas diarias durante las cuales un terminal (de tipo A) de un laboratorio esta siendo utilizado es una variable con funcion de distribucion:

\[ F(x) = \left \{ \begin{matrix} \frac{e^{x-2}}{1+e^{x-2}} & \mbox{si }x \geq 0 \\ 0 & \mbox{en otro caso }\end{matrix}\right. \] Se pide :

a. Dar un algoritmo para simular dicha variable

Solucion

Resolveremos el problema por el metodo de Inversion, primero hallaremos la funcion inversa de \(F(x)\) \[ F^{-1}(u) = 2-log(\frac{-(u-1)}{u}) \]

densidad_fx <- function(x)
{
 exp(x+2)/(exp(x)+exp(2))^2 
}

rdist <- function(){
# Simulaci?n por inversi?n
# Doble exponencial
  U <- runif(1)
  x<- 2-log((1-U)/U)
  return(x)
}


rdistn <- function(n = 1000) {
# Simulaci?n n valores de doble exponencial
    x <- numeric(n)
    for(i in 1:n) x[i]<-rdist()
    return(x)
}

system.time(x<- rdistn (1000))

##    user  system elapsed 
##    0.02    0.00    0.02

hist(x, breaks = "FD", freq = FALSE)
curve(densidad_fx(x), add = TRUE, col="green")

Conclusiones

El tiempo de ejecucion es de 0.01 segundos, Resulta sencillo simular esta distribucion por el metodo de inversion, la aproximacion que resulta a ladensidad teorica graficada en color verde ,es muy cercana.

b. Si para otro tipo de terminales de tipo B el timepo de utilizacion tiene distribucion \(exp(\frac{1}{2})\) , detallar un algoritmo para aproximar, por simulacion , la probabilidad de que un terminal de tipo B sea mas utilizado que uno de Tipo A.

Solucion

Dado que la simulacion para una fucion exponencial ya se encuentra realizada en R basta hacer lo siguiente

a <- rdistn(1000)
b<- rexp(10000,1/2)
 indice <-(a<b)
  mean (indice)

## [1] 0.4409

la variable indice, cuenta las veces que en 1000 generaciones de valores aleatorios, la expresion booleana resulto verdadera, es decir las veces que el teminal de tipo B fue mas utilizado que la de A, con la funcion mean, hallamos la media que correspodne a la probabilidad deseada.

Simulacion Variables Aleatorias Discretas (Seccion 4.5)

Ejercicio 1

Cuando se produce una sobre carga electrica en una cadena productiva el numero de componentes averiadas esta descrito mediante una variable aleatoria, X, que puede tomar los valores 1,2,3,4 y 5, con probabilidades: \[p_1 =0.11, \,\, p_2=0.22 , \,\, p_3= 0.45 , \,\, p_4 = 0.13 , \,\,p_5=0.09\] Para simular dichas variables se usa el algoritmo: arbol binario

¿cual es el numero medio de comparaciones para simular un valor de X?

El numero medio de comparaciones \(N\) es \(E[N]\)

EN <- 1*0.11+2*0.22+3*0.45+4*0.13+5*0.09
EN

## [1] 2.87

Construir un algoritmo basado en un arbol binario optimo. ¿Cual es ahora el numero medio de comparaciones?

Solucion

Mostramos el siguiente codigo, donde se implementa el algoritmo planteado en el ejercicio.

arbol_binario <-function ()
{
  u<- runif(1)
  
  if(u>0.34)
    {   if (u>0.59){x=3}
        else{x=2}
    }
  else if(u>0.21)  {x=4}
   else if(u>0.11) {x=5} else {x=1}
 
  return (x)
  
}

arbinn <- function(n = 1000) {

    x <- numeric(n)
    for(i in 1:n) x[i]<-arbol_binario()
    return(x)
}

system.time(x<-arbinn())

##    user  system elapsed 
##    0.01    0.00    0.01

hist(x, freq = FALSE)
lines(density(x),col="red",lwd=1)

Conclusiones

El tiempo de ejecucion de este algoritmo es 0.03 segundos, para una muestra de 1000 valores, la simulacion fue bastante cercana la funcion de probabilidad teorica.

Luego el numero de comparaciones medio para encontrar una variable es:

ENb <- 1*0.45+2*0.22+3*0.13+4*0.11+5*0.09
ENb

## [1] 2.17

Se puede ver que el numero de comparaciones medio disminuye cuando se toma la estructura de un arbol binario pues este tiene las etiquetas prdenadas descendentemente.

Dar un algoritmo basado en el m?todo de la tabla guia.Calcular , tambien en este caso el numero medio de comparaciones.

x <- c(1 ,2 ,3, 4, 5)
fmp <- c(0.11,0.22,0.45,0.13,0.09)
n<-5
nsim<-1000
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 + i
  }
  return(X)
}

system.time( rx <- rfmp.tabla(x, fmp, n-1, nsim) )

##    user  system elapsed 
##    0.01    0.00    0.02

rx<-rfmp.tabla(x,prob=fmp,4,nsim = 1000)
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")

# Comparación teórica
plot(as.matrix(res), type="h")
points(x, fmp, pch=4)

luego, el numero medio de comparaciones es:

N<- ncomp/nsim
N

## [1] 5.642

Ademas podemos notar, con el grafico de la densidad teorica que la simulacion ene ste caso no es mu cercana especialmente, cuando x=1.

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: \[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)=0.04\] Dar un algoritmo para simular esta variable mediante búsqueda en forma de arbol 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.

x<-c(0,1,2,3,4,5,6)
prob<-c(0.33,0.22,0.15,0.11,0.09,0.06,0.04)
nsim<-1000

rfmp.bin<-function(nsim)
{
  U<-runif(nsim)
  X<-numeric(nsim)
  for(i in 1:nsim)
  {
    ifelse(U[i]>0.41,
           ifelse(U[i]>0.59,X[i]<-0,ifelse(U[i]>0.26,X[i]<-3,X[i]<-2)),
           ifelse(U[i]>0.19,X[i]<-1,ifelse(U[i]>0.1,X[i]<-4,ifelse(U[i]>0.04,X[i]<-5,X[i]<-6))))
    
  }
  return(X)
}
rx<-rfmp.bin(nsim)
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")
# Comparaci?n te?rica
plot(as.matrix(res), type="h")
points(x, prob, pch=4)

#B) ###
rfmp.tabla <- function(x, prob = 1/length(x), m, nsim = 1000) {
  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)
  Fx<-cumsum(prob)
  for(j in 1:nsim) {
    i <- g[floor(U[j]*m)+1]
    while (Fx[i] < U[j]) i <- i + 1
    X[j] <- x[i]
  }
  return(X)
}


rx<-rfmp.tabla(x,prob,6,nsim = 1000)
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")

# Comparación teórica
plot(as.matrix(res), type="h")
points(x, prob, pch=4)

Conclusiones

En referencia al ultimo grafico, se puede aseverar que la simulacion fue bastante cercana la funcion de probabilidad teorica, pues coinciden de manera cercana las marcas de las x con lo esperado.

Ejercicio 3

El numero 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\) y \(P(X=5)= 0.04\).

Dar un algoritmo para simular esta variable mediante busqueda en forma de arbol de Huffman.

Solucion

arbol binario

arbol_binario <-function ()
{
  u<- runif(1)
  
  if(u>0.40)
    {   if (u>0.79){x=1}
        else{x=2}
    }
  else if(u>0.21)  {x=3}
   else if(u>0.13) {x=4} 
      else if (u>0.09){x=5}
          else{x=0}
 
  return (x)
  
}

arbinn <- function(n = 1000) {

    x <- numeric(n)
    for(i in 1:n) x[i]<-arbol_binario()
    return(x)
}

system.time(x<-arbinn())

##    user  system elapsed 
##    0.02    0.00    0.02

x<-arbinn()

hist(x, freq = FALSE)
lines(density(x),col="red",lwd=1)

fmp <- c(0.09,0.21,0.39,0.19,0.08,0.04)
x1 <-c(0,1,2,3,4,5)
rx<-rfmp.tabla(x1,fmp,5,nsim = 1000)
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")

# Comparación teórica
plot(as.matrix(res), type="h")
points(x1, fmp, pch=5)

Calcular el numero medio de comparaciones para simular cada valor de la variable \(X\) .

EN <- 0*0.09+1*0.21+2*0.39+3*0.19+4*0.08+5*0.04
EN

## [1] 2.08

Conclusion

Podemos notar que el numero medio de comparaciones menor es de 2.17 ,el cual se otiene con la estructura de un arbol binario , esto es debido a que su estructura ya mantiene un orden que facilita las comparaciones.

¿Que ventajas se obtienen para este metodo de simualcion en comparacion con la implementacion mediante el metodo secuencial(tanto en su version directa como en la del etiquetado mas eficiente)?

n <- 6
nsim <- 1000

x <-c(0,1,2,3,4,5)
fmp <- c(0.09,0.21,0.39,0.19,0.08,0.04)

ncomp <- 0
rfmp <- function(x, prob = fmp, 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( rx <- rfmp(x, fmp, nsim) )

##    user  system elapsed 
##       0       0       0

rx<-rfmp.tabla(x,prob=fmp,5,nsim = 1000)
res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")

# Comparación teórica
plot(as.matrix(res), type="h")
points(x, fmp, pch=4)

# N?mero de comparaciones
ncomp/nsim

## [1] 3.08

Con esto vemos que la funcion de probabilidad simulada es bastante cercana a la probabilidad teorica, no varia mucho.

El numero de comparaciones sin etiquetas ordenadas es:

ncomp/nsim

## [1] 3.08

Ahora, usando las etiquetas ordenadas tenemos

n <- 6
nsim <- 1000

x <-c(0,1,2,3,4,5)
fmp <- c(0.39,0.21,0.19,0.09,0.08,0.04)
ncomp <- 0

 system.time( rx <- rfmp(x, fmp, nsim) )

##    user  system elapsed 
##       0       0       0

# N?mero de comparaciones
ncomp/nsim

## [1] 2.304

#graficas

res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")

# Comparación teórica
plot(as.matrix(res), type="h")
points(x, fmp, pch=4)

El numero de comparaciones es:

ncomp/nsim

## [1] 2.304

Conclusion

Claramente hay ventaja con respecto al numero medio de comparaciones .Cuando las etiquetas ya estan ordenadas el numero de comparaciones baja de 6 a 4, tambien el tiempo de ejecucion para el segundo caso es menor.
Ademas notemos que esta simulacion tiene un buen comportamiento, ya que es bastante parecida a la probabilidad teorica, los cuales son los valores punteados con una “x”en la grafica.

Ejercicio 4

El codigo de error que puede suministrar una rutina implementada en un lenguaje de alto nivel puede describirse mediante una variable aleatoria C, que puede tomar los valores 0,1,3,4 y 7, con probabilidades: \[p_{0}=0.28, \ \ p_{1}=0.35,\ \ p_{3}=0.16, \ \ p_{4}=0.12, \ \ p_{7}=0.09 \] 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 arbol de Huffman y en el método de las tablas guía. Estudiar la eficiencia de ambos algoritmos y compararla.

s <- c(0,1,3,4,7)
p <- c(0.28,0.35,0.16,0.12,0.09)

simN <- 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 <- simN(nsim))

##    user  system elapsed 
##    0.03    0.00    0.03

res <- as.data.frame(table(u)/nsim)
names(res) <- c("u", "psim")
plot(as.matrix(res),type="h")
points(x, fmp, pch=4)

ncomp <<- 0
rfmp.tabla <- function(x, prob = 1/length(x), m, nsim = 1000) {
  # Simulacion v.a. discreta a partir de funcion de masa de probabilidad
  # por tabla guia de tamanioo m
  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,4,nsim))

##    user  system elapsed 
##    0.02    0.00    0.01

## N medio de comparaciones =  1

Se puede notar que el metodo de las tablas Guia es mas eficiente, pese a eficiencia entre los algoritmos es muy simular, con una diferencia de apenas 0.001s en los tiempos de computo entre el metodo de arboles binarios y el metodo delas Tablas Guía.

Ejercicio 5

Considerese la variable aleatoria discreta, \(X\), numero de fallos mensuales en el suministro electrico de un equipo informatico, que tiene la siguiente masa de probabilidad.

\[\begin{array}{ | l | l | l | l | l | l | l | l | l | } \hline x & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ \hline P(X=x) & 0.35 & 0.20 & 0.16 & 0.12 & 0.08 & 0.05 & 0.03 & 0.01 \\ \hline \end{array}\]

Construir una tabla guia y, basandose en ella, dar un algoritmo para simular esta variable.¿Cual es el numero medio exacto de comparaciones para generar un valor mediante este algoritmo ?

Solucion

x<-c(0,1 ,2 ,3, 4, 5,6,7)
fmp <- c(0.35,0.20,0.16,0.12,0.08,0.05,0.03,0.01)
n<-7
nsim<-1000
ncomp <- 0


system.time( rx <- rfmp.tabla(x, fmp, n-1, nsim) )

##    user  system elapsed 
##       0       0       0

res <- as.data.frame(table(rx)/nsim)
names(res) <- c("x", "psim")

# Comparación teórica
plot(as.matrix(res), type="h")
points(x, fmp, pch=4)

luego, el numero medio de comparaciones es:

N<- ncomp/nsim
N

## [1] 1

Ejercicio 6

Dando por supuesto un algoritmo que permita simular la variable aleatoria del ejercicio anterior, deducir un metodo 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. ¿Seria sencillo obtener teoricamente la masa de probabilidad dela variable aleatoria numero de fallos anuales?

probabilidad <-function (){

x<-c(0,1 ,2 ,3, 4, 5,6,7)
fmp <- c(0.35,0.20,0.16,0.12,0.08,0.05,0.03,0.01)
nsim<-12
n<- 8
ncomp <- 0

 X<- rfmp.tabla (x, prob= fmp, m=n-1, nsim = 12)
 indice<-(sum(X)>30)
  return(indice)
}
muestra <-1000
pro <-replicate(muestra,probabilidad ())
probabilidad<-sum(pro)/muestra

probabilidad

## [1] 0.062

7.3 EJERCICIOS (Libro Ruben Casal)

Ejercicio 2.3

Escribir el codigo necesario para generar, por el motodo de inversi0n, una muestra de \(n\) observaciones de una distribucion de Cauchy.

Solucion

Simularemos valores aleatorios para una distribucion de Cauchi \(C \sim (x,0,1)\), cuya funci?n de distribucion inversa es es: \[F^{-1}(u)=tan(\pi(u- \frac{1}{2}))\]

Generar una muestra de \(10^4\) observaciones y obtener el tiempo de CPU

rdcauchy <- function(lambda = 1){
# Simulacion por inversion
#Cauchy
  U <- runif(1)
  
    return(tan(pi*(U-1/2)))
 
}

rdcauchyn <- function(n = 1000) {
# Simulacion n valores de cauchy
    x <- numeric(n)
    for(i in 1:n) x[i]<-rdcauchy()
    return(x)
}

nsim<- 10^4

system.time(x <- rdcauchyn(nsim))

##    user  system elapsed 
##    0.05    0.00    0.05

Representar el histograma (limitar el rango, e.g. xlim = c(-10, 10)) y compararlo con la densidad teorica (dcauchy).

x<- rdcauchyn(nsim)
hist(x, breaks = "FD", xlim= c(-10, 10), freq = FALSE)
# lines(density(x), col = 'blue')
curve(dcauchy(x,0,1),xlim=c(-10, 10), add = TRUE, col="red")

Obtener conclusiones sobre la existencia de una media teorica a partir de la media muestral aproximada por simulacion (estudiar la convergencia de la media muestral). Suponiendo que el vector x contiene las simulaciones, estudiar la convergencia de la media muestral mediante el grafico:plot(1:nsim, cumsum(x)/(1:nsim), type=“l”, ylab=“Media muestral”, xlab=“N de simulaciones”)

mean(x)

## [1] 1.85602

summary (dcauchy(x,0,1))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.04543 0.16132 0.15946 0.27232 0.31831

plot(1:nsim, cumsum(x)/(1:nsim), type="l", ylab="Media muestral", xlab=10^4)

Ejercicio 2.4

El tiempo de respuesta (en centésimas de segundo) de un servidor de bases de datos es una variable con función de densidad: \[f(x)=xe^{-x}, \ \ si \ \ x \geq 0.\] Escribir el código necesario para generar, por el método de aceptación-rechazo, una muestra de n observaciones de esta distribución empleando como densidad auxiliar una exponencial: \[g(x)=λe^{−λx}, \ \ si \ \ x \geq 0.\] ### Solucion a. Aproximar numéricamente el parámetro óptimo (\(λ_{opt} < 1\)) y la cota óptima (\(c_{opt}\)) de la densidad auxiliar y compararlos con los valores teóricos: \(λ_{opt}=1/2\) y \(c_{opt}=4/e\).

ddensidad <- function(x){
  x*exp(-x)
}
fopt <- function(lambda) {
  # Obtiene c fijado lambda
  optimize(f = function(x){ddensidad(x)/dexp(x,lambda)}, maximum=TRUE, interval=c(0,5))$objective
}

res <- optimize(f=function(x){fopt(x)}, interval=c(0,5))
lambda.opt2 <- res$minimum
c.opt2 <- res$objective

lambda.opt2

## [1] 0.5000098

# valor teórico λopt=1/2

c.opt2

## [1] 1.471518

#valor teórico copt=4/e

Generar una muestra de 1000 observaciones de la distribución de interés (tomando como semilla inicial el nº de grupo multiplicado por 100). Obtener el tiempo de CPU que tarda en generar la secuencia y calcular el número medio de generaciones de la distribución auxiliar.

nsim <- 10^3
ngrupo<-1
set.seed(ngrupo*100)
ngen <- 0
rnormAR <- function() {
  while (TRUE) {
    U <- runif(1)
    X <- rexp(1,lambda.opt2)
    ngen <<- ngen+1 
    if (c.opt2 * U * dexp(X, lambda.opt2) <= ddensidad(X)) return(X)
  }
}
rnormARn <- function(n=1000) {
# Simulacion n valores N(0,1)
    x <- numeric(n)
    for(i in 1:n) x[i]<-rnormAR()
    return(x)
}

ngen <- 0
system.time(x <- rnormARn(nsim))

##    user  system elapsed 
##    0.03    0.00    0.03

{
cat("\nNº de generaciones = ", ngen)
cat("\nNº medio de generaciones = ", ngen/nsim)
cat("\nProporción de rechazos = ", 1-nsim/ngen, "\n")
}

## 
## Nº de generaciones =  1504
## Nº medio de generaciones =  1.504
## Proporción de rechazos =  0.3351064

c.Representar el histograma y compararlo con la densidad teórica.

# Grafico

hist(x, breaks="FD", freq=FALSE)
curve(ddensidad(x), add=TRUE, col = "red")

Metodo Aceptacion - Rechazo(Distribucion Poisson de parametro λ = 3 a partir de una Distribucion Geometrica)

# DISTRIBUCION POISSON

f <- function ( x ){
(exp(1)^( -3 ))*( 3^ x / factorial ( x ))/( 2.13*((3/4)^ x )*1/4)
}
Alg <- function (){
while ( TRUE ){
u <- runif (2)
y <- floor (( log ( u [1 ]))/( log(3/4 )))
if( u [2] <= f (y )) return ( y )
}
}
Alg_2 <- function( n ) replicate (n , Alg ())
x <- Alg_2(10^4)
barplot ( table ( x )/ length ( x ))
y<-rpois(10^4,3)
barplot(table ( y )/ length ( y ),border = "red", add = TRUE)

Simulacion de una Normal mediante el metodo de Marsaglia

El metodo de Marsaglia es un metodo de muestreo de numeros pseudo aleatorios para generar una aproximacion de variables aleatorias normales estandar independientes. El metodo polar funciona seleccionando puntos aleatorios \(( x , y )\) en el cuadrado \(-1<x <1\) , \(-1 <y <1\) hasta \(0 < s = x ^ {2} + y ^ {2} <1\) y luego devolver el par requerido de variables aleatorias normales como:

\[x \sqrt{\frac{-2ln(s)}{s}}\, , \, y\sqrt{\frac{-2ln(s)}{s}}\] donde \(\frac{x}{\sqrt{s}}\) y \(\frac{y}{\sqrt{s}}\) representa el coseno y el seno del angulo que el vector \((x,y)\) hace con el eje x.

tenemos las siguiente implementacion:

Marsaglia <-function (){
u=runif(2)
v1<-2*u[1]-1
v2<-2*u[2]-1
s<-v1^{2}+v2^{2}
while(s>1){
  u=runif(2)
  v1<-2*u[1]-1
  v2<-2*u[2]-1
  s<-v1^{2}+v2^{2}
}
x=v1*sqrt(-2*log(s)/s)
y=v2*sqrt(-2*log(s)/s)
}
normal_Mars <- function(n) replicate(n,Marsaglia())
n=10^4
hist(normal_Mars(n),prob = T)
curve(dnorm(x), add=TRUE, col = "red")

Simulacion de Variables Continuas y Discretas

Ivan Solorzano y Genesis Moreno

2 de mayo de 2019

Simulacion Variables Aleatorias Continuas (Seccion 3.3)

Ejercicio 1

Solucion

Conclusiones:

Ejercicio 2

Solucion

Conclusiones:

Ejercicio 3

Solucion

Conclusiones

Ejercicio 4

Solucion

Conclusiones:

Ejercicio 5

Solucion

Conclusiones

Ejercicio 6

Solucion

Conclusiones:

Ejercicio 7

Solucion

Conclusiones

Ejercicio 8

Solucion

Concluciones:

Ejercicio 9

a. Dar un algoritmo para simular dicha variable

Solucion

Conclusiones

b. Si para otro tipo de terminales de tipo B el timepo de utilizacion tiene distribucion \(exp(\frac{1}{2})\) , detallar un algoritmo para aproximar, por simulacion , la probabilidad de que un terminal de tipo B sea mas utilizado que uno de Tipo A.

Solucion

Simulacion Variables Aleatorias Discretas (Seccion 4.5)

Ejercicio 1

Solucion

Conclusiones

Ejercicio 2

Conclusiones

Ejercicio 3

Solucion

Conclusion

Conclusion

Ejercicio 4

Ejercicio 5

Solucion

Ejercicio 6

7.3 EJERCICIOS (Libro Ruben Casal)

Ejercicio 2.3

Solucion

Ejercicio 2.4

Metodo Aceptacion - Rechazo(Distribucion Poisson de parametro λ = 3 a partir de una Distribucion Geometrica)

Simulacion de una Normal mediante el metodo de Marsaglia