Unidad03-Ejercicios en clase

# Simulación de un proceso de Poisson
simula.pp=function(lambda,neventos=NULL,tmax=NULL){
  # se lanza la simulación hasta conseguir neventos o hasta un instante tmax
  # mensaje de error si no se introduce alguna regla de parada
  if(is.null(tmax) & is.null(neventos)){
    return(cat("Introduce regla de parada: neventos o tmax"))
  }
  # si no se introduce tmax, la regla de parada es neventos
  else if(is.null(tmax)){
    # simulación de los tiempos entre eventos sucesivos
    t=rexp(neventos,lambda)
  }
  # si no se introduce neventos, la regla de parada es tmax
  else if(is.null(neventos)){
    t=c()
    tt=0 # inicialización del reloj
    i=1
    t[1]=rexp(1,lambda)
    # si la primera llegada se produce después de tmax, contabilizamos 0 eventos
    if(t[1]>tmax)
      return(data.frame(eventos=0,t=0,ttotal=0))
    # en otro caso, continuamos hasta llegar a tmax
    else{
    while(tt<=tmax){
      i=i+1
      t[i]=rexp(1,lambda)
      tt=tt+t[i]
        }
    # se seleccionan sólo los eventos que se produjeron antes de tmax
    t=t[cumsum(t)<=tmax]
  }}
return(data.frame(eventos=1:length(t),t,ttotal=cumsum(t)))
}

# Función para obtener la estimación MC
montecarlo=function(datosim,type="mean",z=NULL,alpha=0.95){
  # datosim es un vector con las simulaciones
  # type="mean" -> estimar una media
  # type="prob" -> estimar probabilidad acumulada hasta un punto z
  # z -> Pr(x<=z), z en el espacio de estados de X
  # alpha -> nivel de confianza (por defecto 95%)
  nsim=length(datosim)
  
if(type=="mean") {
  estim=mean(datosim)
  error=sqrt(sum((datosim-estim)^2)/(nsim^2))
  }
else if(type=="prob"){
    Iprob=(datosim<=z)*1
    estim=mean(Iprob)
    error=sqrt(sum((Iprob-estim)^2)/(nsim^2))
  }
else{
    cat("Introduce type='mean' o type='prob'.")
    break
  }
  
ic_low=estim-qnorm((1+alpha)/2)*error
ic_up=estim+qnorm((1+alpha)/2)*error
# return(c(estim=estim,error=error,ic_low=ic_low,ic_up=ic_up))
return(cat("Estimación MC:",estim,"[",ic_low,",",ic_up,"]"))
}

Ejercicios Básicos

Para hacer en clase B: 1,2,3,6,8,9 Simulación : 5

Ejercicio B3.1

En una tienda de animales entran clientes a razón de 8 por hora.

1. ¿Cuál es la probabilidad de que lleguen al menos 4 clientes durante los próximos 30 minutos?
2. ¿Cuál es la probabilidad de que en las 8 horas en que permanece abierta la tienda entren al menos 70 clientes?

lambda=8 # tasa de clientes por hora
# N(t)= número de clientes que entran en una hora PP(lambda)
# -> N(t) ~ Po (t*lambda)
# 1. Pr(N(0.5)>=4)
t=0.5 # media hora -> 
1-ppois(3,lambda*t)

## [1] 0.5665299

# 2. Pr(N(8)>=70) 
t=8  # 8 horas 
1-ppois(69,lambda*t)

## [1] 0.2424104

Ejercicio B3.2

En una estación de bomberos el tiempo entre llamadas por avisos sigue una distribución exponencial con una media de 32 minutos.

1. Acaba de llegar una llamada. ¿Cuál es la probabilidad de que la próxima llamada se produzca en menos de media hora?
2. ¿Cuál es la probabilidad de que se produzcan exactamente dos llamadas durante la próxima hora?

# T= tiempo entre llamadas T ~ Exp(lambda); E(T)=1/lambda
lambda=1/32
# 1. Pr(T < 30)
pexp(30,lambda)

