¿Por qué la probabilidad en el filtro Kalman se calcula utilizando los resultados del filtro en lugar de los resultados del suavizado?

Question

Estoy utilizando el filtro Kalman de una manera muy estándar. El sistema está representado por la ecuación de estado $x_{t+1}=Fx_{t}+v_{t+1}$ y la ecuación de observación $y_{t}=Hx_{t}+Az_{t}+w_{t}$ .

Los libros de texto enseñan que tras aplicar el filtro de Kalman y obtener las "previsiones de un paso adelante" $\hat{x}_{t|t-1}$ (o "estimación filtrada"), debemos utilizarlas para calcular la función de verosimilitud:

$f_{y_{t}|\mathcal{I}_{t-1},z_{t}}\left(y_{t}|\mathcal{I}_{t-1},z_{t}\right)=\det\left[2\pi\left(HP_{t|t-1}H^{\prime}+R\right)\right]^{-\frac{1}{2}}\exp\left\{ -\frac{1}{2}\left(y_{t}-H\hat{x}_{t|t-1}-Az_{t}\right)^{\prime}\left(HP_{t|t-1}H^{\prime}+R\right)^{-1}\left(y_{t}-H\hat{x}_{t|t-1}-Az_{t}\right)\right\}$

Mi pregunta es: ¿Por qué se calcula la función de verosimilitud utilizando la "estimación filtrada" $\hat{x}_{t|t-1}$ y no la "estimación suavizada" $\hat{x}_{t|T}$ ? ¿No es $\hat{x}_{t|T}$ una mejor estimación del vector de estado?

Answer 1

En general, por la regla del producto, la probabilidad exacta se puede escribir $$ f(y_1,\dots,y_n)=f(y_1)\prod_{i=2}^n f(y_i|y_1,\dots,y_{i-1}). $$ De la suposición del modelo de espacio de estados, se deduce que el vector de expectativas y la matriz de varianzas de cada $y_i$ condicionada a las observaciones pasadas puede expresarse como \begin{align} E(y_i|y_1,\dots,y_{i-1}) &= E(Hx_{t}+Az_{t}+w_{t}|y_1,\dots,y_{i-1}) \\&= HE(x_{t}|y_1,\dots,y_{i-1})+Az_{t}+Ew_{t} \\&= H\hat x_{t|t-1}+Az_{t}, \end{align} y \begin{align} \mathrm{Var}(y_i|y_1,\dots,y_{i-1}) &= \mathrm{Var}(Hx_{t}+Az_{t}+w_{t}|y_1,\dots,y_{i-1}) \\&= H\mathrm{Var}(x_{t}|y_1,\dots,y_{i-1})H'+ \mathrm{Var}w_t \\&= HP_{t|t-1}H'+R. \end{align} Así que esto le da la probabilidad exacta sin calcular ninguna estimación suavizada.

Aunque, por supuesto, se podrían utilizar las estimaciones suavizadas, que efectivamente son mejores estimaciones de los estados desconocidos, esto no daría la función de probabilidad. En efecto, se utilizaría el valor observado de $y_i$ para estimar su propio valor esperado, por lo que parece probable que esto provoque algún sesgo en las estimaciones resultantes.

Answer 2

Para responder a tu pregunta: puedes utilizar la densidad de suavizado. Pero no es necesario. La respuesta de Jarle Tufto tiene la descomposición que estás utilizando. Pero hay otras.

Uso de las Recursiones de Kalman

Aquí estás evaluando la probabilidad como $$ f(y_1, \ldots, y_n) = f(y_1)\prod_{i=2}^nf(y_i|y_1, \ldots, y_{i-1}). $$

Sin embargo, las medias y las varianzas no siempre definen completamente las distribuciones de probabilidad en general. La siguiente es la descomposición que se utiliza para pasar de filtrar distribuciones $f(x_{i-1}|y_1,\ldots,y_{i-1})$ a las probabilidades condicionales $f(y_i|y_1,\ldots,y_{i-1})$ :

$$ f(y_i|y_1, \ldots, y_{i-1}) = \iint f(y_i|x_i)f(x_i|x_{i-1})f(x_{i-1}|y_1, \ldots, y_{i-1})dx_{i} dx_{i-1} \tag{1}. $$

