Reduccion_de_Varianza

Diseño de experimentos de simulación.

Ejercicios Propuestos.

Ejercicio 1.

Considérese el algoritmo encontrado en el ejercicio 3, propuesto en el tema 3. Si posteriormente se quiere proceder de forma que los impactos cercanos y leganos al centro estén equilibrados. ¿cómo se podrían simular ahora las distancias al centro de 10 impactos?.

Recordando el ejercicio 3, del tema 3, tenemos:

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.07    0.00    0.08

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

Comenzaremos planteandonos los siguientes estratos:

Valores Menores de la distribucion ($50\%$)
Valores Mayores de la distribucion ($50\%$)

Asi pues consideramos el siguiente psudocódigo:

Generar $U_{i} \sim U[0,1]$, para $i = 1,2,...10$

Sí $i \leq 5$, entonces $U_{i} = 0.5*U_{i} + 0.5$

Sí $5 < i \leq 10$, entonces $U_{i} = 0.5*U_{i}$

Para $i=1,..,10$, devolver: \[\begin{equation*} X_{i} = U_{i} ^{\frac{1}{3}} \end{equation*}\]

nsim <- 10 
U <- runif(10)

for (isim in 1:10) {
  if(isim <=5){U[isim]=0.5*U[isim] + 0.5}
  else{U[isim]=0.5*U[isim]}
}

x <- U**(1/3)
hist(x,freq = FALSE)

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

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

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

Usando el test de kolmogorov-smirnov tenemos que la muestra generada satisface la distribución además de cumplir con la condición requerida.

Ejercicio 2.

Una nueva componente llega cada 40 segundos a una cadena de ensamblado. El tiempo necesario para ensamblar dicha componente a la pieza matriz se supone una variable aleatoria con media de 30 segundos (resulta muy importante haber estimado su distribución: $N(30,1)$, $N(30,5)$, $\Gamma(2,60)$, etc…). Sí un tiempo de ensamblado es superior a 40 segundos, las comonentes que van llegando se acumulan hasta que se les vaya dando salida. Dar un algoritmo que permita responder a las siguientes preguntas ¿Cuál es la probabilidad de que una componente que llega tenga que esperar?, ¿Cuál es la probabilidad de que estén más de tres piezas esperando?

Consideremos el siguiente psudocódigo.

Debemos notar que deberiamos tener una estimación de la distribución que sigue la variable tiempo de ensamblaje.

Una vez tenemos la muestra (generada con cualquiera de las distribuciones), tenemos:

Calcular la probabilidad de tener un tiempo de ensamblaje mayor a 40 segundos tal que: \[ \begin{equation*} p = \frac{\text{Número de exitos}}{\text{Tamaño de la muestra}} \end{equation*} \] Ahora bien, tenemos que la variable aleatoria tiempo de ensamblaje mayor a 40 sigue una distribucion $Bi(n,p)$, donde $n$ es el tamaño de la muestra.
Calculamos la probabilidad de obtener 3 o mas tiempos de espera mayor a 40, para ello podemos calcular mas facilmente $1$-la probabilidad de obtener menos de 3 tiempos de espera mayor a 40. \[ \begin{equation*} \sum_{i=0}^{2} \left( \begin{array} &n\\ i\end{array} \right) p^{n} (1-p)^{i} \end{equation*} \]

Ejercicio 3.

El número de toneladas de pan producidad diariamente por una empresa panificadora tiene una distribución de Pareto con párametros 3 y 2. Suponiendo que la demanda diaria de pan (en toneladas) tiene distribución $N(3,0.5)$ y que es independiente de la producción, calcular la probabilidad de que un día concreto la demanda no sea satisfecha. Suponiendo una pérdida de 0.05 euros por cada kilo de pan no vendido y una penalización de 0.2 euros por cada kilo de pan demandado por clientes y no entregado (por falta de existencias), dar un algoritmo para calcular las pérdidas medias diarias en concepto de excendentes o demandas insatisfechas.

Comenzaremos definiendo el siguiente estadistico. \[ \begin{equation*} A = X - Y \end{equation*} \]

donde $X \sim Pareto(3,2)$ y $Y \sim N(3,0.5)$

library(rmutil )

## Warning: package 'rmutil' was built under R version 3.5.2

## 
## Attaching package: 'rmutil'

## The following object is masked from 'package:stats':
## 
##     nobs

## The following objects are masked from 'package:base':
## 
##     as.data.frame, units

Realizaremos un experimento de simulación para estimar la probabilidad, que consistira en realizar 100 simulaciones con tamaños de muestra 50.

n <- 50
nsim <- 100
probabilidad <- numeric(nsim)

for (isim in 1:nsim) {
  s <<- 0
  x <- rpareto(n,3,2) #Produccion
  y <- rnorm(n,3,0.5) #Demanda
  A <- x-y
  for (jsim in 1:n) {
    if(A[jsim]<0){s <<- s+1}
    
  }
  probabilidad[isim] <- s/n
  
}
plot(1:nsim, cumsum(probabilidad)/(1:nsim),col='red', type="l", ylab="Media muestral", xlab="Nº de simulaciones")

