¿Cómo asegurar las propiedades de la matriz de covarianza cuando se ajusta el modelo normal multivariado usando la máxima probabilidad?

22

Supongamos que tengo el siguiente modelo

yi=f(xi,θ)+εi

donde yiRK , xi es un vector de variables explicativas, θ es los parámetros de la función no lineal f y εiN(0,Σ) , donde Σ naturalmente es K×K matriz.

El objetivo es el habitual para estimar θ y Σ . La elección obvia es el método de máxima verosimilitud. Diario de probabilidad para este modelo (suponiendo que tenemos una muestra (yi,xi),i=1,...,n ) se parece

l(θ,Σ)=n2log(2π)n2logdetΣi=1n(yif(xi,θ))Σ1(yf(xi,θ)))

Ahora, esto parece simple, se especifica la probabilidad de registro, se colocan datos y se usa algún algoritmo para la optimización no lineal. El problema es cómo garantizar que Σ sea ​​positivo definido. Usar, por ejemplo, optimen R (o cualquier otro algoritmo de optimización no lineal) no me garantizará que Σ sea ​​positivo definido.

Entonces, la pregunta es cómo garantizar que mantenga positivo definido. Veo dos posibles soluciones:Σ

  1. Reparametrise como R R donde R es una matriz triangular superior o simétrica. Entonces Σ siempre será positivo-definido y R puede estar sin restricciones.ΣRRRΣR

  2. Usa la probabilidad de perfil. Derivar las fórmulas para θ ( Σ ) y Σ ( θ ) . Start con algunos θ 0 y iterate Σ j = Σ ( θ j - 1 ) , θ j = θ ( Σ j - 1 ) hasta la convergencia.θ^(Σ)Σ^(θ)θ0 0Σ^j=Σ^(θ^j-1)θ^j=θ^(Σ^j-1)

¿Hay alguna otra manera y qué pasa con estos 2 enfoques, funcionarán, son estándar? Esto parece un problema bastante estándar, pero la búsqueda rápida no me dio ningún indicio. Sé que la estimación bayesiana también sería posible, pero por el momento no me gustaría participar en ella.

mpiktas
fuente
Tengo el mismo problema en un algoritmo de Kalman, pero el problema es mucho más complicado y no es tan fácil usar el truco de Hamilton. Me pregunto si una cosa más simple sería usar . De esta manera, fuerzo el código a no dar un error y no cambio la solución. Esto también tiene el beneficio de obligar a este término a tener el mismo signo que la parte final de la probabilidad. ¿Algunas ideas? log(detΣ+1)
econ_pipo

Respuestas:

6

Suponiendo que al construir la matriz de covarianza, usted se está ocupando automáticamente del problema de simetría, su probabilidad de registro será cuando Σ no es definitivo positivo debido al término log d e t Σ en el modelo, ¿verdad? Para evitar un error numérico si d e t Σ < 0 , precalcularía d e t Σ y, si no es positivo, haría que la probabilidad de registro sea igual -Inf, de lo contrario, continúe. De todos modos, debe calcular el determinante, por lo que esto no le está costando ningún cálculo adicional. Σlogdet Σdet Σ<0det Σ

Macro
fuente
5

Como resultado, puede usar la máxima probabilidad de perfil para garantizar las propiedades necesarias. Se puede demostrar que para determinado θ , l ( θ , Σ ) se maximizaθ^l(θ^,Σ)

Σ^=1ni=1nε^iε^i,

dónde

ε^i=yif(xi,θ^)

Entonces es posible demostrar que

i=1n(yif(xi,θ^))Σ^1(yf(xi,θ^)))=const,

por lo tanto solo necesitamos maximizar

lR(θ,Σ)=n2logdetΣ^.

Naturally in this case Σ will satisfy all the necessary properties. The proofs are identical for the case when f is linear which can be found in Time Series Analysis by J. D. Hamilton page 295, hence I omitted them.

mpiktas
fuente
3

An alternative parameterization for the covariance matrix is in terms of eigenvalues λ1,...,λp and p(p1)/2 "Givens" angles θij.

That is, we can write

Σ=GTΛG

where G is orthonormal, and

Λ=diag(λ1,...,λp)

λ1...λp0

Gp(p1)/2 angles, θij, where i=1,2,...,p1 and j=i,...,p1.[1]

(details to be added)

[1]: Hoffman, Raffenetti, Ruedenberg. "Generalization of Euler Angles to N‐Dimensional Orthogonal Matrices". J. Math. Phys. 13, 528 (1972)

charles.y.zheng
fuente
The matrix G is actually orthogonal, because Σ is a symmetric matrix. This is the approach I was going to recommend - Basically amounts to rotating the yi vector and the model function f(xi,θ) so that the errors are independent, then applying OLS to each of the rotated components (I think).
probabilityislogic
2

Along the lines of charles.y.zheng's solution, you may wish to model Σ=Λ+CC, where Λ is a diagonal matrix, and C is a Cholesky factorization of a rank update to Λ. You only then need to keep the diagonal of Λ positive to keep Σ positive definite. That is, you should estimate the diagonal of Λ and the elements of C instead of estimating Σ.

shabbychef
fuente
Can below diagonal elements in this settings be anything I want as long as the diagonal is positive? When simulate matrices this way in numpy not all of them are positive definite.
sztal
Λ is a diagonal matrix.
shabbychef