LogLikelihood Parameter Estimation for Linear Gaussian Kalman Filter

13

He escrito un código que puede hacer el filtrado de Kalman (utilizando varios filtros de tipo Kalman diferentes [Filtro de información et al.]) Para el Análisis de espacio de estado lineal gaussiano para un vector de estado n-dimensional. Los filtros funcionan muy bien y estoy obteniendo una buena salida. Sin embargo, la estimación de parámetros a través de la estimación de loglikelihood me confunde. No soy un estadístico sino un físico, así que por favor sea amable.

Consideremos el modelo lineal de espacio de estado gaussiano

yt=Ztαt+ϵt,
αt+1=Ttαt+Rtηt,

donde yt es nuestro vector de observación, αt nuestro vector de estado en el paso de tiempo t . Las cantidades en negrita son las matrices de transformación del modelo de espacio de estado que se establecen de acuerdo con las características del sistema en consideración. También tenemos

ϵtNID(0,Ht),
ηtNID(0,Qt),
α1NID(a1,P1).

donde . Ahora, he derivado e implementado la recursividad para el filtro de Kalman para este modelo de espacio de estado genérico al adivinar los parámetros iniciales y las matrices de varianza H 1 y Q 1. Puedo producir gráficos comot=1,,nH1Q1

Filtro Kalman

donde los puntos son los niveles de agua del río Nilo en enero durante más de 100 años, la línea es el estado estimado de Kalamn, y las líneas discontinuas son los niveles de confianza del 90%.

Ahora, para este conjunto de datos 1D, las matrices y Q t son simplemente escalares σ ϵ y σ η respectivamente. Así que ahora quiero obtener los parámetros correctos para estos escalares usando la salida del filtro de Kalman y la función de verosimilitudHtQtσϵση

logL(Yn)=np2log(2π)12t=1n(log|Ft|+vtTFt1vt)

Donde es el error de estado y F t es la varianza del error de estado. Ahora, aquí es donde estoy confundido. Del filtro de Kalman, tengo toda la información que necesito para resolver L , pero esto parece no acercarme más a poder calcular la probabilidad máxima de σ ϵ y σ η . Mi pregunta es ¿cómo puedo calcular la probabilidad máxima de σ ϵ y σ η usando el enfoque de loglikelihood y la ecuación anterior? Un desglose algorítmico sería como una cerveza fría para mí en este momento ...vtFtLσϵσησϵση

Gracias por tu tiempo.


Nota. Para el caso 1D, y H t = σ 2 η . Este es el modelo univariante de nivel local.Ht=σϵ2Ht=ση2

Caballero de la Luna
fuente

Respuestas:

11

Cuando ejecuta el filtro de Kalman como lo ha hecho, con valores dados de y σ 2 η , obtiene una secuencia de innovaciones ν t y sus covarianzas F t , por lo tanto, puede calcular el valor de log L ( Y n )σϵ2ση2νtFtlogL(Yn) usando La fórmula que le das.

En otras palabras, puede considerar el filtro de Kalman como una forma de calcular una función implícita de y σ 2 η . Lo único que debe hacer es empaquetar este cálculo en una función o subrutina y manejar esa función en una rutina de optimización, como en R. Esa función debe aceptar como entradas σ 2 ϵ y σ 2 η y devolver el registro L ( S n ) .σϵ2ση2optimσϵ2ση2logL(Yn)

Algunos paquetes en R (p dlm. Ej. ) Lo hacen por usted (vea, por ejemplo, la función dlmMLE).

F. Tusell
fuente
Gracias por su respuesta. Aprecio que parece que tengo todos los componentes necesarios para calcular la verosimilitud explícitamente, sin embargo, todas las referencias parecen sugerir que uso y σ η como incógnitas en la función de verosimilitud y maximizo esto usando un método de tipo Newton. Esto es lo que me confunde; "loglikelihood se maximiza numéricamente con respecto al vector de estado desconocido" - ¿cómo? σϵση
MoonKnight
El cálculo de la probabilidad no es tan explícito, ya que y σ η no aparecen explícitamente en la expresión de log L ( Y n ) . Por el contrario, influyen en la probabilidad a través de ν t y F t . Por lo tanto, debe ejecutar el filtro de Kalman para calcular el registro L ( Y n ) para cada par de valores de σ ϵ y σ η . Una vez que codifica eso en forma de una función, lo maneja a una función de maximización de tipo Newton (o cualquier propósito general) y eso es todo.σϵσηlogL(Yn)νtFtlogL(Yn)σϵση
F. Tusell
1
Resulta que tengo un código detallado (en R) que muestra cómo hacer esto precisamente para los datos del Nilo. Lo uso como ilustración para mis alumnos. Desafortunadamente está en español, pero espero que el código sea bastante claro (y puedo traducir los comentarios si no). Puede obtener este ejemplo de et.bs.ehu.es/~etptupaf/N4.html .
F. Tusell
Esto es enormemente útil. Muchas gracias por tu tiempo. Tu comentario ha ayudado mucho! A veces es difícil "ver la madera de los árboles" y tener algo simple explicado explícitamente es todo lo que se requiere ... Gracias de nuevo.
MoonKnight
También me gustaría preguntar si podría echar un vistazo a la página donde atraviesas la recursión de suavizado de estado. ¿Tu suavidad se ve mejor que la mía y no estoy seguro de por qué? He intentado encontrarlo en su sitio web pero no puedo encontrar la página requerida ...
MoonKnight