Métodos universales para la generación de variables continuas

Ejercicio 1.

Dar un algoritmo para simular la variable aleatoria continua con densidad

\[\begin{eqnarray*} f(x)&=&\left(\frac{2x}{e^{x}}\right)^{2}, \text{ } si \text{ } x \geq 0. \end{eqnarray*}\]

utilizando, de forma auxiliar, el generador de una exponencial con parámetro uno. Comentar su eficiencia

Utilizando el método de aceptación-rechazo con funcion auxiliar la de una exponencial con \(\lambda=1\), comenzaremos encontrando el valor óptimo de \(c\), de tal forma que para todo \(x \geq 0\) se tiene:

\[\begin{align*} f(x) \leq c*dexp(x,\lambda=1) \end{align*}\]

#-----------------------------------------------------------------------------------
c <- optimize(f=function(x){f(x)/dexp(x)}, maximum=TRUE, interval=c(0,4))
c.opt <- c$objective
#Como observacion el intervalo donde la función optimize encuentra el máximo debe ser
#cuidadosamente selecionado en base a donde esta definida la funcion f.

## c óptimo =  2.165365

Asi pues el algotimo que proponemos basado en el método de aceptación-rechazo es el siguiente.

lambda.opt <- 1
fR <- function() {
  # Simulación por aceptación-rechazo
  # f a partir de una exponencial con parámetro lambda=1
  while (TRUE) {
    U <- runif(1)
    X <- rexp(1)
    ngen <<- ngen+1 
    if (c.opt * U * dexp(X, lambda.opt) <= f(X)) return(X)
  }
}

fRn <- function(n=1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-fR()
  return(x)
}
#-----------------------------------------------------------------------------------

Para evaluar el funcionamiento y la eficiencia del algoritmo propuesto comenzaremos generando una muestra de tamaño \(nsim=10^4\), e interprateremos el número medio de comparaciones que el algoritmo realiza.

system.time(v <- fRn(nsim))

##    user  system elapsed 
##    0.29    0.02    0.33

## Análisis del histograma de x.
hist(v, breaks="FD", freq=FALSE)
curve(f(x), col='blue',add=TRUE)

#-----------------------------------------------------------------------------------

## Nº de generaciones =  21816
## Nº medio de comparaciones =  2.1816
## Proporción de rechazos =  0.5416208

Notamos que el número medio de comparaciones es:

## 2.1816

Asi pues tenemos que no es considerablemente eficiente, ya que la eficiencia es mejor mientras este valor sea mas cercano a 1.
Para finalizar el analisis del funcionamiento de este algoritmo usaremos un test de Kolmogorov-Smirnov. Para ello comenzaremos determinando la función de distribución a partir de f.

