Saya ingin menerapkan algoritma EM manual dan kemudian membandingkannya dengan hasil normalmixEM
dari mixtools
paket. Tentu saja, saya akan senang jika keduanya menghasilkan hasil yang sama. Referensi utama adalah Geoffrey McLachlan (2000), Finite Mixture Models .
Saya memiliki kerapatan campuran dari dua Gaussians, dalam bentuk umum, log-likelihood diberikan oleh (McLachlan halaman 48):
Langkah E sekarang, perhitungan ekspektasi bersyarat:
Saya mencoba menulis kode R (data dapat ditemukan di sini ).
# EM algorithm manually
# dat is the data
# initial values
pi1 <- 0.5
pi2 <- 0.5
mu1 <- -0.01
mu2 <- 0.01
sigma1 <- 0.01
sigma2 <- 0.02
loglik[1] <- 0
loglik[2] <- sum(pi1*(log(pi1) + log(dnorm(dat,mu1,sigma1)))) +
sum(pi2*(log(pi2) + log(dnorm(dat,mu2,sigma2))))
tau1 <- 0
tau2 <- 0
k <- 1
# loop
while(abs(loglik[k+1]-loglik[k]) >= 0.00001) {
# E step
tau1 <- pi1*dnorm(dat,mean=mu1,sd=sigma1)/(pi1*dnorm(x,mean=mu1,sd=sigma1) +
pi2*dnorm(dat,mean=mu2,sd=sigma2))
tau2 <- pi2*dnorm(dat,mean=mu2,sd=sigma2)/(pi1*dnorm(x,mean=mu1,sd=sigma1) +
pi2*dnorm(dat,mean=mu2,sd=sigma2))
# M step
pi1 <- sum(tau1)/length(dat)
pi2 <- sum(tau2)/length(dat)
mu1 <- sum(tau1*x)/sum(tau1)
mu2 <- sum(tau2*x)/sum(tau2)
sigma1 <- sum(tau1*(x-mu1)^2)/sum(tau1)
sigma2 <- sum(tau2*(x-mu2)^2)/sum(tau2)
loglik[k] <- sum(tau1*(log(pi1) + log(dnorm(x,mu1,sigma1)))) +
sum(tau2*(log(pi2) + log(dnorm(x,mu2,sigma2))))
k <- k+1
}
# compare
library(mixtools)
gm <- normalmixEM(x, k=2, lambda=c(0.5,0.5), mu=c(-0.01,0.01), sigma=c(0.01,0.02))
gm$lambda
gm$mu
gm$sigma
gm$loglik
Algoritma tidak berfungsi, karena beberapa pengamatan memiliki kemungkinan nol dan log ini -Inf
. Di mana kesalahan saya?
sumber
Jawaban:
Anda memiliki beberapa masalah dalam kode sumber:
Seperti yang ditunjukkan oleh @Pat, Anda seharusnya tidak menggunakan log (dnorm ()) karena nilai ini dapat dengan mudah menuju infinity. Anda harus menggunakan logmvdnorm
Saat Anda menggunakan jumlah , berhati-hatilah untuk menghapus nilai tak terbatas atau hilang
Anda mengulang variabel k salah, Anda harus memperbarui loglik [k + 1] tetapi Anda memperbarui loglik [k]
Nilai awal untuk metode dan mixtools Anda berbeda. Anda menggunakan dalam metode Anda, tetapi menggunakan untuk mixtools (yaitu standar deviasi, dari manual mixtools).σΣ σ
Data Anda tidak terlihat seperti campuran normal (lihat histogram yang saya rencanakan di bagian akhir). Dan satu komponen campuran memiliki sd yang sangat kecil, jadi saya sewenang-wenang menambahkan baris untuk mengatur dan agar sama untuk beberapa sampel ekstrim. Saya menambahkannya hanya untuk memastikan kode dapat bekerja.τ 2τ1 τ2
Saya juga menyarankan Anda memasukkan kode lengkap (misalnya bagaimana Anda menginisialisasi loglik []) dalam kode sumber Anda dan membuat indentasi kode untuk membuatnya mudah dibaca.
Setelah semua, terima kasih telah memperkenalkan paket mixtools , dan saya berencana untuk menggunakannya dalam penelitian masa depan saya.
Saya juga menaruh kode kerja saya untuk referensi Anda:
Historgram
sumber
loklik <- rep(NA, 100)
yang akan mengalokasikan loglik [1], loglik [2] ... loglik [100]. Saya mengajukan pertanyaan itu karena dalam kode asli Anda, saya tidak menemukan penghapusan loglik, mungkin kode tersebut terpotong saat menempel?Saya terus mendapatkan kesalahan ketika mencoba membuka file .rar Anda, tapi itu mungkin saya melakukan sesuatu yang konyol.
Saya tidak dapat melihat kesalahan yang jelas dalam kode Anda. Kemungkinan alasan Anda mendapatkan nol adalah karena presisi titik apung. Ingat, ketika Anda menghitung , Anda mengevaluasi . Tidak perlu perbedaan yang sangat besar antara dan untuk ini dibulatkan menjadi 0 ketika Anda melakukannya di komputer. Ini terlihat dua kali lipat dalam model campuran, karena beberapa data Anda tidak akan "ditugaskan" untuk masing-masing komponen campuran sehingga dapat berakhir sangat jauh darinya. Secara teori, poin-poin ini juga harus berakhir dengan nilaiexp ( - 0,5 ( y - μ ) 2 / σ 2 ) μ y τf( y; θ ) exp( - 0,5 ( y- μ )2/ σ2) μ y τ ketika Anda mengevaluasi kemungkinan log, menangkal masalah - tetapi berkat kesalahan floating point, kuantitasnya telah dievaluasi sebagai -Jika pada tahap ini, jadi semuanya rusak :).
Jika itu masalahnya, ada beberapa solusi yang mungkin:
Pertama adalah memindahkan ke dalam logaritma. Jadi alih-alih mengevaluasiτ
evaluasi
Secara matematis sama, tetapi pikirkan tentang apa yang terjadi ketika dan adalah . Saat ini Anda mendapatkan:f( y| θ ) ≈ 0τ ≈ 0
tetapi dengan tau pindah Anda dapatkan
dengan asumsi R mengevaluasi (Saya tidak tahu apakah itu benar atau tidak karena saya cenderung menggunakan matlab)00= 1
Solusi lain adalah untuk memperluas hal-hal di dalam logaritma. Dengan asumsi Anda menggunakan logaritma natural:
Secara matematis sama, tetapi harus lebih tahan terhadap kesalahan floating point karena Anda telah menghindari menghitung kekuatan negatif yang besar. Ini berarti Anda tidak dapat lagi menggunakan fungsi evaluasi norma bawaan, tetapi jika itu bukan masalah, ini mungkin jawaban yang lebih baik. Misalnya, katakanlah kita memiliki situasi di mana
Evaluasi itu seperti yang saya sarankan, dan Anda mendapatkan -800. Namun, dalam matlab jika kita mengekspansi pengambilan log, kita mendapatkan .log( exp( - 800 ) ) = log( 0 ) = - In f
sumber