Aquí $f(x_i|x_{i-1})$ es la densidad de transición de estado...parte del modelo, y $f(y_i|x_i)$ es la densidad de observación... parte del modelo de nuevo. En tu pregunta los escribes como $x_{t+1}=Fx_{t}+v_{t+1}$ y $y_{t}=Hx_{t}+Az_{t}+w_{t}$ respectivamente. Es lo mismo.

Cuando se obtiene la distribución de predicción de estado de un paso adelante, eso es computar $\int f(x_i|x_{i-1})f(x_{i-1}|y_1, \ldots, y_{i-1}) dx_{i-1}$ . Cuando integras de nuevo, obtienes (1) por completo. Si escribes esa densidad completamente en tu pregunta, es lo mismo.

Aquí sólo estás usando descomposiciones de distribuciones de probabilidad, y suposiciones sobre el modelo. Este cálculo de probabilidad es un cálculo exacto. No hay nada discrecional que puedas usar para hacerlo mejor o peor.

Uso del algoritmo EM

Que yo sepa, no hay otra forma de evaluar la probabilidad directamente en este tipo de modelo de espacio de estados. Sin embargo, se puede hacer una estimación de máxima verosimilitud evaluando una función diferente: se puede utilizar el algoritmo EM. En el paso de la expectativa (paso E) se calcularía $$ \int f(x_1, \ldots, x_n|y_1,\ldots y_n) \log f(y_1,\ldots,y_n,x_1, \ldots,x_n) dx_{1:n} = E_{smooth}[\log f(y_1,\ldots,y_n,x_1, \ldots,x_n)]. $$ Aquí $f(y_1,\ldots,y_n,x_1, \ldots,x_n)$ es la probabilidad de los "datos completos", y estás tomando la expectativa del logaritmo de eso con respecto a la densidad de suavización conjunta. Lo que suele ocurrir es que, al tomar el logaritmo de esta verosimilitud de los datos completos, los términos se dividen en sumas y, debido a la linealidad del operador de expectativas, se toman las expectativas con respecto a las distribuciones marginales de suavizado (las que mencionas en tu pregunta).

Otras cosas

He leído en algunos sitios que el EM es una forma "más estable" de maximizar la probabilidad, pero nunca he visto realmente este punto bien argumentado, ni he visto esta palabra "estable" definida en absoluto, pero tampoco he examinado realmente esto más a fondo. Ninguno de estos algoritmos evita el problema de los máximos locales/globales. Personalmente, tiendo a utilizar el Kalman más a menudo sólo por costumbre.

Es cierto que las estimaciones suavizadas del estado tienen una varianza más pequeña normalmente que el filtrado, así que supongo que tienes razón al tener alguna intuición sobre esto, pero no estás usando realmente los estados. La probabilidad que estás tratando de maximizar no es una función de los estados.

Answer 3

Creo que una mejor respuesta a "por qué" no se utiliza la distribución de suavizado (normalmente) es la eficiencia. En principio, es sencillo calcular la probabilidad marginal (de suavizado) en un sentido de dejar de lado lo siguiente. Elimine la observación j, ejecute el suavizador de Kalman en los datos restantes. A continuación, evalúe la probabilidad de la observación y(j) no vista. Repita esto para todas las j. Sume las probabilidades logarítmicas. Las versiones más rápidas de esto funcionan con bloques (aleatorios) de muestras retenidas (como k-fold CV). Obsérvese que este esquema requiere una implementación más general del filtro/suavizador de Kalman que puede omitir arbitrariamente las actualizaciones de las mediciones cuando sea necesario. El paso de retroceso/suavizado no accede a las mediciones (algoritmo RTS de todos modos) y sigue siendo el mismo.

Si la serie temporal es "lo suficientemente larga", es probable que no resulte muy útil hacer esto, ya que la probabilidad de filtrado "quema" su transitorio inicial. Pero si el conjunto de datos es corto, la probabilidad de suavizado más cara puede merecer la pena. Un alisador de retardo fijo podría ser una solución intermedia.