Bagaimana cara memastikan properti dari matriks kovarians saat memasang model normal multivariat menggunakan kemungkinan maksimum?

22

Misalkan saya memiliki model berikut

yi=f(xi,θ)+εi

di mana ysayaRK , xsaya adalah vektor dari variabel penjelas, θ adalah parameter fungsi non-linear f dan εsayaN(0,Σ) , di mana Σ alami adalah K×K matriks.

Tujuannya adalah untuk memperkirakan θ dan Σ . Pilihan yang jelas adalah metode kemungkinan maksimum. Log-kemungkinan untuk model ini (dengan asumsi kita memiliki sampel (ysaya,xsaya),saya=1,...,n ) terlihat seperti

l(θ,Σ)=-n2log(2π)-n2logdetΣ-saya=1n(ysaya-f(xsaya,θ))Σ-1(y-f(xsaya,θ)))

Sekarang ini tampak sederhana, kemungkinan log ditentukan, dimasukkan ke dalam data, dan menggunakan beberapa algoritma untuk optimasi non-linear. Masalahnya adalah bagaimana memastikan bahwa pasti positif. Menggunakan misalnya dalam R (atau algoritma optimasi non-linear lainnya) tidak akan menjamin saya bahwa pasti positif.ΣΣoptimΣ

Jadi pertanyaannya adalah bagaimana memastikan bahwa tetap positif pasti? Saya melihat dua solusi yang mungkin:Σ

  1. Reparametrise sebagai mana adalah matriks segitiga-atas atau simetris. Maka akan selalu positif-pasti dan dapat tidak dibatasi.R R R Σ RΣRRRΣR

  2. Gunakan kemungkinan profil. Turunkan rumus untuk dan \ hat {\ Sigma} (\ theta) . Mulailah dengan beberapa \ theta_0 dan beralih \ hat {\ Sigma} _j = \ hat \ Sigma (\ hat \ theta_ {j-1}) , \ hat {\ theta} _j = \ hat \ theta (\ hat \ Sigma_ {j -1}) hingga konvergensi. Σ (θ)θ0 Σ j= Σ ( θ j-1) θ j= θ ( Σ j-1)θ^(Σ)Σ^(θ)θ0Σ^j=Σ^(θ^j-1)θ^j=θ^(Σ^j-1)

Apakah ada cara lain dan bagaimana dengan 2 pendekatan ini, apakah akan berhasil, apakah itu standar? Ini sepertinya masalah standar, tetapi pencarian cepat tidak memberi saya petunjuk. Saya tahu bahwa perkiraan Bayesian juga mungkin, tetapi untuk saat ini saya tidak ingin terlibat di dalamnya.

mpiktas
sumber
Saya memiliki masalah yang sama dalam algoritma Kalman, tetapi masalahnya jauh lebih rumit dan tidak semudah menggunakan trik Hamilton. Lalu saya bertanya-tanya apakah hal yang lebih sederhana untuk dilakukan adalah menggunakan . Dengan cara ini saya memaksakan kode untuk tidak memberikan kesalahan dan tidak mengubah solusi. Ini juga menguntungkan memaksa istilah ini memiliki tanda yang sama dengan bagian akhir dari kemungkinan. Ada ide? log(detΣ+1)
econ_pipo

Jawaban:

6

Dengan asumsi bahwa dalam membangun matriks kovarians, Anda secara otomatis menangani masalah simetri, kemungkinan log Anda akan ketika tidak pasti positif karena istilah dalam model yang tepat? Untuk mencegah kesalahan numerik jika Saya akan menghitung ulang dan, jika tidak positif, maka buat kemungkinan log sama dengan -Jika, jika tidak lanjutkan. Anda harus menghitung faktor penentu, jadi ini tidak dikenakan biaya perhitungan tambahan. Σ log d e t Σ d e t Σ < 0 d e t Σ-Σlogdet Σdet Σ<0det Σ

Makro
sumber
5

Ternyata Anda dapat menggunakan kemungkinan maksimum profil untuk memastikan properti yang diperlukan. Anda dapat membuktikan bahwa untuk diberikan θ , l ( θ , Σ ) dimaksimalkan olehθ^l(θ^,Σ)

Σ^=1nsaya=1nε^sayaε^saya,

dimana

ε^saya=ysaya-f(xsaya,θ^)

Maka dimungkinkan untuk menunjukkan itu

saya=1n(ysaya-f(xsaya,θ^))Σ^-1(y-f(xsaya,θ^)))=cHainst,

maka kita hanya perlu memaksimalkan

lR(θ,Σ)=-n2logdetΣ^.

Secara alami dalam hal ini akan memenuhi semua properti yang diperlukan. Buktinya identik untuk kasus ketika f adalah linier yang dapat ditemukan dalam Time Series Analysis oleh JD Hamilton halaman 295, maka saya menghilangkannya.Σf

mpiktas
sumber
3

Sebuah parameterisasi alternatif untuk matriks kovarians adalah dalam hal nilai eigen sudut p dan p ( p - 1 ) / 2 "Memberikan" θ i j .λ1,...,λhalhal(hal-1)/2θsayaj

Artinya, kita bisa menulis

Σ=GTΛG

di mana adalah ortonormal, danG

Λ=dsayaSebuahg(λ1,...,λhal)

dengan .λ1...λhal0

Sementara itu, dapat diparameterisasi unik dalam hal p ( p - 1 ) / 2 sudut, θ i j , di mana saya = 1 , 2 , . . . , P - 1 dan j = i , . . . , p - 1. [1]Ghal(hal-1)/2θsayajsaya=1,2,...,hal-1j=saya,...,hal-1

(detail yang akan ditambahkan)

[1]: Hoffman, Raffenetti, Ruedenberg. "Generalisasi Euler Angles ke Matriks Orthogonal N-Dimensi". J. Math. Phys 13, 528 (1972)

charles.y.zheng
sumber
Matriks sebenarnya ortogonal, karena Σ adalah matriks simetris. Ini adalah pendekatan yang akan saya rekomendasikan - Pada dasarnya sama dengan memutar vektor y i dan fungsi model f ( x i , θ ) sehingga kesalahannya independen, kemudian menerapkan OLS ke masing-masing komponen yang diputar (saya pikir). GΣysayaf(xsaya,θ)
probabilityislogic
2

Di sepanjang garis solusi charles.y.zheng, Anda mungkin ingin memodelkan , di mana Λ adalah matriks diagonal, dan C adalah faktorisasi Cholesky dari pembaruan peringkat ke Λ . Anda hanya perlu menjaga diagonal Λ positif untuk menjaga Σ pasti positif. Artinya, Anda harus memperkirakan diagonal Λ dan elemen C alih-alih memperkirakan Σ .Σ=Λ+CCΛCΛΛΣΛCΣ

shabbychef
sumber
Dapatkah elemen diagonal di bawah dalam pengaturan ini menjadi apa pun yang saya inginkan selama diagonal positif? Ketika mensimulasikan matriks dengan cara ini di numpy tidak semuanya pasti positif.
sztal
Λ