## [1] 0.6083944

# 2. # Pr(N(1)=2)
t=1 # la próxima hora N(1) ~Po(lambda * 1)
dpois(2,lambda*30)

## [1] 0.1720923

Ejercicio B3.3

Una empresa que controla la seguridad de una ciudad ha observado que los intentos de entrar en domicilios ajenos para robar (para los que tienen contratada la seguridad con ellos) ocurren con un PP de tasa 2.2 por día. El sistema opera 24 horas al día.

1. ¿Cuál es la probabilidad de que mañana se produzcan 4 intentos de robo?
2. ¿Cuál es la probabilidad de que no se produzca ningún intento de robo durante la noche (entre las 12 de la noche y las 8 de la mañana)?
3. Si ahora es medianoche y el último intento de robo se produjo a las 10:30pm, ¿cuál es la probabilidad de que el siguiente intento ocurra antes de las 5:30am?

lambda=2.2 # intentos de robo por día
# {Nt;t} : PP(2.2)  --> N(t) ~ Po(lambda*t)

# 1. N(1) ~ Po(lambda*1) ->  Pr(N(1)=4)
dpois(4,lambda)

## [1] 0.1081513

# 2. t= 8 horas=8/24 días ->  Pr(N(t)=0)
t=8/24 # días
dpois(0,lambda/3)

## [1] 0.4803053

# 3. Dos posibilidades son
t=7/24 # duración del periodo entre las 10:30 y las 5:30
# 3.1 Que el siguiente intento sea antes de las 5.30 significa que el número de intentos de robo en un periodo de 7 horas sea al menos de 1
# Pr(N(t)>=1)=1-Pr(N(t)=0)
1-dpois(0,lambda*t)

## [1] 0.4735857

# 3.2 también implica que el tiempo hasta el siguiente intento de robo sea menor a 7 horas
# T tiempo entre intentos de robo ~ Exp(lambda)  -> Pr(T<t)
pexp(t,lambda)

## [1] 0.4735857

Ejercicio B3.4

En el caso de la empresa de seguridad del Ejercicio B3.3 resulta que además se sabe que el 10% de los intentos de entrar en la casa resultan en robo efectivo.

1. ¿Cuál es la probabilidad de que se produzcan 3 intentos de robo que acaben en robo?
2. Ahora mismo son las 6am del lunes. El último intento de robo que acabó en robo ocurrió a medianoche. ¿Cuál es la probabilidad de que no se hayan producido robos durante dicho periodo?

lambda=2.2 # tasa de intentos de robo por día
# {N(t); t} PP(lambda) : número de intentos de robo en un día
# p= proporción de intentos de robo que acaban en robo
p=0.1
# {M(t); t} número de robos efectivos en un día: PP(lambda*p)
# --> M(t) ~Po (lambda*p*t)

# número esperado de robos en un día
lambda*p*1

## [1] 0.22

# 1. Pr(M(1)=3) en un día 
t=1
dpois(3,lambda*p*t)

## [1] 0.001424203

# 2. Pr(M(t)=0), t=6 horas=6/24 días
t=6/24
dpois(0,lambda*p*t)

## [1] 0.9464851

Ejercicio B3.5

En una máquina hay dos tipos comunes de fallos críticos: en la componente electrónica A o en la B. Si cualquiera de estas componentes falla, la máquina se para. La componente A falla conforme a un PP con tasa 1.1 fallos por turno (la fábrica trabaja 24/7 en turnos de 8 horas). La componente B falla según un PP con tasa 1.2 fallos por día.

# Fallos_A : PP(lambda.a)
# Fallos_B : PP(lambda.b)
lambda.a=1.1;lambda.b=1.2/3  # tasas por turno (8 horas)
# Nt=número de fallos (A o B) : PP(lambda.a+lambda.b)
lambda=lambda.a+lambda.b

1. ¿Cuál es la probabilidad de que se produzcan exactamente 5 fallos de la máquina durante un día?

