Modelo de Markov semioculto con parámetros de las probabilidades de emisión que dependen de los regresores

Question

Estoy intentando implantar un Modelo de Markov oculto para detectar agotamientos pasados en las ventas (una ruptura de existencias se produce cuando el minorista se queda sin un producto). Probablemente sea mejor decir un Semi - Modelo de Markov oculto déjame explicártelo. Supongamos que tengo un año de datos horarios sobre las ventas de un producto en un minorista.

VARIABLES DEL MODELO:

• $Y_t \in Z^+$ número de artículos vendidos en el momento t.
• Estados latentes $X_t \in \big\{0,1\big\}$ con $0$ que indica la indisponibilidad del producto por falta de existencias, $1$ su disponibilidad.

Así que tengo $n=2$ estados latentes y tienen $m_1 \in Z^+$ por lo que a $\infty$ número de símbolos de observación posibles para el estado 1 (todos los enteros positivos) , mientras que $m_0= 0 $ (el único símbolo observable es 0 cuando no me queda producto). En realidad podría elegir establecer los símbolos observables para el estado $1 \in \big\{sale / nosale\big\}$ para simplificar el modelo y tener un número finito de símbolos de observación para el estado $1$ .

PARÁMETROS DEL MODELO:

• la probabilidad inicial $P(X_0)$ distribución, (puedo fijar el inicio de la serie con una venta, de modo que $P(X_0=1)=1$ . )
• la matriz de probabilidad de transición de estado $P(X_{t+1}=j|X_t=i), i,j \in \big\{0,1\big\}$ ( Aquí considero un HMM de primer orden, pero no es imprescindible. )
• la emisión probabilidad $P(Y_t=m|X_t=i), m \in Z^+,i,j \in \big\{0,1\big\}$ .

Además, sé que $P(Y_t=m|X_t=0)=0 $
$\forall m \ne 0$ , ya que no puedo observar las ventas de un producto no disponible y eso podría ser una restricción útil añadida en la secuencia de estados que deseo encontrar, por ejemplo:

$Y_t=0,Y_{t+1}=3,Y_{t+2}=1,Y_{t+3}=0,\dots$ me hace imponer: $X_t \in \big\{0,1\big\}, X_{t+1}=1, X_{t+2}=1, X_{t+3} \in \big\{0,1\big\},\dots$

TENGO TAMBIÉN UN CONJUNTO DE REGRESORES (hora del día, día de la semana, promociones y precios de los otros productos,etc.) y estaba pensando en una distribución lineal/poisson dado el conjunto de regresores (piense en esto: hay una gran promoción para otro producto y obtengo 0 ventas para el producto de interés aunque esté disponible, por lo que me gustaría utilizar esta información adicional) para $P(Y_t=m|X_t=1) = Z_t*\beta $ donde $Z_t$ es el conjunto de regresores en el momento $t$ y $\beta$ son las estimaciones de los coeficientes de regresión OLS. Es una idea, también podría decir simplemente eso: $Y_t|X_t=1 \sim Poi(\lambda)$ donde $\lambda$ Deseo estimar.