mean(probabilidad)

## [1] 0.7338

Para calcular las perdidas diarias realizaremos nuevamente un experimento de simulación con 100 simulaciones con tamaño de muestra 50.

n <- 50
nsim <- 100
perdida_demanda <- numeric(nsim)
perdida_excedentes <- numeric(nsim)

for (isim in 1:nsim) {
  s <<- 0
  k <<- 0
  x <- rpareto(n,3,2) #Produccion
  y <- rnorm(n,3,0.5) #Demanda
  A <- x-y
  for (jsim in 1:n) {
    ifelse(A[jsim]<0,s <<- s-(A[jsim]*0.2), k<<- k+(A[jsim]*0.05) )
    
  }
  perdida_demanda[isim] <- s
  perdida_excedentes[isim] <- k
  
}
plot(1:nsim, cumsum(perdida_demanda)/(1:nsim),col='orange', type="l", ylab="Media muestral", xlab="Nº de simulaciones")

{
 cat("Perdida media por demandas insatisfechas =", mean(perdida_demanda))
}

## Perdida media por demandas insatisfechas = 15.08178

plot(1:nsim, cumsum(perdida_excedentes)/(1:nsim),col='blue', type="l", ylab="Media muestral", xlab="Nº de simulaciones")

{
 cat("Perdida media por excedentes =", mean(perdida_excedentes))
}

## Perdida media por excedentes = 3.891495

Ejercicio 4.

En una cadena productiva, el tiempo que tarda una componente en llegar a una máquina ensambladora desde que llegó la pieza anterior sigue una distribución exponencial de media 10 segundos. El tiempo (en segundos) que emplea la máquina ensambladora con cada componente se ajusta a una distribución $\Gamma(3,20)$. Cada vez que una nueva componente llega a la fase de ensamblaje y la máquina está ocupada, dicha pieza se desvía a otra linea de producción distinta. Dar un algoritmo, lo más oreciso posible, para poder aproximar mediante simulación el porcentaje de tiempo que la máquina ensambladora mediante simulación el porcentaje de tiempo que la máquina ensambladora perderá esperando la llegada de una nueva pieza. Justificar el funcionamiento de dicho algoritmo.

Proponemos el siguiente pseudo-código:

Generar una muestra X de tamaño $n$ tal que $X \sim exp(10)$

** Notemos que si la máquina ensambladora se encuentra trabajando, la nueva componente se desvía a otra línea de producción** Así pues, no necesariamente (por optimización computacional), se tendrá una muestra Y de tamaño $n$ tal que $Y \sim \Gamma(3,20)$, por facilidad computacional, asumiremos que si se necesitará una muestra de $Y$ de tamaño $n$.

Generar una muestra Y de tamaño $n$ tal que $Y \sim \Gamma(3,20)$
Inicializar t_espera$=0$, $i=1$, $j=1$, $Z=Y$.

Tomamos una muestra auxiliar Z=Y para facilitar las operaciones computacionales.
Mientras $i < n $ y $j<n$:

4.1. Sí $X_{i} \geq Y_{j}$

  Hacer t_resta = t_resta+x[i]-y[j].
  Hacer x[i+1]=x[i+1]-x[i]+y[j].
  Hacer i=i+1
  Hacer j=j+1

4.2. Caso contrario.

  Hacer y[j]=y[j]-x[i]
  Hacer i=i+1

Devolver \[ \begin{equation*} p = \frac{\text{t_espera}}{\text{t_espera}+\sum_{k=1}^{j}Y_{j}} *100 \end{equation*} \]

Realizaremos un experiemnto para probar el funcionamiento del algoritmo usando un tamaño de muestra igual a 50.

n <- 50
x <- rexp(n,10)
y <- rgamma(n,3,20)
z <- y
t_espera <- 0 
i <- 1
j <- 1


while(i < n && j <n){
  if(x[i]>=y[j]){
    t_espera <- t_espera +x[i]-y[j]
    x[i+1] <- x[i+1]-(x[i]-y[j])
    i <- i+1
    j <- j+1
  }
  else{
    y[j]<-y[j]-x[i]
    i<-i+1
  }
  
}

p <- (t_espera/(t_espera+sum(z[1:j])))*100

{
  cat("Porcentaje de tiempo esperando =", p)
}

## Porcentaje de tiempo esperando = 34.6114

Ejercicio 5.

Una máquina produce tiras de goma de longitud aleatoria con distribución exp(2) (en metros). Despues de fabricarlas, las tiras de goma se pasan a otra máquina que las estira hasta que rompen en ds (para comprobar su elasticidad). Suponiendo que el lugar por donde se rompe cada tira es aleatorio y con distribución uniforme a lo largo de su longitud, descrbir detalladamente un algoritmo que simule las longidutes de los dos trozos en los que se rompe cada goma.