# 1. Pr(Nt=5) t=3 turnos
dpois(5,lambda*3)

## [1] 0.1708269

2. ¿Cuál es la probabilidad de que la máquina no se pare más de una vez durante el siguiente turno?

# 2. Pr(Nt <= 1) t=1 turno
ppois(1,lambda)

## [1] 0.5578254

3. Ahora mismo es mediodía y el último fallo se produjo hace 4 horas. ¿Cuál es la probabilidad de que el próximo parón se produzca antes de las 6pm?

# 3. T=tiempo a fallo ~Exp(lambda). Pr(T<10horas) 10horas=10/8 turnos
pexp(10/8,lambda)

## [1] 0.846645

4. Simula los parones de la máquina, mostrando qué componente falla en cada ocasión y asumiendo que el tiempo de reparación es despreciable y la máquina continua trabajando inmediatamente. Estima la probabilidad de que la máquina no se pare más de 1 vez durante un turno.

Simulación: algoritmo con Poisson compuesto, sin identificación del tipo de avería

1. Simular el número total de fallos por turno (t=1), n ~ Po(lambda_A+lambda_B)
2. Calcular un indicador para la probabilidad solicitada (n<=1)*1
3. Repetir 1-2 nsim veces/turnos para calcular la estimación MC.

nsim=5000
# Nt= nº total de fallos (A y B) hasta t: Po((lambda_A+lambda_B)t)
lambda_A=1.1
lambda_B=1.2/3
lambda=lambda_A+lambda_B
t=1 # 1 turno
set.seed(123)
n=rpois(nsim,lambda)
montecarlo(n<=1)

## Estimación MC: 0.557 [ 0.5432313 , 0.5707687 ]

Simulación: algoritmo con Poisson, identificando el tipo de avería

1. Simular el nº de fallos tipo A (nA) en un turno según Po(lambda_A).
2. Simular el nº de fallos tipo B (nB) en un turno según Po(lambda_B).
3. Construir un indicador que vale 1 si hay a lo sumo 1 fallo en un turno (nA+nB).

Repetir pasos 1-3 nsim veces para tener nsim turnos simulados, y el número de fallos en cada uno. Calcular la aprox MC.

# Asumimos independencia entre los dos tipos de avería: 
# las componentes A y B fallan independientemente la una de la otra
#4 Tr=0 tiempo de reparación
# t=1 -> 1 turno
# N_A(1)~Po(lambda_A) 
lambda_A=1.1
# N_B(1)~Po(lambda_B)
lambda_B=1.2/3
nsim=5000

set.seed(123)
# número de fallos (A o B)
nA=rpois(nsim,lambda_A)
nB=rpois(nsim,lambda_B)
montecarlo(((nA+nB)<=1)*1)  

## Estimación MC: 0.5596 [ 0.5458398 , 0.5733602 ]

Simulación: algoritmo con Exponencial

1. Inicializo sendas simulaciones de tiempo hasta fallo por A y por B, tfallos.
2. Elijo tt el mínimo de los dos (la avería que se produce antes), esto es, el tiempo entre averías consecutivas. Actualizo el reloj sumando este tiempo. Si excedo el tiempo de simulación paro y no registro el fallo. Si no lo excedo, registro el tiempo de fallo y el tipo (componente).
3. Actualizo tiempos hasta fallo: simulo un nuevo tiempo de fallo para el tipo de componente que acaba de fallar; actualizo el tiempo hasta fallo para la otra restando el último tiempo entre fallos tt.
4. Repito 2-3 hasta llegar al tiempo máximo de simulación.

La primera avería la diferencio del resto, pues si el primer tiempo hasta fallo excede el tiempo de simulación, he de devolver un dataframe vacío (NULL).

# simulación de un turno. Averías sin reparaciones.
simula_turno_sinrep = function(){
# Tasas de fallo (por turno)
lambda_A=1.1
lambda_B=1.2/3
# tasas de fallo A y B
lambdas=c(lambda_A,lambda_B)

# tipos de averías
tipos=c("A","B")

# parámetros reparación U(rlow,rup)
rlow=5/480
rup=1/8

# info a almacenar
tipo=c() #tipo de fallo
t=c() #  instantes de fallo

# inicialización
reloj=0 # reloj
TMAX=1 # tiempo máximo (1 turno)

# Primera avería
# simulamos tiempos de fallo
tA=rexp(1,lambdas[1]) #fallo por A
tB=rexp(1,lambdas[2]) #fallo por B
# tiempos hasta fallo para cada componente
tfallos=c(tA,tB)

# identifico qué avería se da primero
minimo=which.min(tfallos)
# identifico el tiempo entre averías consecutivas
tt=tfallos[minimo]

  # Si excedemos TMAX salimos sin averías
  if(tt>TMAX){
    return(data.frame(t=NULL,tipo=NULL))
  }
  # en otro caso, identificamos el tipo de fallo 
  # y actualizamos reloj
  else{
  i=1 # primer fallo 
  t[i]=tt # instante en que se produce (reloj)
  tipo[i]=tipos[minimo] # tipo de fallo
  }

# Siguientes averías (>1) hasta que lleguemos a TMAX
while(t[i]<TMAX){
  # simulo un nuevo tiempo hasta fallo en la componente averiada
  tfallos[minimo]=rexp(1,lambdas[minimo])
  # actualizo el tiempo hasta fallo en la componente no averiada
  tfallos[-minimo]=tfallos[-minimo]-tt
  # identifico qué avería se da antes
  minimo=which.min(tfallos)
  # calculo el tiempo hasta fallo
  tt=tfallos[minimo]
  # si sobrepaso TMAX, salgo sin registrar avería
  if((t[i]+tt)>TMAX)
    break
  # en otro caso actualizo iteración, instante de fallo (reloj) y tipo de avería
  else{
   i=i+1 # actualizamos la iteración
   t[i]=t[i-1]+tt
   tipo[i]=tipos[minimo]
  }
}  

return(data.frame(t,tipo))
}

# Aprox MC
nsim=5000
n_averias=c()
for(i in 1:nsim){
  turno=simula_turno_sinrep()
  n_averias[i]=nrow(turno)
}
montecarlo((n_averias<=1)*1)  

## Estimación MC: 0.5558 [ 0.5420275 , 0.5695725 ]

5. Simula los parones de la máquina, mostrando qué componente falla en cada ocasión y asumiendo que el tiempo de reparación de cada componente se distribuye uniforme entre 5 y 60 minutos. Estima la probabilidad de que la máquina no se pare más de 1 vez durante un turno. Estima la probabilidad de que no haya más de una avería tipo A durante un turno. Haz lo mismo con las averías tipo B.

Simulación: algoritmo con Exponencial

1. Inicializo sendas simulaciones de tiempo hasta fallo por A y por B, tfallos.
2. Elijo tt el mínimo de los dos (la avería que se produce antes), esto es, el tiempo entre averías consecutivas. Actualizo el reloj sumando este tiempo. Si excedo el tiempo de simulación paro y no registro el fallo. Si no lo excedo, registro el tiempo de fallo y el tipo (componente).
3. Simulo un tiempo de reparación. Actualizo el reloj sumando este tiempo y verifico que no excedo el periodo de simulación. Si lo excedo paro sin registrar, si no, sigo.
4. Actualizo tiempos hasta fallo: simulo un nuevo tiempo de fallo para el tipo de componente que acaba de fallar; actualizo el tiempo hasta fallo para la otra restando el último tiempo entre fallos tt y añadiendo el último tiempo de reparación.
5. Repito 2-3 hasta llegar al tiempo máximo de simulación.

La primera avería la diferencio del resto, pues si el primer tiempo hasta fallo excede el tiempo de simulación, he de devolver un dataframe vacío (NULL).

