SDP

#Numero de epocas
Tmax = 3

#Estados 
estados = c('Excelente','Bueno','Promedio','Malo')

#Decisiones
decisiones = c("Reemplazar",'No reemplazar')

#Retornos inmediatos

datos_retornos = c(-1e7,100,
                   -100,80,
                   -100,50,
                   -100,10)

retorno_inmediato = matrix(datos_retornos,byrow = T, nrow = 4, ncol = 2)
dimnames(retorno_inmediato) = list(estados,decisiones)


#Estado futuro - Matrices de probabilidad

#Decisión de reemplazar
datos_Reemplazar = c(0  ,0  ,0,0,
                     0.7,0.3,0,0,
                     0.7,0.3,0,0,
                     0.7,0.3,0,0)

p_Reemplazar = matrix(datos_Reemplazar, byrow = T, ncol = 4, nrow = 4)
dimnames(p_Reemplazar) = list(estados,estados)


#Decisión de NO reemplazar

datos_Noreemplazar = c(0.7,0.3,0,0,
                       0,0.7,0.3,0,
                       0,0,0.6,0.4,
                       0,0,0,1)

p_NO_reemplazar = matrix(datos_Noreemplazar, byrow = T, ncol = 4, nrow = 4)
dimnames(p_NO_reemplazar) = list(estados,estados)


#Matriz con las ecuaciones de Bellmann
Ft_optimo = matrix(0,nrow = length(estados), ncol = Tmax)
dimnames(Ft_optimo) = list(estados,seq(1:Tmax))


#Matriz con las decisiones optimas
Mat_Dec_Optimo = matrix(0,nrow = length(estados), ncol = Tmax)
dimnames(Mat_Dec_Optimo) = list(estados,seq(1:Tmax))

#---------------------Recursión hacia atras de la semana 3-----------------------------

#Crear las funciones de Bellman f_3(i)
f_i = matrix(0, nrow = length(estados), ncol = length(decisiones))
dimnames(f_i) = list(estados,decisiones)

#En la ultima epoca la ecuación de Bellman es el retorno inmediato
f_i = retorno_inmediato


#Escoger la mejor decisión para cada estado
Ft_optimo[,Tmax] = apply(f_i,1,max)
#apply es una función que aplica una función especificada a las filas 
#o columnas de una matriz o DataFrame.

#f_i es la matriz o DataFrame a la que se aplica la función.

#1 indica el margen a lo largo del cual se aplica la función. 
#En este caso, 1 indica que la función debe aplicarse a las filas. 
#Si fuera 2, la función se aplicaría a las columnas.

#max es la función que se va a aplicar.

#Por lo tanto, apply(f_i, 1, max) aplicará la función 
#max a cada fila de la matriz o DataFrame f_i. 
#El resultado será un vector que contiene el valor máximo de cada fila.

Mat_Dec_Optimo[,Tmax] = decisiones[apply(f_i,1,which.max)]

#----------------------Recursion hacia atras------------------------------#


for(t in (Tmax-1):1){
  
  # Crear las funciones de Bellman - guarda la ecuacion
  f_i<- matrix(0,nrow=length(estados),ncol=length(decisiones))
  dimnames(f_i) <- list(estados,decisiones)
  
  f_i[,"Reemplazar"] <- retorno_inmediato[,"Reemplazar"] + p_Reemplazar%*%Ft_optimo[,t+1]
  f_i[,"No reemplazar"] <- retorno_inmediato[,"No reemplazar"] + p_NO_reemplazar%*%Ft_optimo[,t+1]
  
  # Maximo
  Ft_optimo[,t] <- apply(f_i, 1, max)
  Mat_Dec_Optimo[,t]<- decisiones[apply(f_i, 1, which.max)]
  
}


# Mapa de Calor -----------------------------------------------------------


library(magrittr)
library(dplyr)
library(tidyr)
library(ggplot2)

Mat_Dec_Optimo <- data.frame(Mat_Dec_Optimo)
dimnames(Mat_Dec_Optimo) <- list(estados,seq(1:Tmax))

# Unir todas las épocas y decisiones
Mat_Dec_Optimo <- gather(Mat_Dec_Optimo, Epocas, Decision)

# rep genera una lista de valores iguales a los estados
Mat_Dec_Optimo$Estado <- rep(estados,Tmax)

# Convertir a número
Mat_Dec_Optimo$Epocas <- as.numeric(Mat_Dec_Optimo$Epocas)

# Graficar
ggplot(data = Mat_Dec_Optimo, aes(x = Epocas, y = Estado, fill = Decision)) + geom_tile()+ ylim("Malo","Promedio","Bueno","Excelente") + labs(x = "Epocas")