# 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
cat("\n Estimación MC de una ",ifelse(type=="mean","media","probabilidad"),estim,"[",ic_low,",",ic_up,"] \n")
return(c(estim=estim,error=error,ic_low=ic_low,ic_up=ic_up))
}Un proceso de fabricación tiene una tasa de defectos del 20% y los artículos se colocan en cajas de cinco. Un inspector toma muestras de dos artículos de cada caja. Si uno o los dos artículos seleccionados son defectuosos, la caja se rechaza. Si un cliente pide 10 cajas, ¿cuál es el número esperado de artículos defectuosos que recibirá el cliente? Da una aproximación del error y una banda de confianza.
Algoritmo
Repetir los pasos 2 a 4 hasta tener 10 cajas comercializadas (llegan al cliente):
Repetir el proceso nsim veces para obtener una aproximación MC con ndefectos.
nsim=10000
p=0.2
ncajas=10
npiezas=5
ninsp=2
ndefectos=c() # número de defectos en ncajas
for(i in 1:nsim){
ncajas_com=0 # inicialización num.cajas comercializadas
ndefectos[i]=0
while(ncajas_com<=ncajas){
caja=rbinom(npiezas,1,p)
inspector=sample(caja,2)
comercializada=ifelse(sum(inspector)>0,FALSE,TRUE)
if(comercializada){
ncajas_com=ncajas_com+1
ndefectos[i]=ndefectos[i]+sum(caja)
}
} # ncajas comercializadas
} # nsim simulaciones
simulacion=tibble(id=1:nsim,ndefectos)
montecarlo(simulacion$ndefectos)##
## Estimación MC de una media 6.5933 [ 6.54782 , 6.63878 ]## estim error ic_low ic_up
## 6.59330000 0.02320451 6.54781999 6.63878001Un vendedor a domicilio vende ollas y sartenes. Sólo entra en el 50% de las casas que visita. De las casas en las que entra, 1/6 de los propietarios no están interesados en comprar nada, 1/2 de ellos acaban haciendo un pedido de 60 dólares, y 1/3 de ellos acaban haciendo un pedido de 100 dólares. Estima el promedio de ventas por visita. Estima el promedio de ventas (en dólares) conseguidas para 25 visitas (junto con una medida del error).
Información
Casas que visita
Algoritmo
Para estimar una venta media, repetimos nsim veces (una venta) y calculamos la aprox.MC.
nsim=5000
nvisitas_dia=25
p=c(0.5,0.5*1/6,0.5*1/2,0.5*1/3)
valor_venta=c(0,0,60,100)
pacum=cumsum(p)
# Aproximación venta por visita: simulación de nsim ventas
ventas_visita=c() # vector para registras las nsim ventas
for(i in 1:nsim){
u=runif(1)
# por la transformada inversa simulamos una venta a partir de u
ventas_visita[i]=valor_venta[min(which(pacum>=u))]
}
montecarlo(ventas_visita)##
## Estimación MC de una media 32.676 [ 31.57617 , 33.77583 ]## estim error ic_low ic_up
## 32.6760000 0.5611486 31.5761689 33.7758311# Aproximación venta total por jornada/día (25 visitas):
# simulación de nsim jornadas (cada una con 25 visitas)
ventas_dia=c()
for(i in 1:nsim){
ventas_dia[i]=0
for(j in 1:nvisitas_dia){
u=runif(1)
venta_visita=valor_venta[min(which(pacum>=u))]
ventas_dia[i]=ventas_dia[i]+venta_visita
}
}
# venta media por día (25 visitas)
montecarlo(ventas_dia)##
## Estimación MC de una media 788.804 [ 783.3583 , 794.2497 ]## estim error ic_low ic_up
## 788.804000 2.778455 783.358329 794.249671# SIN BUCLES. Ventas por jornada
# Función para simular una venta a partir de un valor uniforme
# por el método de la transformada inversa
fun=function(u){
pacum=cumsum(p)
venta=valor_venta[which.min(pacum>u)]
return(venta)
}
# número de jornadas
nsim=5000
nvisitas_dia=25
# se simulan a la vez todas las uniformes
u=runif(nsim*nvisitas_dia)
simulaciones=data.frame(u=u)
# y de ellas todas las ventas
simulaciones$venta=apply(simulaciones,1,fun)
# creamos un identificador de jornada (simulación)
simulaciones$dia=as.factor(rep(1:nsim,nvisitas_dia))
# agrupamos por jornada (día) y acumulamos las ventas de las 25 visitas
ventas_dia=simulaciones %>%
group_by(dia) %>%
summarise(ventas=sum(venta))
# estimamos por MC las ventas diarias
montecarlo(ventas_dia$ventas)##
## Estimación MC de una media 0 [ 0 , 0 ]## estim error ic_low ic_up
## 0 0 0 0Un contratista de techos independiente ha determinado que el número de trabajos que se solicitan para el mes de septiembre es bastante variable. A partir de la experiencia anterior, las probabilidades de obtener 0, 1, 2 o 3 trabajos se han determinado como 0.1, 0.35, 0.30 y 0.25, respectivamente. El beneficio obtenido por cada trabajo es de 300 dólares. ¿Cuál es el beneficio esperado para el mes de septiembre?
Algoritmo
Repetir nsim veces la simulación de un mes (septiembre) de trabajo
Calcular la aprox.MC a partir del vector de beneficios almacenado.
p = c(0.1, 0.35, 0.30, 0.25) # Probabilidad de encargar trabajos
ntrabajos = 0:3
beneficio = 300 # Beneficio por trabajo
# Acumulados para transformación inversa
pacum <- cumsum(p)
nsim=5000
beneficio_mes=c()
for(i in 1:nsim){
u=runif(1)
# Trabajos solicitados
ntrabajos_solic <- ntrabajos[min(which(u<= pacum))]
# Beneficio mes de septiembre
beneficio_mes[i]=ntrabajos_solic*beneficio
}
montecarlo(beneficio_mes)##
## Estimación MC de una media 514.08 [ 506.1355 , 522.0245 ]## estim error ic_low ic_up
## 514.080000 4.053388 506.135506 522.024494Un sistema de visión está diseñado para medir el ángulo en el que el brazo de un robot se desvía de la vertical. Sin embargo, el sistema de visión no es totalmente preciso, y el resultado de las observaciones es una variable aleatoria con una distribución uniforme. Si la medición indica que el rango del ángulo está entre 9.7 y 10.5 grados, aproxima por simulación la probabilidad de que el ángulo real esté entre 9.9 y 10.1 grados y da una medida del error de dicha aproximación.
Algoritmo
# Fijamos simulaciones
nsim <- 10000
# Generamos datos uniformes
angulo <- runif(nsim, 9.7, 10.5)
# creamos el indicador
indicador <- (angulo <=10.1 & angulo >9.7)*1
estim=mean(indicador)
error=sd(indicador)/sqrt(nsim)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de una probabilidad",estim,"[",ic_low,",",ic_up,"] \n")##
## Estimación MC de una probabilidad 0.4959 [ 0.4861 , 0.5057 ]En una operación de soldadura automatizada, la posición en la que se coloca la soldadura es muy importante. La desviación del centro de la placa es una variable aleatoria normal con una media de 0 pulgadas y una desviación estándar de 0.01 pulgadas. Una desviación positiva indica una desviación a la derecha del centro y una desviación negativa indica una desviación a la izquierda del centro.
¿Cuál es la probabilidad de que en una placa dada, la ubicación real de la soldadura se desvíe menos de 0.005 pulgadas (en valor absoluto) del centro?
¿Cuál es la probabilidad de que en una placa dada, la ubicación real de la soldadura se desvíe más de 0.02 pulgadas (en valor absoluto) del centro?
Da aproximaciones del error cometido al aproximar por simulación.
Algoritmo
Simular nsim desviaciones de la soldadura: D ~ N(0,0.01)
nsim=5000
d=rnorm(nsim,0,0.01)
# Pr(|D|<0.005)
indicador=(abs(d)<0.005)*1
estim=mean(indicador)
error=sd(indicador)/sqrt(nsim)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de Pr(|D|<0.005)",estim,"[",ic_low,",",ic_up,"] \n")##
## Estimación MC de Pr(|D|<0.005) 0.389 [ 0.3754854 , 0.4025146 ]# Pr(|D|>0.02)
indicador=(abs(d)>0.02)*1
estim=mean(indicador)
error=sd(indicador)/sqrt(nsim)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de Pr(|D|>0.02)",estim,"[",ic_low,",",ic_up,"] \n")##
## Estimación MC de Pr(|D|>0.02) 0.0474 [ 0.04150952 , 0.05329048 ]Un departamento de compras percibe que el 75% de sus pedidos especiales se reciben a tiempo. De los pedidos que se reciben a tiempo, el 80% cumple totalmente las especificaciones; de los pedidos que llegan con retraso, el 60% cumple con las especificaciones. Responde a las siguientes preguntas utilizando simulación, dando también una medida del error.
¿Cuál es la probabilidad de que un pedido llegue a tiempo y cumpla con las especificaciones?
¿Cuál es la probabilidad de que un pedido cumpla con las especificaciones?
Si se han recibido cuatro pedidos, ¿cuál es la probabilidad de que los cuatro pedidos cumplan con las especificaciones?
Información
PEDIDOS
Algoritmo
Repetir nsim veces
Los conteos de las probabilidades implican generar un indicador 0/1 que vale 1 si se da la condición, promediar dicho indicador (para obtener la estimación MC de la probabilidad) y calcular su desviación típica para estimar el error MC e IC.
nsim=4*5000
p_pedido=c(0.75,0.25)
pacum_pedido=cumsum(p_pedido)
p_espec_llega=c(0.8,0.2)
pacum_espec_llega=cumsum(p_espec_llega)
p_espec_nollega=c(0.6,0.4)
pacum_espec_nollega=cumsum(p_espec_nollega)
simulaciones=data.frame(u=runif(nsim),llega=NA,espec=NA)
for(i in 1:nsim){
llegada=ifelse(simulaciones$u[i]<pacum_pedido[1],"llega","no llega")
simulaciones$llega[i]=llegada
v=runif(1)
if(llegada=="llega")
espec=ifelse(v<pacum_espec_llega[1],"cumple","no cumple")
else
espec=ifelse(v<pacum_espec_nollega[1],"cumple","no cumple")
simulaciones$espec[i]=espec
}
head(simulaciones)## u llega espec
## 1 0.3454662 llega no cumple
## 2 0.9287886 no llega no cumple
## 3 0.6085737 llega cumple
## 4 0.4896798 llega cumple
## 5 0.6509288 llega cumple
## 6 0.4769310 llega cumple# Pr(pedido=a tiempo & especificaciones=cumple)
probs = simulaciones %>%
mutate(indi1=(llega=="llega" & espec=="cumple")*1,
#Pr(especificaciones=cumple)
indi2=(espec=="cumple")*1) %>%
summarise(estim1=mean(indi1),sd1=sd(indi1),
estim2=mean(indi2),sd2=sd(indi2))
probs## estim1 sd1 estim2 sd2
## 1 0.60015 0.4898796 0.7473 0.4345712# Pr(4pedidos cumplan especificaciones)
p4=simulaciones %>%
# agrupamos en grupos de 4 pedidos (pedido4) y creamos un indicador que vale 1 si cada pedido cumple especificaciones
mutate(pedido4=gl(5000,4,length=4*5000), cumple=(espec=="cumple")*1)%>%
# agrupamos en grupos de 4 pedidos y contamos cuántos cumplen espec en cada grupo
group_by(pedido4) %>%
summarise(cumplen=sum(cumple))%>%
# creamos un indicador que vale 1 si cumplen espec los 4 pedidos de cada grupo
mutate(indicador=(cumplen==4)*1)
p4## # A tibble: 5,000 × 3
## pedido4 cumplen indicador
## <fct> <dbl> <dbl>
## 1 1 2 0
## 2 2 4 1
## 3 3 2 0
## 4 4 4 1
## 5 5 4 1
## 6 6 2 0
## 7 7 3 0
## 8 8 2 0
## 9 9 1 0
## 10 10 3 0
## # … with 4,990 more rows# aproximamos por MC la Pr(cumplan 4 pedidos)
estim=mean(p4$indicador)
error=sd(p4$indicador)/sqrt(5000)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de Pr(cumplan espec 4 pedidos)",estim,"[",ic_low,",",ic_up,"] \n")##
## Estimación MC de Pr(cumplan espec 4 pedidos) 0.3096 [ 0.2967839 , 0.3224161 ]Una empresa manufacturera tiene tres operarios para una máquina que produce cierto tipo de componentes. El operario A tiene una tasa de defectos del 5%, el operario B del 3%, y el operario C del 2%. Los tres operarios producen el mismo número de componentes. Si un componente elegido al azar resulta defectuoso, ¿cuál es la probabilidad de que el componente haya sido producido por A, B, o C? Responde a la pregunta con simulación y da una medida del error de la aproximación.
nsim=10000
operarios=c("A","B","C")
oper=sample(operarios,size=nsim,replace=TRUE)
head(oper)## [1] "B" "A" "A" "B" "C" "A"p=c(1/3,1/3,1/3)
pacum=cumsum(p)
u=runif(nsim)
oper=c()
defectos=c()
for(i in 1:nsim){
oper[i]=operarios[min(which(pacum>=u[i]))]
if(oper[i]=="A"){
pdefectos=0.05}
else if(oper[i]=="B"){
pdefectos=0.03}
else{
pdefectos=0.02}
#pdefectos=ifelse(oper[i]=="A",0.05,ifelse(oper[i]=="B",0.03,0.02))
defectos[i]=rbinom(1,1,pdefectos)
}
simulaciones=data.frame(oper,defectos)
table(defectos)## defectos
## 0 1
## 9656 344# Pr(defecto)
estim=mean(simulaciones$defectos)
error=sd(simulaciones$defectos)/sqrt(nsim)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de Pr(defecto)",estim,"[",ic_low,",",ic_up,"]")##
## Estimación MC de Pr(defecto) 0.0344 [ 0.0308277 , 0.0379723 ]#PR(A)
indicador=(simulaciones$oper=="A")*1
estim=mean(indicador)
error=sd(indicador)/sqrt(nsim)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de Pr(operador)",estim,"[",ic_low,",",ic_up,"]")##
## Estimación MC de Pr(operador) 0.3353 [ 0.3260466 , 0.3445534 ]# Pr(A|Defectuoso)
table(simulaciones)## defectos
## oper 0 1
## A 3170 183
## B 3224 97
## C 3262 64selec=which(simulaciones$defectos==1)
productor=simulaciones$oper[selec ]
indicador=(simulaciones$oper[selec]=="A")*1
estim=mean(indicador)
error=sd(indicador)/sqrt(length(indicador))
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("\n Estimación MC de Pr(A|defectos)",estim,"[",ic_low,",",ic_up,"]")##
## Estimación MC de Pr(A|defectos) 0.5319767 [ 0.479171 , 0.5847825 ]El 1% de los préstamos que hace cierta empresa financiera no son saldados (es decir, la cantidad prestada no es devuelta en su totalidad). La compañía efectúa un estudio rutinario de las posibilidades crediticias de los solicitantes. Encuentra que el 30% de los préstamos no saldados se hicieron a clientes de alto riesgo, el 40% a clientes de riesgo moderado y el restante 30% a clientes de bajo riesgo. De los préstamos que fueron saldados, el 10% se hicieron a clientes de alto riesgo, el 40% a clientes de riesgo moderado y el 50% a clientes de bajo riesgo. Responde a las siguientes preguntas utilizando simulación, dando también una medida del error.
¿Cuál es la probabilidad de que un préstamo de alto riesgo no sea saldado?
¿Cuál es la probabilidad de que una deuda no sea saldada, dado que el riesgo del cliente es moderado?
Información
PRÉSTAMOS
Algoritmo 1. Calcular las probabilidades marginales del árbol para cada combinación Saldo/Cliente: p1,p2,…,p6 2. Simular pares de combinaciones Saldo/Cliente con las probabilidades anteriores. 3. \(1\) Pr(NS|RA)=Pr(RA|NS)Pr(NS)/Pr(RA)=(0.3x0.01)/0.102, 1. Pr(RA)=Pr(RA|NS)Pr(NS)+Pr(RA|S)Pr(S)=0.3*0.01+0.1*0.99=0.102
2. Aproximar por MC con casos favorables Saldo==NS y casos posibles Cliente==RA