Contoh Kasus untuk EM Algorithm dengan MATLAB


Dari suatu eksperimen diketahui bahwa 1/3 populasi memiliki gen A, 1/3 populasi memiliki gen B dan sisanya gen C. Misalkan peluang pasutri memiliki gen yang sama adalah phi, dan gen A mendominasi gen B dan C, serta gen B mendominasi C dan seorang anak akan mendapatkan gen dominan dari kedua orang tuanya. Sehingga diketahui informasi sebagai berikut:

Ibu

Ayah

Anak

Peluang

A

A

A

(1/3) phi

A

B

A

(1/6)(1-phi)

A

C

A

(1/6)(1-phi)

B

A

A

(1/6)(1-phi)

B

B

B

(1/3) phi

B

C

B

(1/6)(1-phi)

C

A

A

(1/6)(1-phi)

C

B

B

(1/6)(1-phi)

C

C

C

(1/3) phi

Dari 100 sampel anak, 45 anak memiliki gen A, 32 anak memiliki gen B dan 23 memiliki gen C. Estimasi nilai phi! (petunjuk gunakan algoritma EM)

Pembahasan

P(anak gen A) = P(y1) = (1/3)π + (2/3)(1-π)

P(anak gen B) = P(y2) = (1/3)π + (1/3)(1-π)=1/3

P(anak gen C) = P(y3)= (1/3)π

Sehingga:

Misal nA, nb, dan nc adalah jumlah anak yang masing-masing memiliki gen A, B, dan C.

Misal kita pecah kejadian A ke dalam dua kelompok baru, D dan E, dengan peluang masing-masing adalah (1/3)π dan (2/3)π, dan nA = nD + nE . Sehingga fungsi diatas menjadi:

dimana y1 = x1 + x2; y2 = x3; dan y3 = x4.

Dengan EM algorithm:

  1. E step. E step mengestimasi sufficient statistics dari complete data

    x, dengan syarat data teramati y. Dalam contoh ini (x3, x4, x5) diketahui, sehingga sufficient statistics yang harus diestimasi hanya untuk x1 dan x2, dimana x1+x2= y1=45. Dengan demikian estimasi dari x1 dan x2 adalah:

  2. M step. M step kemudian menggunakan hasil estimasi complete data (x1(p), x2(p), x3, x4) untuk mengestimasi π menggunakan maksimum likelihood dengan memperlakukan hasil estimasi complete data sebagai data pengamatan.

, dengan K adalah suatu konstanta

Dengan memaksimumkan fungsi diatas, diperoleh:

Algoritma EM dalam contoh ini yaitu dengan melakukan berulang-ulang langkah diatas, sehingga nilai π konvergen pada presisi yang ditentukan.

MATLAB code dari kasus di atas adalah sebagai berikut:

n=[ 45 32 23];

pi0=0.5;

prec=0.001;

pi=pi0;

diff=1;

while(diff>prec)

piold=pi;

d=n(1)*pi/(pi+2*(1-pi));

e=n(1)-d;

pi=(n(3)+d)/(n(3)+d+e);

diff=abs(pi-piold);

end

pi

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: