# Complementaria Simulacion -----------------------------------------------
library(markovchain)
## Package: markovchain
## Version: 0.8.6
## Date: 2021-05-17
## BugReport: https://github.com/spedygiorgio/markovchain/issues
library(formattable)
## Warning: package 'formattable' was built under R version 4.1.3
library(plotly)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.1.3
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:formattable':
##
## style
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
estados<-c(1:3)
matP<-matrix(c(0.5,0.3,0.2,
0.3,0.4,0.3,
0.2,0.3,0.5),nrow = 3,ncol = 3,byrow = TRUE,dimnames = list(estados,estados))
matP
## 1 2 3
## 1 0.5 0.3 0.2
## 2 0.3 0.4 0.3
## 3 0.2 0.3 0.5
#Muestra aleatoria de pasar a un estado empezando en estado 1
estado_actual=1
sample(x = estados,size = 1,prob = matP[estado_actual,])
## [1] 1
#Simular varios dias
set.seed(0) #Definir semillas para valores aleatorios
dias<-10 #Dias a simular el nivel de ICA
simulacion_estados<-c() #Inicalizar lista de simulacion de nivel ICA en proximos 20 dias
for (i in 1:dias) { #Simular dias
simulacion_estados<-c(simulacion_estados,estado_actual)
estado_actual<-sample(x = estados,size = 1,prob = matP[estado_actual,])
}
simulacion_estados #Simulacion del Nivel ICA
## [1] 1 3 3 3 2 1 1 3 1 2
#Graficar replica
plot_ly(x=c(1:dias), y=simulacion_estados, type="scatter", mode="lines", name="Data")%>%
layout(title="Simulacion Nivel ICA",
xaxis=list(title="Dias"),
yaxis=list(title="Nivel del ICA"))
#Simular 10 replicas
escenarios<-10
matriz_calidad<-matrix(0,nrow = dias,ncol=escenarios)
for (e in 1:escenarios) {
estado_actual=1
for (i in 1:dias) {
matriz_calidad[i,e]<-estado_actual
estado_actual<-sample(x = estados,size = 1,prob = matP[estado_actual,])
}
}
#Graficar replica
fig<-plot_ly(x=c(1:dias), y=matriz_calidad[,1],type="scatter", mode="lines")%>%
layout(title="Simulacion Nivel ICA",
xaxis=list(title="Dias"),
yaxis=list(title="Nivel del ICA"))
for (e in 2:escenarios) {
fig<-fig%>%add_trace(y=matriz_calidad[,e],mode = "lines")
}
fig
# Ejercicio ---------------------------------------------------------------
# Literal A ---------------------------------------------------------------
#Tasa de la demanda
lambda<-15
#Modelación del manejo de inventarios de Wattenspharma con cadenas de Markov
#Crear los estados
estados<-c(0:100)
#Crear y llenar la matriz P de la política de inventarios
matrizP<-matrix(0,nrow = length(estados), ncol = length(estados),dimnames = list(estados,estados))
#Para la Politica -> si i<=70 solicita 30 cajas
for (i in estados) {
for (j in estados) {
if(i<=70 & j>0){
matrizP[i+1,j+1]= dpois(30+i-j,lambda=lambda)
}else if(i<=70 & j==0){
matrizP[i+1,j+1]=ppois(30+i-1,lambda = lambda, lower.tail = F)
}else if(i>70 & j>0){
matrizP[i+1,j+1]=dpois(i-j,lambda=lambda)
}else if(i>70 & j==0){
matrizP[i+1,j+1]=ppois(i-1,lambda = lambda, lower.tail = F)
}
}
}
# Crear la cadena usando el paquete markovchain
cmtd<-new("markovchain", states=as.character(estados), transitionMatrix=matrizP)
# Literal B ---------------------------------------------------------------
#B.1 Valor esperado de inventario proximas semanas
# Vector para guardar el valor esperado del inventario en cada semana
vEsperadoInv <- c()
# Recorrer las semanas
for(n in 1:30){
estado_inicial <- "0" # Definir el estado inicial
P_t <- cmtd^n # Calcular la matriz de probabilidades en el transitorio
probs <- P_t[estado_inicial,] # Calcular probabilidades en el transitorio dado el estado inicial
vEsperadoInv <- c(vEsperadoInv, probs%*%estados) # Calcular valor esperado y agregar al vector
}
plot_ly(y=vEsperadoInv,type="scatter", mode="lines")%>%
layout(title="Valor esperado del inventario",
xaxis=list(title="Dias"),
yaxis=list(title="Inventario"))
#B.2 Costo esperado ordenar durante 30 sem
# Calcula el valor esperado del costo de ordenar de las siguientes 30 semanas
cOrdenar <- 5000
vEsperadoOrd <- 0
for(n in 1:30){
estado_inicial <- "0" # Definir el estado inicial
P_t <- cmtd^n # Calcular la matriz de probabilidades en el transitorio
probs <- P_t[estado_inicial,] # Calcular probabilidades en el transitorio dado el estado inicial
# Calcular el valor esperado de ordenar y sumarlo al costo hasta el momento
vEsperadoOrd <- vEsperadoOrd + sum(probs[1:71])*cOrdenar
}
vEsperadoOrd
## [1] 84249.97
# C. Simulacion inventario ------------------------------------------------
escenarios <- 8
semanas <- 30
# Inicializar matriz para almacenar la simulación del inventario para cada escenario
inve_escenario <- matrix(0, nrow=semanas, ncol=escenarios)
# Simular la cadena de Markov
set.seed(0)
for(j in 1:escenarios){
# Definir el estado inicial (inventario inicial)
estado <- 0
# Simular las transiciones por 30 semanas
for(i in 1:semanas){
# Guardar el inventario actual
inve_escenario[i,j] <- estado
# Obtener estado futuro
estado <- sample(estados, size=1, prob=matrizP[as.character(estado),])
}
}
inve_escenario
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 0 0 0 0 0 0 0
## [2,] 9 17 13 15 19 21 9 11
## [3,] 23 35 21 29 23 33 23 18
## [4,] 40 47 35 44 41 51 39 34
## [5,] 52 65 48 63 59 68 48 54
## [6,] 61 80 65 84 75 81 66 76
## [7,] 77 60 84 74 56 69 87 56
## [8,] 56 79 68 64 69 84 72 73
## [9,] 78 69 81 77 87 72 62 62
## [10,] 67 84 71 60 71 57 73 69
## [11,] 79 65 56 70 55 71 65 76
## [12,] 65 78 77 82 67 55 83 63
## [13,] 79 58 64 71 79 69 64 77
## [14,] 64 77 74 58 65 78 81 63
## [15,] 83 67 61 72 81 61 67 80
## [16,] 70 85 78 46 70 81 75 68
## [17,] 90 73 61 58 78 60 59 76
## [18,] 78 63 70 72 60 73 71 64
## [19,] 59 79 91 58 78 58 57 81
## [20,] 63 62 78 71 66 75 67 67
## [21,] 80 73 68 49 71 56 84 80
## [22,] 70 62 75 61 59 73 74 71
## [23,] 92 75 58 67 78 55 58 58
## [24,] 76 66 69 78 60 65 72 73
## [25,] 65 79 82 65 74 86 60 60
## [26,] 80 63 69 78 58 73 74 72
## [27,] 64 79 89 64 69 60 59 59
## [28,] 81 65 75 79 82 69 77 78
## [29,] 67 82 56 60 67 88 65 67
## [30,] 84 70 71 76 78 69 80 85
#Graficar replica
fig<-plot_ly(x=c(1:semanas), y=inve_escenario[,1],type="scatter", mode="lines")%>%
layout(title="Simulacion Inventario",
xaxis=list(title="Semana"),
yaxis=list(title="Inventario"))
for (e in 2:escenarios) {
fig<-fig%>%add_trace(y=inve_escenario[,e],mode = "lines")
}
fig
# D. Simulacion inventario ------------------------------------------------
escenarios <- 1000
semanas <- 30
# Inicializar matriz para almacenar la simulación del inventario para cada escenario
inve_escenario <- matrix(0, nrow=semanas, ncol=escenarios)
# Simular la cadena de Markov
set.seed(0)
for(j in 1:escenarios){
# Definir el estado inicial (inventario inicial)
estado <- 0
# Simular las transiciones por 30 semanas
for(i in 1:semanas){
# Guardar el inventario actual
inve_escenario[i,j] <- estado
# Obtener estado futuro
estado <- sample(estados, size=1, prob=matrizP[as.character(estado),])
}
}
plot_ly(x = inve_escenario[30,], type = "histogram")%>%
layout(title="Histograma Inventario",xaxis=list(title="Inventario"),yaxis=list(title="Frecuencia"))
#Valor esperado inventario en dia 30 con simulacion
Inventario30<-mean(inve_escenario[30,])
Inventario30
## [1] 70.567
#Valor esperado inventario en dia 30 con transitorio
vEsperadoInv[30]
## [1] 70.49893
# E. Simulacion costo ordenar ------------------------------------------------
escenarios <- 1000
semanas <- 30
# Inicializar matriz para almacenar la simulación del inventario para cada escenario
cOrdenar <- c()
# Simular la cadena de Markov
set.seed(0)
for(j in 1:escenarios){
# Definir el estado inicial (inventario inicial)
estado <- 0
costo<-0
# Simular las transiciones por 30 semanas
for(i in 1:semanas){
# Evaluar si costo
if(estado<=70){
costo=costo+5000
}
# Obtener estado futuro
estado <- sample(estados, size=1, prob=matrizP[as.character(estado),])
}
cOrdenar<-c(cOrdenar,costo)
}
plot_ly(x = cOrdenar, type = "histogram")%>%
layout(title="Histograma costo ordenar",xaxis=list(title="Costo"),yaxis=list(title="Frecuencia"))
#Valor esperado costo ordenar con simulacion
Inventario30<-mean(cOrdenar)
Inventario30
## [1] 86550
#Valor esperado cOrdenar con transitorio simulacion
vEsperadoOrd
## [1] 84249.97