\[\begin{align*} F(x) &= \int^{x}_{-\infty} f(t)dt\\ &= \int_{0}^{x} \left(\frac{2t}{e^t} \right)^2 dt\\ &= \int_{0}^{x} 4t^{2}e^{-2t} dt\\ &= 4\int_{0}^{x} t^{2}e^{-2t} dt\\ &\text{Tomando } u=t^{2} \text{ y } dv=e^{-2t}\\ &\text{Entnces } du=2tdt \text{ y } v=-\frac{e^{-2t}}{2}\\ &= 4\left( -\frac{t^{2} e^{-2t}}{2} -\int_{0}^{x} - e^{-2t} t dt\right)\\ &= 4\left( -\frac{t^{2} e^{-2t}}{2} +\int_{0}^{x} e^{-2t} t dt\right)\\ &\text{Tomando } u=t \text{ y } dv=e^{-2t}\\ &\text{Entnces } du=dt \text{ y } v=-\frac{e^{-2t}}{2}\\ &= 4\left( -\frac{t^{2} e^{-2t}}{2} +\left( -\frac{te^{-2t}}{2} +\int_{0}^{x} \frac{e^{-2t}}{2} dt\right)\right)\\ &= \left. 4\left( -\frac{t^{2} e^{-2t}}{2} +\left( -\frac{te^{-2t}}{2} -\frac{e^{-2t}}{4}\right)\right)\right|^{x}_{0}\\ &= \left. \left( -2t^{2} e^{-2t} +\left( -2te^{-2t} -e^{-2t}\right)\right)\right|^{x}_{0}\\ &= -2x^{2}e^{-2x} - 2x^{-2x}-e^{-2x}+1 \end{align*}\] Asi pues definimos: \[\begin{align*} F(x)= \left\{ \begin{array}{lcc} -2x^{2}e^{-2x} - 2x^{-2x}-e^{-2x}+1 & si & x \geq 0 \\ \\ 0 & \text{Caso contrario} & \\ \end{array} \right.. \end{align*}\]

#Test de Kolmogorov-Smirnov
ks.test(v,'dF')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  v
## D = 0.0071393, p-value = 0.6879
## alternative hypothesis: two-sided

Con lo cual, tenemos que en efecto el algoritmo funciona para generar variables cuya densidad esta dada por \(f\).

Ejercicio 2

Dada una variable aleatoria continua con función de densidad

\[\begin{align*} f(x)= \left\{ \begin{array}{lcc} 6x(1-x) & \text{si } x \in [0,1] \\ \\ 0 & \text{Caso contrario} & \\ \end{array} \right.. \end{align*}\]

dar un algoritmo para simular valores de la misma. Analizar la eficiencia de dicho algoritmo.

Utilizando el método de aceptación-rechazo con función auxiliar la de una U(0,1), comenzaremos encontrando el valor Óptimo de \(c\), de tal forma que para todo \(x \in [0,1]\) se tiene:

\[\begin{align*} f(x) \leq c*dunif(x,0,1) \end{align*}\]

#-----------------------------------------------------------------------------------
##Usaremos la funcion de distribucion Uniforme(0,1)
#-----------------------------------------------------------------------------------
c <- optimize(f=function(x){df(x)/dunif(x,0,1)}, maximum=TRUE, interval=c(0,1))
c.opt <- c$objective
#-----------------------------------------------------------------------------------

## c Óptimo =  1.5

Asi pues el algotimo que proponemos basado en el método de aceptación-rechazo es el siguiente.

#-----------------------------------------------------------------------------------
rfR <- function() {
  # Simulación por aceptación-rechazo
  # f a partir de la U(0,1)
  while (TRUE) {
    U <- runif(1)
    X <- runif(1,0,1)
    ngen <<- ngen+1 
    if (c.opt * U * dunif(X,0,1) <= df(X)) return(X)
  }
}

rfRn <- function(n=1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rfR()
  return(x)
}
#-----------------------------------------------------------------------------------

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

##    user  system elapsed 
##     0.2     0.0     0.2

#-----------------------------------------------------------------------------------

#Análisis del histograma de x
hist(x, breaks="FD", freq=FALSE)
curve(df(x), col='yellow',lwd=2,add=TRUE)

#-----------------------------------------------------------------------------------

## Nº de generaciones =  14873
## Nº medio de generaciones =  1.4873
## Proporción de rechazos =  0.3276407

Notamos que el número medio de comparaciones es:

## 1.4873

Asi pues tenemos que es considerablemente eficiente, ya que la eficiencia es mejor mientras este valor sea mas cercano a 1.

Para finalizar el analisis del funcionamiento de este algoritmo usaremos un test de Kolmogorov-Smirnov. Para ello comenzaremos determinando la funcion de distribucion a partir de f.

\[\begin{align*} F(x) &= \int_{0}^{x} f(t)dt\\ &= \int_{0}^{x} 6t(1-t) dt\\ &= 6 \int_{0}^{x}t -t^{2} dt\\ &= 6 \left.\left( \frac{t^{2}}{2}-\frac{t^{3}}{3}\right) \right|^{x}_{0}\\ &= \left.3t^{2} - 2t^{3} \right|^{x}_{0}\\ &= 3x^{2}-2x^{3} \end{align*}\]

Luego, definiendo:

\[\begin{align*} F(x)= \left\{ \begin{array}{lcc} 0 & \text{si } x < 0 \\ 3x^{2} - 2x^{3} & \text{si } x \in [0,1] \\ 1 & \text{si } x>1 \end{array} \right.. \end{align*}\]

# Kolmogorov-Smirnov Tests
ks.test(x,'pF')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.009042, p-value = 0.387
## alternative hypothesis: two-sided

#-----------------------------------------------------------------------------------

Con lo cual, tenemos que en efecto el algoritmo funciona para generar variables cuya densidad esta dada por \(f\).

Ejercicio 3

Al bombardear una lamina circular de radio 1cm, hecha de plata, con particulas \(\alpha\), la distancia de cada impacto al centro del círculo resulta ser una variable aleatoria continua con función de densidad dada por \(f(x)=3x^{2}\) , si \(0\leq x \leq 1\). Detallar un algoritmo, lo más sencillo posible, para simular la distancia al centro de la lámina en sucesivos impactos.

Para ello, primero determinamos la función \(F\) de distribución.

\[\begin{align*} F(x) &= \int_{0}^{x}f(t)dt\\ &= \int_{0}^{x} 3t^{2}dt\\ &= 3 \left.\left(\frac{t^{3}}{3}\right)\right|^{x}_{0}\\ &= x^{3} \end{align*}\]

Entonces, definimos:

\[\begin{align*} F(x)= \left\{ \begin{array}{lcc} 0 & \text{si } x < 0 \\ x^{3} & \text{si } x \in [0,1] \\ 1 & \text{si } x>1 \end{array} \right.. \end{align*}\]

Dado que la funcion \(F\) es invertible en el intervalo \([0,1]\), podemos hacer uso del método de inversión. Asi pues, el algoritmo que proponemos es el siguiente.

##### Usando el método de inversión. 
#-----------------------------------------------------------------------------------
rfR <- function(){
  U <- runif(1)
  return((U)^(1/3))
  
}
#-----------------------------------------------------------------------------------
rfRn <- function(n = 1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rfR()
  return(x)
}
#-----------------------------------------------------------------------------------

Para evaluar el funcionamiento y la eficiencia del algoritmo propuesto comenzaremos generando una muestra de tamaño \(nsim=10^4\), además realizaremos el contraste de Kolmogorov-Smirnov,

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

##    user  system elapsed 
##    0.06    0.00    0.07

#Tenemos que el algoritmo es eficiente.

hist(x, breaks = "FD", freq = FALSE)
curve(f(x),col='brown',lwd=2, add = TRUE)

#-----------------------------------------------------------------------------------
#Con lo cual tenemos que en efecto los datos generados siguen una distribucion F.

#-----------------------------------------------------------------------------------
#Kolmogorov-Smirnov Tests
ks.test(x,'pF')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.0085815, p-value = 0.453
## alternative hypothesis: two-sided

#-----------------------------------------------------------------------------------

Asi pues, aceptamos que las datos generados siguen una distribucion F.

Ejercicio 4.

El tiempo de respuesta en centesimas de segundo de un servidor infórmatico es una variable con función de distribución dada por \(F(x)=1-(x+1)e^{-x}\), si \(x \geq 0\) y cero en otro caso. Dar un ejemplo, lo más detallado posible, para simular valores de dicho tiempo de respuesta.

Dado que \(F\), no es invertible en \([0,\infty]\), usaremos el método de aceptación-rechazo, usando como funcion axuliar la de una gamma con parametros k=2, lambda=sc. Comenzaremos determinando la funcion de densidad f, a partir de F.

\[\begin{align*} f(x) &= \frac{d}{dx} F(x)\\ &= \frac{d}{dx} 1-(x+1)e^{-x}\\ &= -e^{-x} +(x+1)e^{-x}\\ &= -e^{-x} +xe^{-x}+e^{-x}\\ &= xe^{-x} \end{align*}\]

Entonces definimos:

\[\begin{align*} f(x)= \left\{ \begin{array}{lcc} xe^{-x} & \text{si } x \geq 0 \\ 0 & \text{caso contrario } \\ \end{array} \right.. \end{align*}\]

Ahora determinaremos los valores de lambda y c optimos tales que:

\[\begin{align*} f(x) \leq c*dgamma(2,lambda.opt) \end{align*}\]

#-----------------------------------------------------------------------------------
  fopt <- function(sc) {
    # Obtiene c fijado k
    optimize(f = function(x){df(x)/dgamma(x,2,sc)},
             maximum=TRUE, interval=c(0,10))$objective
  }
  
  # Encontar lambda que minimiza
  res <- optimize(f=function(x){fopt(x)}, interval=c(0,10))
  sc.opt <- res$minimum
  c.opt <- res$objective
#-----------------------------------------------------------------------------------

## c optimo = 1.000021
##  lamda optimo= 0.9999895

Asi pues el algotimo que proponemos basado en el método de aceptación-rechazo es el siguiente.

#===================================================================================
## Usaremos la distribución gamma, con parametros shape=2 , scale = sc.opt
#===================================================================================
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
rfR <- function() {
  # Simulación por aceptación-rechazo
  while (TRUE) {
    U <- runif(1)
    X <- rgamma(1,2,sc.opt)
    ngen <<- ngen+1 
    if (c.opt * U * dgamma(X,2,sc.opt) <= df(X)) return(X)
  }
}

rfRn <- function(n=1000) {
  # Simulación n valores N(0,1)
  x <- numeric(n)
  for(i in 1:n) x[i]<-rfR()
  return(x)
}
#-----------------------------------------------------------------------------------

Para evaluar el funcionamiento y la eficiencia del algoritmo propuesto comenzaremos generando una muestra de tamaño \(nsim=10^{4}\), e interprateremos el número medio de comparaciones que el algoritmo realiza.

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

##    user  system elapsed 
##    0.25    0.00    0.25

#-----------------------------------------------------------------------------------
# Analisis del histograma de x
hist(x, breaks="FD", freq=FALSE)
curve(df(x), col='blue', add=TRUE)

# Tenemos que x sigue una distribución F
#-----------------------------------------------------------------------------------

## Nº de generaciones =  10000
## Nº medio de generaciones =  1
## Proporción de rechazos =  0

Es intersante notar que virtualmente el algoritmo funciona perfectamente, pero debemos considerar que la distribucion \(gamma(k,\lambda)\), no es sencilla de simular, ademas que en base a la graficas parece que hemos caido en la solucion trivial de \(f=g\) al momento de buscar g. Lo cual ha sucedido por que al ver la forma de la curva de densidad de la funcion \(f\), esta era muy similar a la de una \(gamma(2,1)\). MÃ¡s adelante usaremos como funcion auxilar otra diferente a la gamma y veremos que sucede.

Para finalizar el análisis del funcionamiento del algoritmo realizaremos el test de Kolmogorov-Smirnov

#Kolmogorov-Smirnov Tests
ks.test(x,'pF')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.011377, p-value = 0.1502
## alternative hypothesis: two-sided

Con lo cual aceptamos, que los datos generados siguen una distribucion F.

Ejercicio 5.

La densidad de probabilidad de una variable aleatoria viene dada por.

\[\begin{align*} f(x)= \left\{ \begin{array}{lcc} (2x^{2} +1)e^{-2x} & \text{si } x \geq 0 \\ 0 & \text{caso contrario } \\ \end{array} \right.. \end{align*}\]

Encontrar un algoritmo para simular valores procedentes de discha distribución y comentar su eficiencia.

Comenzaremos definiendo la función F, puesto que la función de densidad en este ejercicio es una variante de la vista en el ejercicio 1, obviaremos los calculos para determinar F.

Asi pues, definimos: \[\begin{align*} F(x)= \left\{ \begin{array}{lcc} -x^{2}e^{-2x} - x^{-2x}-e^{-2x}+1 & Si & x \geq 0 \\ \\ 0 & \text{Caso contrario} & \\ \end{array} \right.. \end{align*}\]

#===================================================================================
# Usando como distribución auxiliar una auxiliar con parámetro lambda. 
#===================================================================================
#-----------------------------------------------------------------------------------

# Obtención de valores c y lambda óptimos aproximados
fopt <- function(lambda) {
  # Obtiene c fijado lambda
  optimize(f = function(x){df(x)/dexp(x,lambda)},
           maximum=TRUE, interval=c(0,10))$objective
}

# Encontar lambda que minimiza
res <- optimize(f=function(x){fopt(x)}, interval=c(0,10))
lambda.opt <- res$minimum
c.opt <- res$objective

#-----------------------------------------------------------------------------------

## c óptimo =  1.150589
##  lambda óptimo = 0.7932166

Asi el algoritmo que proponemos es el siguiente.

rfR <- function() {
  # Simulación por aceptación-rechazo
  # f a partir de la exponencial
  while (TRUE) {
    U <- runif(1)
    X <- rexp(1,lambda.opt)
    ngen <<- ngen+1 
    if (c.opt * U * dexp(X, lambda.opt) <= df(X)) return(X)
  }
}


rfRn <- function(n=1000) {
  # Simulación n valores
  x <- numeric(n)
  for(i in 1:n) x[i]<-rfR()
  return(x)
}
#-----------------------------------------------------------------------------------

Para ver el funcionamiento y eficiencia del algoritmo generaremos una muestra de tamaño \(nsim=10^{4}\), ademas de estudiar el número medio de generaciones.

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

##    user  system elapsed 
##    0.22    0.00    0.21

hist(x, breaks="FD", freq=FALSE)
curve(df(x),col='blue',lwd=2, add=TRUE)

#Tenemos que x sigue una distribucion F
#-----------------------------------------------------------------------------------

## Nº de generaciones =  11590
## Nº medio de generaciones =  1.159
## Proporción de rechazos =  0.1371872

Asi pues el algoritmo es bastante eficiente.

Para terminar el análisis del algoritmo realziaremos un ks.test.

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.0075984, p-value = 0.6106
## alternative hypothesis: two-sided

Con lo cual tenemos que x sigue una distribución dada por F.

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.

Comenzamos determinando \(F\), la funcion de distribución, tal que:

\[\begin{align*} F(x) &= \int_{0}^{x} f(t)dt\\ &= \int_{0}^{x} \frac{t^{3} -12t+20}{48}dt\\ &= \frac{1}{48} \left.\left(\frac{t^{4}}{4} -6t^{2} + 20t \right)\right|^{x}_{0}\\ &=\frac{1}{48}\left(\frac{x^{4}}{4} -6x^{2} + 20x \right) \end{align*}\]

Asi pues definimos:

\[\begin{align*} F(x)= \left\{ \begin{array}{lcc} 0 & \text{si } x<0 \\ \\ \frac{1}{48}\left(\frac{x^{4}}{4} -6x^{2} + 20x \right) & \text{si} \in [0,4] \\ \\ 1 & \text{si} x>4 \end{array} \right.. \end{align*}\]

#------------------------------------------------------------------------------------
#====================================================================================
###Usaremos una distribución U(0,4) 
#====================================================================================
#------------------------------------------------------------------------------------
# Obtención del valor c. 
c <- optimize(f=function(x){dej6(x)/dunif(x,0,4)}, maximum=TRUE, interval=c(0,4))
c.opt <- c$objective

## c óptimo =  2.999819

Asi pues el algoritmo que proponemos es el siguiente.

rej6R <- function() {
  # Simulación por aceptación-rechazo
  # F a partir de la U(0,4)
  while (TRUE) {
    U <- runif(1)
    X <- runif(1,0,4)
    ngen <<- ngen+1 
    if (c.opt * U * dunif(X,0,4) <= dej6(X)) return(X)
  }
}

rej6Rn <- function(n=1000) {
  # Simulación n valores
  x <- numeric(n)
  for(i in 1:n) x[i]<-rej6R()
  return(x)
}

Para estudiar el funcionamiento y eficiencia del algoritmo comezaremos generando una muestra de tamaño \(nsim=10^{4}\), ademas de estudiar el número medio de generaciones.

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

##    user  system elapsed 
##    0.61    0.00    0.64

#Análisis del Histograma
hist(x, breaks="FD", freq=FALSE)
curve(dej6(x),col='blue', add=TRUE)

#x sigue una distribución F

## Nº de generaciones =  29806
## Nº medio de generaciones =  2.9806
## Proporción de rechazos =  0.6644971

Asi pues, el algoritmo no es muy eficiente pero es el optimo.

#------------------------------------------------------------------------------------
# Prueba de Bondad de Agujste ks.test
#------------------------------------------------------------------------------------
ks.test(x,pej6(x),min=0,max=4)

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and pej6(x)
## D = 0.8856, p-value = 0.413
## alternative hypothesis: two-sided

Con lo cual aceptamos que x sigue una distribucion con funcion de densidad \(f\).

Ejercicio 7.

El tiempo (en años) entre dos inspecciones fiscales a una empresa es una variable aleatoria con densidad dada por:

\[\begin{align*} f(x)= \left\{ \begin{array}{lcc} 0.11-0.0018x-0.00003x^{2} & \text{Si} x \in [0,10] \\ \\ 0 & \text{Caso contrario} \\ \end{array} \right.. \end{align*}\]

Dar un algoritmo, lo más sencillo posible, que permita simular valores de esa variable. Comentrar el grado de eficiencia de dicho algoritmo frente a otras alternativas posibles para conseguir el mismo fin.

###Usaremos una distribución U(0,10) 
# Obtención del valor c. 
c <- optimize(f=function(x){dej7(x)/dunif(x,0,10)}, maximum=TRUE, interval=c(0,10))
c.opt <- c$objective

## c óptimo = 1.099999

Asi pues el algoritmo que proponemos es el siguiente:

rej7R <- function() {
  # Simulación por aceptación-rechazo
  # f a partir de una U(0,10)
  while (TRUE) {
    U <- runif(1)
    X <- runif(1,0,10)
    ngen <<- ngen+1 
    if (c.opt * U * dunif(X,0,10) <= dej7(X)) return(X)
  }
}


rej7Rn <- function(n=1000) {
  # Simulación n valores N(0,1)
  x <- numeric(n)
  for(i in 1:n) x[i]<-rej7R()
  return(x)
}

Para estudiar el funcionamiento del algoritmo propuesto generaremos una muestra de tamaÃ±os \(nsim=10^{4}\), ademas estudiaremos el número medio de generaciones.

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

##    user  system elapsed 
##    0.25    0.00    0.25

#Análisis del Histograma de x.
hist(x, breaks="FD", freq=FALSE)
curve(dej7(x),col='blue',lwd=2, add=TRUE)

## Nº de generaciones =  10930
## Nº medio de generaciones =  1.093
## Proporción de rechazos =  0.08508692

Comparacion con otros metodos.

Dado que:

## Nº medio de generaciones =  1.093

Como observamos que es aceptablemente cercano a 1, entonces el algoritmo es bastante eficiente. Dado que sigue una distribucion \(U(0,10)\), notemos que el

## c óptimo =  1.099999

Muy cercano a uno podemos usar los metodos de generacion de numeros aleatorios modificandolos en donde sea necesario para generar valores que sigan una distribucion \(U(0,10)\).

Además por el mismo analisis de que el valor óptimo de \(c \approx 1\), entonces la solución del problema de encontrar la mejor función auxiliar \(g\), es la trivial \(f=g\), con lo cual asumimos que \(x\) sigue una distribución \(U(0,10)\) en casi todas partes.

Por último para terminar el análisis del funcionamiento del algoritmo realizaremos un ks.test. Para ello comenzaremos determinando la función de distribución \(F\), a partir de la función de densidad \(f\).

\[\begin{align*} F(x) &= \int_{0}^{x} f(t)dt\\ &= \int_{0}^{x} 0.11-0.0018t-0.00003t^{2} dt\\ &= 0.11 x - \frac{0.0018}{2}x^{2} - \frac{0.00003}{3}x^{3} \end{align*}\]

Asi pues, definimos:

\[\begin{align*} F(x)= \left\{ \begin{array}{lcc} 0 & \text{si } x<0 \\ \\ 0.11 x - \frac{0.0018}{2}x^{2} - \frac{0.00003}{3}x^{3} & \text{si } \in [0,10] \\ \\ 1 & \text{si } x>10 \end{array} \right.. \end{align*}\]

#Kolmogorov-Smirnov Tests
ks.test(x,'pej7')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.0092969, p-value = 0.3531
## alternative hypothesis: two-sided

Asi pues x sigue una distribución F.

Utilizando un generador congruencial con parámetros \(a=5\), \(c=33\) y \(m=1024\) y semilla \(x_{0} =27\), simular tres tiempos entre inspecciones, por medio del algoritmo encontrado.

rej7R <- function() {
  # Simulación por aceptación-rechazo
  # f a partir de una U(0,10)
  while (TRUE) {
    U <- RANDC() #Generamos con el método congruencial.
    X <- runif(1,0,10)
    ngen <<- ngen+1 
    if (c.opt * U * dunif(X,0,10) <= dej7(X)) return(X)
  }
}

#Generacion de los datos. 
initRANDC(27) 
system.time(u <- rej7Rn(3))

##    user  system elapsed 
##    0.04    0.00    0.03

## 3.078943 5.578786 0.2110733

Describir un método para simular el número de inspecciones realizadas desde una dada hasta dentro de veinte años.

Obtener una muestra de tamaño 20, por cualquiera de los métodos antes propuestos.
Estudiar la convergencia de la media muestral. Mediante:

v <- rej7Rn(20)
plot(1:20, cumsum(v)/(1:20),col='blue', type="l", ylab="Media muestral", xlab="Nº de simulaciones")

#La convergencia de la media muestra nos aproximara al número de inspecciones realizadas desde una dada hasta dentro de veinte años.

Ejercicio 8.

El tiempo, en milisegundos, de anticipación (o retraso) de una señal respecto de otra en un sistema de comunicaciones puede considerar una variable aleatoria con función de distribución dada por \(F(x)=\frac{e^{x}}{1+e^{x}}\).

Dar un algoritmo, lo más detallado y simple posible, para simular valores de dicho tiempo de respuesta. Comparar la eficiencia computacional del método con la del algoritmo clásico para simular la distribución exponencial.

Proponemos el siguiente algoritmo basado en el método de aceptación-rechazo:

#================================================================
#Usando como distribución auxiliar N(0,sd)
#================================================================
#----------------------------------------------------------------
# Optimizamos el parámetro sd talque: 
#----------------------------------------------------------------

# Obtención de valores c y sd óptimos aproximados
fopt <- function(sd) {
  # Obtiene c fijado u
  optimize(f = function(x){dej8(x)/dnorm(x,0,sd)},
           maximum=TRUE, interval=c(-10,10))$objective
}

# Encontar sd que minimiza
res <- optimize(f=function(x){fopt(x)}, interval=c(-10,10))
sd.opt <- res$minimum
c.opt <- res$objective

{
  cat("c óptimo = ", c.opt)
  cat("\n sd.optimo = ", sd.opt)
}

## c óptimo =  1.262125
##  sd.optimo =  2.014061

#----------------------------------------------------------------
rej8R <- function() {
  # Simulación por aceptación-rechazo
  # F a partir de la normal(0,sd)
  while (TRUE) {
    U <- runif(1)
    X <- rnorm(1,0,sd.opt)
    ngen <<- ngen+1 
    if (c.opt * U * dnorm(X,0, sd.opt) <= dej8(X)) return(X)
  }
}


rej8Rn <- function(n=1000) {
  # Simulación n valores N(0,1)
  x <- numeric(n)
  for(i in 1:n) x[i]<-rej8R()
  return(x)
}

#GENERACION DE UNA  MUESTRA DE TAMAÑO  nsim=10^4 
system.time(x <- rej8Rn(nsim))

##    user  system elapsed 
##    0.19    0.00    0.19

#Analisis del Histograma
hist(x, breaks="FD", freq=FALSE)
curve(dej8(x),col='red', lwd=2, add=TRUE)

# x sigue una distribución F

## Nº de generaciones =  12709
## Nº medio de generaciones =  1.2709
## Proporción de rechazos =  0.213156

La eficiencia de el algoritmo es bastante buena.

Método de Inversión

Tenemos el sigueinte algoritmo basados en el método de inversión.

#================================================================
# Método de Inversión
#================================================================

#----------------------------------------------------------------
rdej8 <- function(){
  U <- runif(1)
  return(log(U)-log(1-U))
  
}

rdej8n <- function(n = 1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rdej8()
  return(x)
}

#GENERACIóN DE UNA  MUESTRA DE TAMAÑO  nsim=10^4
set.seed(54321)
system.time(u <- rdej8n(10^4))

##    user  system elapsed 
##    0.06    0.00    0.06

#Análisis del Histograma de u
hist(u, breaks = "FD", freq = FALSE)
curve(dej8(x),col='green', lwd=2, add = TRUE)

# u sigue una distribucion F.

Comparando los tiempos de computo, el metodo de Inversion es mas eficiente al método de aceptación usando la distribución auxiliar normal N(0,sd). Para finalizar el análisis de estos algoritmos realizaremos la prueba de Kolmogorov-Smirnov. Asi pues:

# Para la muestra generada con el método de aceptación-rechazo
ks.test(x,'pej8')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.010483, p-value = 0.2218
## alternative hypothesis: two-sided

# x sigue una distribucion F

#----------------------------------------------------------------
#Para la muestra generadad con el método de inversión.
ks.test(u,'pej8')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  u
## D = 0.0085815, p-value = 0.453
## alternative hypothesis: two-sided

# u sigue una distribución F

Ejercicio 9.

El número de horas diarias durante las cuales un terminal (de tipo A) de un laboratorio está siendo utilizado es una variable con función de distribución:

\[\begin{align*} F(x) &= \left\{\begin{array}{lcc} \frac{e^{x-2}}{1+e^{x-2}} & \text{si } x \geq 0\\ 0 & \text{caso contrario} \end{array}\right. \end{align*}\]

Se pide:

(a) Dar un algoritmo para simular dicha variable.

Proponemos el siguiente algoritmo basado en el método de inversión.

Propondremos un algoritmo basados en el método de aceptación-rechazo usando como función auxiliar una exponencia con parámetro lambda.

#================================================================
#Usando como distribucion auxiliar exponencial(lambda)
#================================================================
#----------------------------------------------------------------
# Optimizamos el parámetro lambda talque: 
#----------------------------------------------------------------

# Obtencion de valores c y sd óptimos aproximados
fopt <- function(lambda) {
  # Obtiene c fijado lambda
  optimize(f = function(x){dej9(x)/dexp(x,lambda)},
           maximum=TRUE, interval=c(0,10))$objective
}

# Encontar lambda que minimiza
res <- optimize(f=function(x){fopt(x)}, interval=c(0,10))
lambda.opt <- res$minimum
c.opt <- res$objective

#----------------------------------------------------------------
rnormAR <- function() {
  # Simulación por aceptación-rechazo
  # Normal estandar a partir de la normal
  while (TRUE) {
    U <- runif(1)
    X <- rexp(1,lambda.opt)
    ngen <<- ngen+1 
    if (c.opt * U * dexp(X,lambda.opt) <= dej9(X)) return(X)
  }
}


rnormARn <- function(n=1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rnormAR()
  return(x)
}

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

##    user  system elapsed 
##    0.36    0.00    0.38

hist(x, breaks="FD", freq=FALSE)
curve(dej9(x),col='orange',lwd=2, add=TRUE)

## Nº de generaciones =  18519
## Nº medio de generaciones =  1.8519
## Proporción de rechazos =  0.460014

#Ks.test
ks.test(x,'pej9')

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.11982, p-value < 2.2e-16
## alternative hypothesis: two-sided

Rechazamos que x siga una distribución F.

Asi pues, si analizamos F, debemos notar que:

Existe un punto de discontinuidad en x=0, ya que:

\[\begin{align*} \lim_{x \to 0^{-}} F(x) &\neq \lim_{x \to 0^{+}} F(x)\\ 0 &\neq \frac{e^{-2}}{1+ e^{-2}} \end{align*}\]

Con lo cual ese punto de discontinudad es lo que origina problemas, asi pues nos quedaremos con estas aproximaciones a los valores.

(b) Si para otro tipo de terminales (tipo B) el tiempo de utilización tiene distribución \(exp \left( \frac{1}{2}\right)\), detallar un algoritmo para aproximar, por simulación, la probabilidad de que un terminal de tipo B sea más utilizado que uno de tipo A.

Generar n muestras de tamaño m que sigan una distribución \(F\).

1.1. Generar n muestras de tamaño m que sigan una distribución \(exp(0.5)\)

Realizar un estudio de contraste, usando \(ks.test(x,v,alternative='less')\)
Estudiar la convergencia de la media muestral de los p-valores originados en el paso anterior.

Asi pues, tenemos:

set.seed(1727)
n <- 50
nsim <- 100
estadistico <- numeric(nsim)
pvalor <- numeric(nsim)

# Realizar contrastes
for(isim in 1:nsim) {
  u <- rnormARn(nsim) 
  v <- rexp(nsim,0.5)  # Generar
  tmp <- ks.test(u,v,alternative = 'less',min=0,max=10)
  estadistico[isim] <- tmp$statistic
  pvalor[isim] <- tmp$p.value
}
plot(1:nsim, cumsum(pvalor)/(1:nsim),col='red', type="l", ylab="Media muestral", xlab="Nº de simulaciones")

Ejercicio 10.

Escribir el código necesario para generar, por el método de inversión, una muestra de n observaciones de una distribución de Cauchy.

Proponemos, el siguiente algoritmo basado en el método de inversión.

#===================================================================================
# Distribucion de Cauchy - Método de Inversión. 
#===================================================================================
#-----------------------------------------------------------------------------------
rcauchy <- function(){
  U <- runif(1)
  return(tan(pi*(U-0.5)))
}

rcauchyn <- function(n = 1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rcauchy()
  return(x)
}

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

set.seed(1727)
nsim <- 10^4
system.time(x <- rcauchyn(nsim))

##    user  system elapsed 
##    0.06    0.00    0.06

#Al ser un algoritmo basado en el método de inversión, tenemos que es eficiente.

(b) Representar el histograma (limitar el rango, es decir, \(xlim=c(-10,10)\) y compararlo con la densidad teórica \((dcauchy)\).

#Análisis del Histograma.
hist(x, breaks = "FD", xlim = c(-10,10), freq = FALSE)
curve(dcauchy(x),col='blue', lwd=2, add = TRUE)

# A primera vista x sigue una distribución de Cauchy.

(c) Obtener conclusiones sobre la existencia de una media teórica a partir de la media muestral aproximada por simulación (estudiar la convergencia de la media muestral). Suponiendo que el vector x contiene las simulaciones, estudiar la convergencia de la media muestral mediante el gráfico:plot(1:nsim, cumsum(x)/(1:nsim), type=“l”, ylab=“Media muestral”, xlab=“Nº de simulaciones”)

plot(1:nsim, cumsum(x)/(1:nsim),col='orange',lwd=2, type="l", ylab="Media muestral", 
     xlab="Nº de simulaciones")

#Aunque no esta definida la media para una distribución de Cauchy, al aumentar el número de simulaciones la suceción parece converger a 0. Lo cual indicaria la existencia de una media muestral aproximada.

Ejercicio 11.

El tiempo de respuesta (en centésimas de segundo) de un servidor de bases de datos es una variable con función de densidad:

\[\begin{align*} f(x) &= x e^{-x}, \text{ si } x \geq 0 \end{align*}\]

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:

\[\begin{align*} g(x) &= \lambda e^{-\lambda x}, \text{ si } x \geq 0 \end{align*}\]

(a) Aproximar numericamente el parámetro optimo \(\lambda_{opt}<1\) y la cota óptima \(c_{opt}\) de la densidad auxiliar y compararlos con los valores teoricos: \(\lambda_{opt} = \frac{1}{2}\) y \(c_{opt}=\frac{4}{e}\)

#===================================================================================
## Usaremos la distribución exponencial, con parámetro lamda
#===================================================================================
#-----------------------------------------------------------------------------------
# Obtención de valores c y lambda óptimos aproximados
fopt <- function(lambda) {
  # Obtiene c fijado lambda
  optimize(f = function(x){dej11(x)/dexp(x,lambda)},
           maximum=TRUE, interval=c(0,10))$objective
}
# Encontar lambda que minimiza
res <- optimize(f=function(x){fopt(x)}, interval=c(0,10))
lambda.opt <- res$minimum
c.opt <- res$objective

## Lambda óptimo =  0.5000085
##  c óptimo =  1.471518

curve(c.opt*dexp(x,lambda.opt), 0, 10, col = 3, lwd = 3)
curve(dej11(x),col='blue',add=TRUE)

#-----------------------------------------------------------------------------------
#lambda.opt
# [1] 0.5000085 esta bien aproximado al teorico 1/2
#c.opt
#[1] 1.471518 esta bien aproximado al teorico de 4/exp(1)
#-----------------------------------------------------------------------------------

Proponemos el siguiente algoritmo:

rnormAR <- function() {
  # Simulación por aceptación-rechazo
  while (TRUE) {
    U <- runif(1)
    X <- rexp(1,lambda.opt)
    ngen <<- ngen+1 
    if (c.opt * U * dexp(X,lambda.opt) <= dej11(X)) return(X)
  }
}

rnormARn <- function(n=1000) {
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rnormAR()
  return(x)
}

#-----------------------------------------------------------------------------------

(b) Generar una muestra de 1000 observaciones de la distribución de interés (tomando como semilla inigial el número 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.

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

##    user  system elapsed 
##    0.03    0.00    0.03

#-----------------------------------------------------------------------------------

## Nº de generaciones =  1474
## Nº medio de generaciones =  1.474
## Proporción de rechazos =  0.3215739

Como vimos en el ejercicio 4 una mejor función auxiliar para este ejercicio es gamma(2,lambda), donde lambda era muy cercano a 1. Aun asi, como se explico antes gamma no es una distribución sencilla de simular, por lo cual, el uso de exp(0.5) podria ser mejor cosdierando la facilidad de implementación.

(c) Representar el histograma y compararlo con la densidad teorica.

hist(x, breaks="FD", freq=FALSE)
curve(dej11(x), col='red',lwd=2,add=TRUE)

#Tenemos que x sigue una distribución con densidad f.

Métodos especificos para la simulación de distribuciones notables

Distribuciones Continuas.

Las distribuciones Gamma y Erlang

#==============================================================
# Distribucion de Erlang
#==============================================================
rErlang <- function(a=1,p=1){
  #Por defecto a=1, p=1.
  U<-numeric(p)
  X<-numeric(p)
  for (i in 1:p) {
    U[i] <- runif(1)
    X[i] <- - log(U[i])/a 
  }
  return(Y<-sum(X))
}

rErlangn <- function(n=1000){
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rErlang()
  return(x)
}

u <- rErlangn()
hist(u,freq=FALSE)
curve(dgamma(x,1,1),col='red',lwd=2,add=TRUE)

#==============================================================
# Distribucion de Erlang Optimizado
#==============================================================
rErlop <- function(a=1,p=1){
  p <- ceiling(p)
  L <- -1/a
  S <- 1
  for (i in 1:p){
    U <- runif(1)
    S <- S*U
  }
  return(Y<-L*log(S))
}

rErlopn <- function(n=1000){
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-rErlop()
  return(x)
}

u <- rErlopn()
hist(u,freq=FALSE)
curve(dgamma(x,1,1),col='blue',lwd=2,add=TRUE)

#==============================================================
# Funcion Auxiliar doble exponencial modificada.
#==============================================================

rdexp <- function(p = 1,lambda = 1){
  # Simulación por inversión
  # Doble exponencial
  U <- runif(1)
  if (U<(exp(p-1)/2)) {
    return(log(2*U)/lambda)
  } else {
    return(-log(2*(1-U))/lambda)
  }
}
#==============================================================

#==============================================================
# Distribucion Gamma(1,p)
#==============================================================
#--------------------------------------------------------------
Tx <- function(x,p,thetha){
  
  (abs(((thetha-1)*x)/(thetha*(p-1)))^(p-1))* exp(-x +  ((abs(x-p+1)+((p-1)*(thetha+1)))/(thetha))   ) 
}


#--------------------------------------------------------------
gamR <- function(p){
  #============================================================
  lambda.opt <- 2/(1+sqrt((4*p) - 3))
  thetha <- 1/lambda.opt
  #============================================================
  X<-rdexp(p,lambda.opt)
  while (X<0) {
    X<-rdexp(p,lambda.opt)
  }
  U<-runif(1)
  ifelse((U<=Tx(X,p,thetha)), return(X),return(gamR(p)))
  
}
#--------------------------------------------------------------
gamRn <- function(n=1000,p){
  # Simulación n valores 
  x <- numeric(n)
  for(i in 1:n) x[i]<-gamR(p)
  return(x)
}

Generacion de Variables Aleatorias Continuas

Amagua Sandobalin Gabriel Sebastian, Diaz Quichimbo Daniel Arturo

Mayo del 2019