Ejercicios 7.3 (Ruben)

1.- 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.

#Usando el método de inversión
rcauchy <- function(){
  U <- runif(1)
  return(tan(pi*(U-0.5)))
}

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

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

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

##    user  system elapsed 
##    0.16    0.04    0.28

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

hist(x, breaks = "FD", xlim = c(-10,10), freq = FALSE)
curve(dcauchy(x),col='red', lwd=2, add = TRUE)

Notemos que parece que sigue una distribucion de Cauchy.

3. 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:10^4, cumsum(x)/(1:10^4),lwd=2, type="l", ylab="Media muestral", xlab="Nº de simulaciones")

Notemos que mientras mas simulaciones la variable converge a 0.

2.- 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) = x e^{-x}, \text{ 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) = \lambda e^{-\lambda x}, \text{ si } x \geq 0 \] 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}\)

#Funcion de densidad
fun_ejer_2 <- function(x){
  ifelse((x>=0),x*exp(-x),0)
}

#Valores c y lambda óptimos
fopt <- function(lambda) {
  # Fijamos lambda
  optimize(f = function(x){fun_ejer_2(x)/dexp(x,lambda)},
           maximum=TRUE, interval=c(0,10))$objective
}
# Encontar lambda
res <- optimize(f=function(x){fopt(x)}, interval=c(0,10))
lambda.opt <- res$minimum
c.opt <- res$objective

cat("\n c óptimo = ", c.opt)

## 
##  c óptimo =  1.471518

#que es proximo a 4/exp(1)

cat("Lambda óptimo = ", lambda.opt)

## Lambda óptimo =  0.5000085

#que es proximo a 1/2

2. 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.

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

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

#fijamos la semilla
set.seed(500)
n<-0
system.time(x <- rfunn(1000))

##    user  system elapsed 
##    0.08    0.00    0.15

cat("Nº de generaciones = ", n)

## Nº de generaciones =  1471

cat("\nNº medio de generaciones = ", n/1000)

## 
## Nº medio de generaciones =  1.471

cat("\nProporción de rechazos = ", 1-n/1000, "\n")

## 
## Proporción de rechazos =  -0.471

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

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