Pido ayuda para seguir construyendo el modelo y algún consejo sobre cómo implementarlo en R/Python Conozco la función depmix del paquete depmixS4, pero no admite diferentes distribuciones de probabilidad de las observaciones (una forma de solucionarlo era considerar los datos diarios, no los horarios). En ese caso, dado que el minorista repone diariamente los artículos agotados, las probabilidades de emisión cuando el estado latente es $0$ (que en el caso diario representaría el evento " Se ha producido una ruptura de existencias ese día/el artículo no está disponible ese día "

$P(Y_t=m|X_t=0) >= 0 $ porque podría tener que el minorista se quedó sin un producto en algún momento durante un día después de vender algunos artículos del mismo.
En ese caso podría utilizar una Poisson para ambas probabilidades de observación, pero una será cero-inflada, mientras que la otra no.

No tengo conocimientos profundos sobre HMM, aunque ya he trabajado con Cadenas de Markov y Distribuciones Condicionales/Variables Omitidas en Estadística Bayesiana y Aplicada. Mi principal referencia era hmm_fundamentals

EDITAR :

Quiero hacer hincapié en el hecho de que me gustaría explotar la información relativa a los regresores, de modo que, cuando el Estado Oculto es $X_t=1$ = producto disponible , la distribución de la Observación Probabilidad $f(Y_t|X_t=1)$ sigue a $Poi(\lambda_t)$ donde $\lambda_t= Z_t*\beta_{OLS}$ con $Z_t$ que indica el conjunto de regresores en el momento $t$ .

EDITAR 2 :

He observado en los datos que algunos productos que suelen venderse mucho no tienen un claro agotamiento de existencias, permítanme explicarlo:

Y a continuación, la semana del 11 al 18 de octubre:

Como se puede ver hay un periodo bastante largo en el que el producto rara vez se vende y además cuando se vende es muy por debajo de su comportamiento habitual (he utilizado este periodo de 3 meses, pero su comportamiento habitual se confirma en un periodo más largo). Yo interpreté esto como: quedan pocos artículos (quizás antiguos artículos almacenados) y pasa el tiempo antes de otra reposición. En este caso el problema cambiaría, ya que:

$P(Y_t>0|X_t=0) \geq 0$ y el problema sería puramente de Modelo de Markov Oculto (sin estados latentes conocidos en la secuencia).

Answer 1

Aquí tienes un código R rápido para empezar. Advertencias importantes: Lo escribí yo mismo, así que podría tener errores, ser estadísticamente incorrecto, estar mal diseñado... ¡úsalo bajo tu propia responsabilidad!

params_are_valid <- function(params) {
    stopifnot("lambdas" %in% names(params))   # Poisson parameters for each hidden state
    stopifnot(all(params$lambas >= 0.0))
    stopifnot(length(params$lambdas) == params$n_hidden_states)
    stopifnot("mu" %in% names(params))  # Initial distribution over hidden states
    stopifnot(length(params$mu) == params$n_hidden_states)
    stopifnot(isTRUE(all.equal(sum(params$mu), 1.0)))  # Probabilities sum to 1
    stopifnot(all(params$mu >= 0.0))
    stopifnot("P" %in% names(params))  # Transition probabilities for hidden state
    stopifnot(nrow(params$P) == params$n_hidden_states && ncol(params$P) == params$n_hidden_states)
    stopifnot(isTRUE(all.equal(rowSums(params$P), rep(1, params$n_hidden_states))))  # Probabilities sum to 1
    stopifnot(all(params$P >= 0))
    return(TRUE)
}

baum_welch_poisson <- function(observed_y, params) {
    ## Baum-Welch algorithm for HMM with discrete hidden x
    ## Discrete observations y (NAs allowed) with y|x ~ Poisson(params$lambdas[x])
    ## Written following Ramon van Handel's HMM notes, page 40, algorithm 3.2
    ## https://www.princeton.edu/~rvan/orf557/hmm080728.pdf
    ## Careful, his observation index is in {0, 1, ... , n} while I use {1, 2, ... , length(observed_y)}
    y_length <- length(observed_y)
    stopifnot(y_length > 1)
    stopifnot(all(observed_y >= 0 | is.na(observed_y)))
    stopifnot(all(observed_y %% 1 == 0 | is.na(observed_y)))  # Integer observations
    stopifnot(params_are_valid(params))
    c <- vector("numeric", y_length)
    if(is.na(observed_y[1])) {
        upsilon <- rep(1, params$n_hidden_states)
    } else {
        upsilon <- dpois(observed_y[1], lambda=params$lambdas)
    }
    stopifnot(length(upsilon) == params$n_hidden_states)
    c[1] <- sum(upsilon * params$mu)
    ## Matrix pi_contemporaneous gives probabilities over x_k conditional on {y_1, y_2, ... , y_k}
    ## Notation in van Handel's HMM notes is pi_k, whereas pi_{k|n} conditions on full history of y
    pi_contemporaneous <- matrix(NA, params$n_hidden_states, y_length)
    pi_contemporaneous[, 1] <- upsilon * params$mu / c[1]
    upsilon_list <- list()
    upsilon_list[[1]] <- upsilon
    for(k in seq(2, y_length)) {
        ## Forward loop
        if(is.na(observed_y[k])) {
            upsilon <- rep(1, params$n_hidden_states)
        } else {
            upsilon <- dpois(observed_y[k], lambda=params$lambdas)
        }
        upsilon_list[[k]] <- upsilon  # Cache for backward loop
        pi_tilde <- upsilon * t(params$P) %*% pi_contemporaneous[, k-1]
        c[k] <- sum(pi_tilde)
        pi_contemporaneous[, k] <- pi_tilde  / c[k]
    }
    beta <- matrix(NA, params$n_hidden_states, y_length)
    beta[, y_length] <- 1 / c[y_length]
    ## Matrix pi gives probabilities over x conditional on full history of y
    ## Notation in van Handel's HMM notes is pi_{k|n}, as opposed to pi_k
    pi <- matrix(NA, params$n_hidden_states, y_length)
    pi[, y_length] <- pi_contemporaneous[, y_length]
    pi_transition_list <- list()  # List of posterior probabilities over hidden x transitions
    for(k in seq(1, y_length - 1)) {
        ## Backward loop
        upsilon <- diag(upsilon_list[[y_length - k + 1]], params$n_hidden_states, params$n_hidden_states)
        pi_matrix <- diag(pi_contemporaneous[, y_length - k],
                          params$n_hidden_states, params$n_hidden_states)
        beta_matrix <- diag(beta[, y_length - k + 1], params$n_hidden_states, params$n_hidden_states)
        beta[, y_length - k] <- params$P %*% upsilon %*% beta[, y_length - k + 1] / c[y_length - k]
        pi_transition_list[[y_length - k]] <- pi_matrix %*% params$P %*% upsilon %*% beta_matrix
        stopifnot(isTRUE(all.equal(sum(pi_transition_list[[y_length - k]]), 1.0)))
        pi[, y_length - k] <- rowSums(pi_transition_list[[y_length - k]])
    }
    loglik <- sum(log(c))
    return(list(loglik=loglik, pi=pi, pi_transition_list=pi_transition_list))
}

## Notice that params0$lambas is zero for first hidden state, i.e. first hidden state represents a stockout
params0 <- list("n_hidden_states"=2,
                "mu"=c(0.10, 0.90),  # Initial distribution over hidden states
                "P"=rbind(c(0.50, 0.50),
                          c(0.10, 0.90)),  # Transition probabilities for hidden state
                "lambdas"=c(0, 2))  # Observe y|x ~ Poisson(lambdas[x])

observations <- c(1, 2, 3, NA, NA, 0, 0, 0, 5, 5, 0, 5, 6, 7, NA, 5, 8, 9, 0, 6)
bw_list <- baum_welch_poisson(observations, params0)

round(bw_list$pi, 3)  # Posterior probabilities over hidden states given observations (first row is stockout state)
bw_list$pi[, which(observations != 0)]  # Sanity check: we're certain we're in hidden state 2 whenever observations > 0
bw_list$pi[, which(observations == 0 |
                   is.na(observations))]  # When observations are zero (or NA), we could be in either state

El código ejecuta Baum-Welch en un HMM simple de dos estados con observaciones de Poisson. En este caso, establezco que la tasa de Poisson (lambda) sea cero para el primer estado oculto (y positiva para el segundo estado oculto), lo que hace que el primer estado oculto sea un estado de "agotamiento" como en su ejemplo.

En realidad, el ejemplo no ajusta los parámetros (los estima a partir de los datos), para lo que habría que escribir, por ejemplo, una función de maximización de expectativas (EM).