# simulación de un turno. Averías sin reparaciones.
simula_turno_conrep = function(){

# Tasas de fallo (por turno)
lambda_A=1.1
lambda_B=1.2/3
# tasas de fallo A y B
lambdas=c(lambda_A,lambda_B)

# tipos de averías
tipos=c("A","B")

# parámetros reparación U(rlow,rup)
rlow=5/480
rup=1/8

# info a almacenar
tipo=c() #tipo de fallo
t=c() #  instantes de fallo
trep=c() # tiempos de reparación

# inicialización
reloj=0 # reloj
TMAX=1 # tiempo máximo (1 turno)

# Primera avería
# simulamos los tiempos hasta la 1ª avería
tA=rexp(1,lambdas[1]) #fallo por A
tB=rexp(1,lambdas[2]) #fallo por B
# tiempos hasta fallo para cada componente
tfallos=c(tA,tB)

# identifico qué avería se da primero
minimo=which.min(tfallos)
# tiempo transcurrido (reloj)
tt=tfallos[minimo]

  # Si excedemos TMAX salimos sin eventos
  if(tt>TMAX){
    return(data.frame(t=NULL,tipo=NULL,trep=NULL))
  }
  # en otro caso, identificamos el tipo de fallo 
  # actualizamos reloj y simulamos tpo. reparación
  else{
  i=1
  t[i]=tt
  tipo[i]=tipos[minimo]
  reloj=reloj+tt
  # simulo tiempo reparación
  trep[i]=runif(1,rlow,rup)
  reloj=reloj+trep[i]
   # si el tiempo de reparación sobrepasa TMAX, devolvemos TMAX-t[i] y salimos
    if(reloj>=TMAX){
    trep[i]=TMAX-t[i]
    return(data.frame(t,tipo,trep))
     }
  }

# Siguientes averías (>1) hasta que lleguemos a TMAX
while(reloj<TMAX){
  # simulo nuevo tiempo hasta fallo en la componente averiada
  tfallos[minimo]=rexp(1,lambdas[minimo])
  # actualizo el tiempo hasta fallo en la componente no averiada
  # durante la reparación la máquina estaba parada
  # y no consumía tiempo de funcionamiento
  tfallos[-minimo]=tfallos[-minimo]-tt+trep[i]
  # identifico qué avería se da antes
  minimo=which.min(tfallos)
  # calculo el tiempo hasta fallo
  tt=tfallos[minimo]
  #actualizo reloj
  reloj=t[i]+tt
  # si sobrepaso TMAX, salgo sin registrar fallo
  if(reloj>TMAX)
    break
  # en otro caso actualizo iteración, reloj, tiempos y tipo de avería
  else{
   i=i+1 # actualizamos la iteración: nuevo fallo
   t[i]=t[i-1]+trep[i-1]+tt
   tipo[i]=tipos[minimo]
   trep[i]=runif(1,rlow,rup)
   reloj=t[i]+trep[i]
   # si el tiempo de reparación sobrepasa TMAX, devolvemos TMAX-tpo.fallo y salimos
    if(reloj>=TMAX){
    trep[i]=TMAX-t[i]
    break
      }
  }
}  

return(data.frame(t,tipo,trep))
}

# probabilidad de que la máquina no se pare más de 1 vez durante un turno

nsim=1000
n_averias=c() # nº total de averías
n_averias_A=c() # nº total de averías A
n_averias_B=c() # nº total de averías B

for(i in 1:nsim){
  turno=simula_turno_conrep()
  n_averias[i]=nrow(turno)
  n_averias_A[i]=sum(turno$tipo=="A")
  n_averias_B[i]=sum(turno$tipo=="B")
}
# Pr(N(1)<=1)
montecarlo(n_averias<=1)  

## Estimación MC: 0.563 [ 0.5322572 , 0.5937428 ]

# Pr(NA(1)<=1)
montecarlo(n_averias_A<=1)  

## Estimación MC: 0.721 [ 0.6932017 , 0.7487983 ]

# Pr(NB(1)<=1)
montecarlo(n_averias_B<=1)   

## Estimación MC: 0.96 [ 0.9478545 , 0.9721455 ]