Se lleva a cabo un análisis de mercado para conocer las preferencias de compras de coches en formato “renting”, según el cual cada año se renueva el coche a cada cliente del servicio en el mes de enero. La empresa está interesada en conocer si los clientes cambian el tipo de vehículo entre las tres opciones posibles (“sedan”, “station wagon”, y “convertible”) de un año al siguiente. Para estudiar este proceso se han utilizado los datos del histórico de los últimos 5 años para 5000 clientes, que están disponibles en el fichero coches. Este sistema se puede modelizar según una \(CMTD\) con espacio de estados \(S = \{s, sw, c\}\).
Obtén la matriz de saltos para cada cliente con createSequenceMatrix(as.character(coches[i,]). Luego, obtén la matriz de saltos total sumando los saltos de todos los clientes. A partir de ella, estima la matriz de probabilidades de transición.
En primer lugar, leemos los datos utilizando el enlace proporcionado, y verificamos la lectura, así como la posible presencia de valores faltantes. No se encuentra ninguna irregularidad al respecto de los posibles valores y no se detectan valores faltantes: todos los clientes tienen un registro completo de 5 años.
# lectura de datos
coches = read_csv("https://raw.githubusercontent.com/UMH1477/data/main/coches.csv")## Rows: 5000 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): V1, V2, V3, V4, V5
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# verificación de lectura
str(coches)## spec_tbl_df [5,000 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ V1: chr [1:5000] "c" "s" "s" "s" ...
## $ V2: chr [1:5000] "s" "sw" "sw" "s" ...
## $ V3: chr [1:5000] "sw" "sw" "sw" "s" ...
## $ V4: chr [1:5000] "s" "sw" "sw" "c" ...
## $ V5: chr [1:5000] "s" "s" "s" "c" ...
## - attr(*, "spec")=
## .. cols(
## .. V1 = col_character(),
## .. V2 = col_character(),
## .. V3 = col_character(),
## .. V4 = col_character(),
## .. V5 = col_character()
## .. )
## - attr(*, "problems")=<externalptr># valores faltantes
sum(is.na(coches))## [1] 0A continuación tratamos de obtener la matriz de saltos entre todos los estados posibles, contando los saltos (cambios de coche) de cada cliente registrado en cada una de las filas de la base de datos. Con la función createSequenceMatrix podemos contar los cambios de coche que cada cliente ha realizado a lo largo de los 5 años de registro, siempre y cuando reconozca como caracteres los valores que encuentra en las filas. Para ello habremos de forzar dicho reconocimiento con la función as.character. Para conseguir el número total de saltos relativos a los cambios de coche habremos de sumar todos los saltos contabilizados para cada uno de los clientes. Planteamos pues, un bucle que recorre las filas de la base de datos, contabilizando en cada fila, el número de cambios para todos los tipos de cambios posibles. Al final, para convertir estos conteos en probabilidades de transición, calculamos las sumas por filas y dividimos las filas por dichas sumas.
# Calculamos para todos los clientes el número total de transiciones
# especificamos todos los posibles estados
estados = c("s", "sw", "c")
# inicializamos a cero la matriz en la que almacenamos el número de saltos
cambios = matrix(rep(0, 9), ncol = 3, nrow = 3)
for(i in 1:nrow(coches)){
cambios = cambios + createSequenceMatrix(as.character(coches[i,]),
possibleStates = estados) }
# convertimos los conteos a probabilidades condicionadas (por filas)
p=cambios/apply(cambios,1,sum);p## c s sw
## c 0.25131349 0.7486865 0.0000000
## s 0.08606268 0.5520508 0.3618865
## sw 0.00000000 0.3960828 0.6039172Otra opción para conseguir los conteos es utilizar una función que hemos programado nosotros, conteo.saltos, en el Anexo.
cambios.bis = matrix(rep(0, 9), ncol = 3, nrow = 3)
for(i in 1:nrow(coches)){
cambios.bis = cambios.bis + conteo.saltos(coches[i,],estados) }
cambios.bis## s sw c
## s 5478 3591 854
## sw 3539 5396 0
## c 855 0 287Como última alternativa, podemos calcular directamente la matriz de transición utilizando la función markovchainFit, que trabaja de modo idéntico, reconociendo cada fila como un registro diferente (y sin necesidad de bucles).
p.bis=markovchainFit(coches)$estimate;
p.bis## Unnamed Markov chain
## A 3 - dimensional discrete Markov Chain defined by the following states:
## c, s, sw
## The transition matrix (by rows) is defined as follows:
## c s sw
## c 0.25131349 0.7486865 0.0000000
## s 0.08606268 0.5520508 0.3618865
## sw 0.00000000 0.3960828 0.6039172En cualquier caso, se muestra la matriz de transición obtenida, y apreciamos que el orden de filas y columnas no es el orden en el que se presentó el espacio de estados en el enunciado, \(S=\{s,sw,c\}\). Este hecho exigirá que tengamos especial cuidado con el orden en que combinemos información en adelante, por lo que recomendamos definir las etiquetas de los estados en ese orden.
¿Qué características tiene este proceso? ¿Es irreducible? ¿Tiene estados recurrentes? ¿Y transitorios?
Para dar respuesta a esta pregunta definimos el proceso CMTD y lo describimos con la función summary.
# Aptdo 2.Propiedades
# definimos las etiquetas de los estados en el orden en que aparecen en la matriz de transición
estados = c("c","s", "sw")
# definimos el proceso
proceso= new("markovchain", states = estados,byrow = TRUE, transitionMatrix = p,name="proceso")
# verificamos que se ha definido correctamente
proceso## proceso
## A 3 - dimensional discrete Markov Chain defined by the following states:
## c, s, sw
## The transition matrix (by rows) is defined as follows:
## c s sw
## c 0.25131349 0.7486865 0.0000000
## s 0.08606268 0.5520508 0.3618865
## sw 0.00000000 0.3960828 0.6039172# de modo alternativo podemos utilizar la matriz de transición obtenida del ajuste, para lo que hemos de seleccionar las filas del proceso para que entienda la matriz de transición
proceso.bis= new("markovchain", states = estados,byrow = TRUE, transitionMatrix = p.bis[1:3],name="proceso")
proceso.bis## proceso
## A 3 - dimensional discrete Markov Chain defined by the following states:
## c, s, sw
## The transition matrix (by rows) is defined as follows:
## c s sw
## c 0.25131349 0.7486865 0.0000000
## s 0.08606268 0.5520508 0.3618865
## sw 0.00000000 0.3960828 0.6039172# extraemos información del proceso
summary(proceso)## proceso Markov chain that is composed by:
## Closed classes:
## c s sw
## Recurrent classes:
## {c,s,sw}
## Transient classes:
## NONE
## The Markov chain is irreducible
## The absorbing states are: NONEUna vez verificado que el proceso se ha definido conforme a la matriz de transición que calculamos anteriormente, resumimos la información que hemos conseguido:
¿Cuál es la probabilidad de que un cliente que hoy contrata un “sedan”, el próximo año cambie de tipo de coche? Si este año 200 clientes contrataron un “sedan”, ¿cuántos de ellos mantendrán el mismo tipo de vehículo al próximo año?
Puesto que nos preguntan por probabilidades de cambio de coche en un año, acudimos a la matriz de transición (de un paso). Si sólo nos interesan los cambios de coche de los clientes que inicialmente han contratado un “sedan”, nos fijamos exclusivamente en las propiedades de transición para la fila “s”. Mirando estas probabilidades, este tipo de cliente cambiará de coche cuando no vuelva a contratar “sedan”, por lo que habremos de acumular las probabilidades para “c” y “sw”.
# Matriz de transición: predice cambios en un año
p## c s sw
## c 0.25131349 0.7486865 0.0000000
## s 0.08606268 0.5520508 0.3618865
## sw 0.00000000 0.3960828 0.6039172# los que inicialmente contratan un sedan son
p[2,]## c s sw
## 0.08606268 0.55205079 0.36188653# la probabilidad de que cambie de coche
prob.cambio=sum(p[2,c(1,3)]);prob.cambio## [1] 0.4479492La probabilidad de que un cliente que ha contratado un “sedan” cambie de vehículo al próximo año es pues de 0.45. Si este año 200 clientes contrataron un “sedan”, para saber cuántos conservarán este tipo de vehículo habremos de multiplicar dicha cifra por la probabilidad de que un cliente no cambie su “sedan” (0.55). El resultado es de 110 clientes que se espera que mantengan el “sedan”.
# De los 200 sedan mantienen coche:
round(200*p[2,2])## [1] 110¿Cuál es la probabilidad de que un cliente, al cabo de 3 años, tenga el mismo tipo de vehículo? ¿La fidelidad a un tipo de vehículo es similar para todos los tipos de vehículo?
Para predecir el comportamiento de los clientes al cabo de los tres años hemos de acudir a la matriz de transición de 3 pasos. Puesto que nos preguntan por la probabilidad de mantener el mismo tipo de vehículo, habremos de fijarnos en la diagonal de dicha matriz.
n=3
# calculamos la matriz de transición de n pasos
pn=ptran.n(p,n);pn## c s sw
## c 0.08382965 0.5348821 0.3812882
## s 0.06148553 0.5004002 0.4381143
## sw 0.04797134 0.4795137 0.4725150# alternativa para la matriz de transición
proceso^3## proceso^3
## A 3 - dimensional discrete Markov Chain defined by the following states:
## c, s, sw
## The transition matrix (by rows) is defined as follows:
## c s sw
## c 0.08382965 0.5348821 0.3812882
## s 0.06148553 0.5004002 0.4381143
## sw 0.04797134 0.4795137 0.4725150# la diagonal da las probabilidades de fidelidad
fidel=diag(pn);fidel## c s sw
## 0.08382965 0.50040016 0.47251495Obtenemos así que los vehículos para los que los clientes son más fieles (a tres años) son los de tipo “sedan”, con una probabilidad de permanencia de 0.5, seguidos de los “state wagon”, con una probabilidad de permanencia de 0.47, y finalmente los “convertibles”, cuya probabilidad de permanencia a tres años es sólo de 0.08. Claramente la fidelidad no es similar para todos los tipos de vehículos.
En función de lo que un cliente contrata la primera vez, ¿cuánto tiempo tarda en cambiar de tipo de vehículo? ¿En cada situación, cuál es el preferido y cuántos años tarda en cambiar a ese tipo?
Puesto que nos preguntan por el primer cambio de vehículo, habremos de acudir a los tiempos esperados de primer paso a otros tipos de vehículo distintos al que contrató inicialmente el cliente. El tiempo a cambio será pues el menor de los tiempos a cambio a cualquiera de los dos tipos de vehículo distintos al vehículo de partida.
# Tiempo de primer paso
meanFirstPassageTime(proceso)## c s sw
## c 0.00000 1.335673 4.416615
## s 22.23570 0.000000 3.080942
## sw 24.76042 2.524724 0.000000Tenemos así que los que contratan inicialmente un “convertible” tardan en media 1.33 años en cambiar de vehículo, y lo hacen a un “sedan”. Los que contratan inicialmente un “sedan” en 3.1 años han cambiado a un “state wagon”. Y finalmente, los que contratan un “state wagon” tardan 2.5 años en cambiar a un “sedan”.
Si a día de hoy la flota de coches contratados es del 40% sedan, 35% state wagon y 25% convertible, en función de las preferencias y migraciones de los clientes, ¿cuál será la composición de la flota el año próximo? ¿Y dentro de 5 años?
Para responder la cuestión sobre la composición de la flota el año próximo, conociendo la flota de inicio (distribución inicial), calculamos la distribución marginal, multiplicando las probabilidades condicionales (de un paso) por la distribución inicial (teorema de probabilidad total). Si queremos el resultado en porcentaje, multiplicamos por 100.
\[Pr(X_1=i)=\sum_j Pr(X_1=i|X_0=j) \cdot Pr(X_0=j)\]
Para responder la cuestión similar relativa a 5 años habremos de utilizar la matriz de transición de 5 pasos y repetir el mismo procedimiento para calcular la distribución marginal a 5 años. \[Pr(X_5=i)=\sum_j Pr(X_5=i|X_0=j) \cdot Pr(X_0=j)\]
# comprobamos ordenación de p
# flota inicial en orden (c,s,sw)
p0=c(0.25,0.4,0.35)
# Distribución marginal dentro de un año
porc.marg1=round(p0%*%p*100,1); porc.marg1## c s sw
## [1,] 9.7 54.7 35.6# Distribución marginal dentro de 5 años
porc.marg5=round(p0 %*% ptran.n(p,5)*100,1);porc.marg5## c s sw
## [1,] 5.8 49.4 44.8Así tenemos que en un año, si nos adaptamos al comportamiento del cliente, la composición de la flota sería de un 9.7% de coches “convertible”, un 54.7% de “sedan” y un 35.6% de “state wagon”. En cinco años sin embargo, habría que reducir aún más la flota de convertibles, pues el cupo de clientes que van a optar por ellos es del 5.8%, mientras que los que elijan “sedan” llegará al 49.4% y los que elijan “state wagon” llegará a un 44.8%.
¿Cuáles son las previsiones de la empresa a 10 años respecto a la proporción de vehículos de cada tipo que tendrá contratada en su flota de vehículos?
Para responder preguntas sobre la composición de la flota a 10 años, siempre que la empresa se adapte al comportamiento del cliente, consideramos que podemos asumir haber llegado a un comportamiento estacionario. Así la proporción de vehículos nos la da la distribución estacionaria.
prob.est=steadyStates(proceso); prob.est## c s sw
## [1,] 0.05666504 0.4929471 0.4503879Calculamos esta distribución y ratificamos el comportamiento que ya apreciábamos con 5 años, en el que “convertible” tiende a desaparecer (o al menos minimizarse), con un cupo del 5.7% en nuestra flota, que está compuesta mayoritariamente por vehículos “sedan”, en un 49.3%, y “state wagon”, con un cupo del 45%.
Si un “sedan” proporciona un beneficio anual de 1200 euros, el “station wagon” de 1500, y el “convertible” de 2500 euros. Por cada tipo de vehículo, ¿cuál es el beneficio esperado acumulado a finales del año próximo? Si este año arranco con una flota de 250 sedan, 400 state wagon y 150 convertibles, ¿cuáles serán los beneficios?
Puesto que nos preguntan por beneficios esperados acumulados (por cada coche) a finales del año próximo, habremos de contabilizar todo lo que se puede ganar durante este año, y lo que se pueda ganar durante el siguiente año, con lo cual estaremos contabilizando el coste de contrato de este año, y el de uno más. Hemos de considerar pues, el tiempo de uso (u ocupación) de estos dos primeros años, que nos viene dado por la matriz de ocupación de un paso. Este tiempo, multiplicado por los beneficios aplicables a cada tipo de vehículo, nos dará el beneficio esperado por vehículo.
# matriz de ocupación de un paso
mocupa=mocupa.proceso(proceso,1);mocupa## c s sw
## c 1.25131349 0.7486865 0.0000000
## s 0.08606268 1.5520508 0.3618865
## sw 0.00000000 0.3960828 1.6039172benef=c(2500,1200,1500)
benef.m=matrix(benef,ncol=1)
# beneficios por cada tipo de vehículo
g=mocupa %*% benef;g## [,1]
## c 4026.708
## s 2620.447
## sw 2881.175# También lo podemos calcular directamente con
expectedRewards(proceso,1,benef)## [1] 4026.708 2620.447 2881.175Tenemos así que el beneficio esperado por un “convertible” durante este año y el próximo es de 4026.71€, por un “sedan” de 2620.45€ y por un “state wagon” de 2881.18€.
Si este año arranco con una flota de 250 sedan, 400 state wagon y 150 convertibles, para calcular el beneficio total de mi flota hasta finales del año próximo simplemente habrá que multiplicar estos beneficios (por vehículo) por el número de vehículos de cada tipo, y sumar todas las cantidades.
#flota inicial
f.ini=c(150,250,400)
beneficio.total=f.ini %*% gSe obtendrán pues unos beneficios totales de 2.4115881^{6}€.
¿Y cuál será el beneficio total esperado por vehículo que se conseguirá dentro de 10 años? Si para entonces la flota total de vehículos es de 1000, ¿cuáles serán los beneficios esperados?
Cuando nos pregunta por el beneficio esperado por vehículo a 10 años vista, pensamos en evaluar beneficios anuales por vehículo en el largo plazo. Esto lo evaluamos con el beneficio esperado que predecimos a través de la distribución estacionaria, que es la que establece con qué probabilidad un cliente va a contratar cada uno de los tipos de vehículos disponibles. Si para entonces nuestra flota de vehículos es de 1000, bastará multiplicar este beneficio esperado por el número de vehículos.
# calculamos beneficios esperados por vehículo y año
g=steadyStates(proceso)%*%benef;g## [,1]
## [1,] 1408.781# beneficios esperados para la flota completa durante un año
g1000=g*1000; g1000## [,1]
## [1,] 1408781Así, por cada vehículo contratado el beneficio esperado durante un año (transcurridos ya 10) será de 1408.78€, cantidad que se multiplica por 1000 si nuestra flota es de 1000 vehículos, y rebasando así el millón de euros: 1.408.780€.
tiempo.pp(proceso, estado).ptran.n(ptran,n)mocupa.proceso(proceso, n)###############################
# Función para calcular los tiempos de primer paso
tiempo.pp <- function(proceso, estado)
{
# estados del proceso
estados <- states(proceso)
numestados <- length(estados)
# posición de los estados deseados
lestat <- length(estado)
pos <- which(estados %in% estado)
# matriz P_N
P_N <- proceso[-pos,-pos]
# vector de unos
vector.1 <- matrix(rep(1, numestados-lestat), ncol=1)
# sistema de ecuaciones
sistema <- diag(numestados-lestat) - P_N
# solución del sistema
solucion <- solve(sistema, vector.1)
return(solucion)
}
###############################
# Matriz de probs. transición de n pasos
ptran.n=function(ptran,n){
# ptran es la matriz de transición de 1 paso
# n son los pasos a dar
i=1
p=ptran
while(i<n){
p=p%*%ptran
i=i+1
}
return(p)
}
###############################
# Función para calcular los tiempos de ocupación
mocupa.proceso <- function(proceso, n)
{
# Número de estados del proceso
nestat <- dim(proceso)
# Estados
nombres<- names(proceso)
# Generamos la matriz de ocupaciones
mocupa <- diag(nestat)
dimnames(mocupa) <- list(nombres, nombres)
# mocupa <- matrix(rep(0, nestat*nestat),
# nrow = nestat, dimnames = list(nombres, nombres))
# Bucle de cálculo de los tiempos de ocupación
P=proceso[1:nestat,1:nestat]
for (i in 1:n)
mocupa <- mocupa + ptran.n(P,i)
return(mocupa)
}createSequenceMatrix para contar los saltos en una secuencia conteo.saltos(secuencia,estados).###############################
# Función que reproduce createSequenceMatrix para contar los saltos en una secuencia
conteo.saltos=function(secuencia,estados){
# cuenta transiciones entre estados a partir de una secuencia
n=length(secuencia)
nestados=length(estados)
pmat=matrix(0,ncol=nestados,nrow=nestados,dimnames=list(estados,estados))
for(i in 1:(n-1)){
for(j in 1:nestados){
for(k in 1:nestados){
pmat[j,k]=pmat[j,k]+1*(secuencia[i]==estados[j]& secuencia[i+1]==estados[k])
}}}
return(pmat)
}
###############################