Apellidos: Amagua Sandobalin

Nombre: Gabriel Sebastián

  1. Describir un algoritmo basado en el método de inversión que permita generar observaciones de una variable aleatoria con distribución exponencial “inflada con ceros” (zero inflated), con función de distribución: \[F(x)=\left\{ \begin{array}{ll} 0 & \text{si } x < 0,\\ \pi + \left( 1-\pi \right)\left(1-e^{-\lambda x}\right) & \text{si } x \ge 0. \end{array} \right.\]

Observación Hay que notar que la función de distribución F, dada es discontinua en el punto \(x=0\), además dicha discontinuidad es inevitable.

Así pues en tal caso, no podriamos aplicar el método de Inversión para generar observaciones de está variable aleatoria, aun así. Procederemos ignorando la observación anterior.

Comenzamos determinando la inversa.

Para \(x \geq 0\), tenemos:

\[\begin{align*} \pi + \left( 1-\pi \right)\left(1-e^{-\lambda x}\right) &= \mu\\ \left( 1-\pi \right)\left(1-e^{-\lambda x}\right) &= \mu - \pi \\ e^{-\lambda x} &= \frac{\mu-\pi}{1-\pi} \\ e^{-\lambda x} &= 1 - \frac{\mu-\pi}{1-\pi} \\ \ln(e^{-\lambda x}) &= \ln\left(1 - \frac{\mu-\pi}{1-\pi}\right)\\ -\lambda x &= \ln\left(1 - \frac{\mu-\pi}{1-\pi}\right)\\ x &= -\frac{1}{\lambda} \ln\left(1 - \frac{\mu-\pi}{1-\pi}\right)\\ \end{align*}\]

Observación Debemos notar que existen problemas en el argumento del logaritmo de la función inversa definida anteriormente. Se debe tener que:

\[\begin{align*} 0 &< 1 - \frac{\mu-\pi}{1-\pi} \\ -1 &< - \frac{\mu-\pi}{1-\pi} \\ 1 &> \frac{\mu-\pi}{1-\pi} \\ 1-\pi &< \mu-\pi \\ 1 &< \mu \end{align*}\]

Observación Además como \(x \geq 0\), tenemos:

\[\begin{align*} 0 &\leq -\frac{1}{\lambda} \ln\left(1 - \frac{\mu-\pi}{1-\pi}\right) \\ 0 &\geq \ln\left(1 - \frac{\mu-\pi}{1-\pi}\right)\\ 1 &\geq 1 - \frac{\mu-\pi}{1-\pi}\\ 0 &\geq - \frac{\mu-\pi}{1-\pi}\\ 0 &\leq \frac{\mu-\pi}{1-\pi}\\ 0 &\geq \mu-\pi \\ \mu &\leq \pi \end{align*}\]

Con lo cual basados en las observaciones anteriores tenemos que; \(\mu \sim U(1,\pi)\)

Asi pues proponemos el siguiente algoritmo:

# Ejercicio 1
# Amagua Sandobalin Gabriel Sebastian.

pF <- function(x,lambda){
  ifelse((x<0),0,
         (pi+((1-pi)*(1-exp(-lambda*x)))))
}

fR <- function(lambda){
  #Generamos un numero aleatorio con distribucion
  # U(1,pi) pues Segun los calculos anexados en las hojas U debe 
  # ser maypor que 1 y menor o igual que pi
  U <- runif(1,1,pi)
  return(x  <- (-((1)/(lambda))*log(1-((U-pi)/(1-pi)))))
}

fRn <- function(lambda,n=10000) {
  x <- numeric(n)
  for(i in 1:n) x[i]<-fR(lambda)
  return(x)
}
set.seed(1727)
u <- fRn(5)
hist(u,freq = FALSE)

  1. Estamos interesados en una variable con función de densidad \(Gamma(s,r)\): \[f(x)=\frac{r^{s}}{\Gamma (s)}x^{s-1}e^{-rx}\text{ si }x\geq 0,\] (siguiendo la notación de la función dgamma(x, shape, rate) de R). Escribir el código necesario para generar, por el método de aceptación-rechazo, una muestra de \(n\) observaciones de una distribución \(Gamma(3,3)\) empleando como densidad auxiliar una exponencial (dexp(x, rate)): \[g(x)=\lambda e^{-\lambda x}\text{ si }x\geq 0.\] (NOTA: No emplear las funciones de densidad implementadas en R).
# Ejercicio 2
# Amagua Sandobalin Gabriel Sebastian.
# Funciones auxiliares para generar, por el método de
# aceptación-rechazo

#====================================================================================
# FUNCIÓN DE DENSIDAD PARA GAMMA(3,3) 
#====================================================================================
df <- function(x){
  (27/gamma(3))*(x^2)*(exp(-3*x))
  
}

#====================================================================================
# FUNCIÓN DE DISTRIBUCIÓN  PARA GAMMA(3,3) 
#====================================================================================

pF <- function(x){
  ifelse((x<0),0,(1-((exp(-3*x)*(9*x^2 + 6*x + 2)/(2)))))
}

#====================================================================================
# FUNCIÓN AUXILIAR DE DENSIDAD PARA EXPONENCIAL(lambda) 
#==================================================================================== 
daux <-function(x,lambda){
  lambda*exp(-lambda*x)
}
#====================================================================================
# FUNCIÓN PARA GENERAR UN VALOR ALEATORIO CON DISTRIBUCIÓN EXPONENCIAL(lambda) 
#==================================================================================== 

auxR <- function(lambda){
  U <- runif(1)
  return(x  <- (-log((1-U)/lambda)))
}
1.  Aproximar numéricamente el parámetro óptimo ($\lambda
    _{opt}$) y la cota óptima ($c_{opt}$) de la densidad auxiliar y
    compararlos con los valores teóricos.
#-----------------------------------------------------------------------------------
# Obtención de valores c y lambda óptimos aproximados
fopt <- function(lambda) {
  # Obtiene c fijado lambda
  optimize(f = function(x){df(x)/daux(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
#-----------------------------------------------------------------------------------
curve(c.opt*daux(x,lambda.opt), 0, 10, col = 3, lwd = 3)
curve(df(x),col='blue',add=TRUE)

#-----------------------------------------------------------------------------------
## Lambda obtenido numéricamente =  0.9999934 
##  C obtenido numéricamente = 1.827026
2.  Generar una muestra de 1000 observaciones de la distribución de
    interés tomando como semilla inicial los cuatro primeros dígitos
    de la CI. Obtener el tiempo de CPU que tarda en generar la
    secuencia y calcular el número medio de generaciones de la
    distribución auxiliar. Realice un contraste de hipótesis de 
    Bondad de Ajuste

Así pues proponemos el siguiente algoritmo:

#-----------------------------------------------------------------------------------
rnormAR <- function() {
  # Generación de un resultado aleatorio con dsitribución Gamma(3,3), usando
  # como función auxiliar exp(lambda.opt)
  while (TRUE) {
    U <- runif(1)
    X <- auxR(lambda.opt)
    ngen <<- ngen+1 
    if (c.opt * U * daux(X,lambda.opt) <= df(X)) return(X)
  }
}

rnormARn <- function(n=1000) {
  # Generación de n resultados aleatorios con distribución Gamma(3,3)
  x <- numeric(n)
  for(i in 1:n) x[i]<-rnormAR()
  return(x)
}
#-----------------------------------------------------------------------------------
set.seed(1727)
nsim <- 10^3
ngen <- 0
system.time(x <- rnormARn(nsim))
##    user  system elapsed 
##    0.04    0.00    0.05
# El tiempo de computo es pracitamente 0, asi pues el 
# algoritmo propuesto es bastante eficiente.
## 
## Nº de generaciones =  1894
## Nº medio de generaciones =  1.894
## Proporcion de rechazos =  0.4720169

El número medio de comparaciones es 1.894, lo cual no es tan cercano a 1 como para afirmar que la función auxiliar usada sea la mas eficiente, asun asi es óptima.

# Kolmogorov-Smirnov Tests
ks.test(x,'pF')
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.020647, p-value = 0.7875
## alternative hypothesis: two-sided
#Como el p-valor es aproximadamente 0.8, aceptamos que x sigue una 
#distribucion Gamma(3,3)
3.  Representar el histograma y compararlo con la densidad teórica.
hist(x, breaks="FD", freq=FALSE)
curve(df(x),col='red', add=TRUE)

# A simple vista el algoritmo propuesto funciona. 
# x sigue una distribución Gamma(3,3)

  1. Una compañia que vende ordenadores quiere determinar cuántas licencias de un determinado sistema operativo comprar a Microsoft para suministrarla con dichos ordenadores durante un año. Como el sistema operativo tiende a cambiar de un año para otro, la compañia querría comprar a Microsoft tantas licencias como ordenadores venda ese año, pero, evidentemente, esta cantidad le es desconocida de antemano.Una licencia del sistema operativo actual le cuesta a la compañia 75 euros y al suministrarla con el ordenador le genera unos ingresos brutos de 120 euros.

    Cuando Microsoft lanza al mercado un nuevo sistema operativo, la compañia puede devolver las licencias sobrantes, recibiendo 60 euros por cada una (de los 75 que pagá). Además, la compañia ha estimado que la cantidad de ordenadores que puede vender en el año se distribuye según la siguiente distribución:

    Demanda Prob.
    800 0.15
    850 0.20
    975 0.25
    1200 0.30
    1300 0.10

    La cuestión que se plantea la compañia es cuántas licencias del sistema operativo solicita a Microsoft.

    Se trata de un problema complicado, por lo que la empresa se plantea comparar las siguientes alternativas de pedido: 975, 1000, 1025, 1050, 1075 y 1100. Desarrolle un algoritmo para simular el beneficio y a partir de esto poder estimar el beneficio esperado de la compañia para cada uno de los tamaños de pedido contemplados.

    Nota: Utilice los cuatro primeros dígitos de su CI, y utilice n=1000 escenarios (años) de simulación.

Como este ejercicio es un caso particular del problema del vendedor de periodicos, proponemos el siguiente algoritmo.

#====================================================================================
# DATOS INICIALES
#====================================================================================
s <- c(800,850,975,1200,1300)
p <- c(0.15,0.20,0.25,0.30,0.10) 
#cp <- c(0.15,0.35,0.60,0.90,1)
#Ed <- sum(s*p)
#cv <- (120-75)/((120-75)+(75-60))
cu <- 120-75
co <- 75-60
ap <- c(975,1000,1025,1075,1100)
#====================================================================================
# GENERACIÓN POR TABLAS GUÍA. 
#====================================================================================
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)
  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)
}
#====================================================================================
# FIJAMOS LA SEMILLA Y EL TAMAÑO DE LA MUESTRA
#====================================================================================

set.seed(1727)
nsim <- 1000

#====================================================================================
#====================================================================================
# ALGORITMO PROPUESTO PARA DETERMINAR EL BENEFICIO PARA UNA ALTERNATIVA
#====================================================================================

benefit <- function(x,Q,nsim,cu){
  beneficio <- numeric(nsim)
  for (isim in 1:nsim) {
    if(x[isim]>=Q ){beneficio[isim]<-(Q*cu)}
    else{beneficio[isim]<-((x[isim])*cu)}
    
  }
  return(beneficio)
}
#====================================================================================
# ALGORITMO PROPUESTO PARA DETERMINAR LA PERDIDA PARA UNA ALTERNATIVA
#====================================================================================

loss <- function(x,Q,nsim,co){
  perdida <- numeric(nsim)
  for (isim in 1:nsim) {
    if(x[isim]>=Q ){perdida[isim]<-(0)}
    else{perdida[isim]<-((Q-x[isim])*co)}
    
  }
  return(perdida)
}


#=======================================================================================
#====================================================================================
# ALGORITMO PROPUESTO PARA DETERMINAR EL BENEFICIO, PERDIDA Y PROFIT
# PARA N ALTERNATIVAs
#====================================================================================


Beneficios_para_cada_alternativa <- function(y,nsim,cu,co){
  beneficio_por_alternativa <- numeric(length(y))
  perdida_por_alternativa <- numeric(length(y))
  profit_por_alternativa <- numeric(length(y))
  x <- rfmp.tabla(s,p,4,nsim)
  for (asim in 1:length(y)) {
    beneficio_por_alternativa[asim] <- (sum(benefit(x,y[asim],nsim,cu)))
    perdida_por_alternativa[asim] <- (sum(loss(x,y[asim],nsim,co)))
    profit_por_alternativa[asim] <- sum(benefit(x,y[asim],nsim,cu)-loss(x,y[asim],nsim,co))
    
  }
  
  return(data.frame("Alternativa"=y,"Beneficio"=beneficio_por_alternativa,"Perdida"=perdida_por_alternativa,"Profit"=profit_por_alternativa))
  
  
}
ejercicio_3 <- Beneficios_para_cada_alternativa(ap,nsim,cu,co)
##   Alternativa Beneficio Perdida   Profit
## 1         975  41560875  771375 40789500
## 2        1000  42005250  998250 41007000
## 3        1025  42449625 1225125 41224500
## 4        1075  43338375 1678875 41659500
## 5        1100  43782750 1905750 41877000