BAB III STATISTIK INFERENSI PADA RANTAI MARKOV 3.1
Pendahuluan
Pada Bab II telah dibahas mengenai rantai Markov berorde-r atau Ō(r) dan matriks peluang transisinya. Pada bagian ini, akan dibahas bagaimana menentukan orde rantai Markov dari barisan hasil pengamatan. Perhatikan barisan hasil pengamatan terhadap cuaca berikut sebagai contoh, Diketahui barisan hasil pengamatan terhadap keadaan cuaca selama 20 hari sebagai berikut. Disini keadaan cuaca dibagi kedalam tiga kategori yaitu hujan (disimbolkan dengan angka 1), berawan (2), dan cerah (3).
t
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Obs
2
2
1
1
3
1
3
2
1
2
3
3
3
2
2
2
2
1
2
1
Tabel. 1 Contoh Hasil Pengamatan
Masalah disini ialah bagaimana menentukan orde dari barisan pengamatan tersebut, dengan kata lain bagaimana menentukan apakah keadaan berawan di hari ke-20 bergantung secara langsung hanya kepada hujan di hari ke-19, atau
28
bergantung pula pada hujan di hari ke-18, atau bahkan bergantung pada keadaan di hari ke-17, dan seterusnya. Dalam melakukan pengujian terhadap rantai Markov, matriks peluang transisinya tidak harus diketahui, melainkan cukup dengan taksirannya. Penentuan penaksir titik sebuah parameter dapat ditempuh dengan menggunakan beberapa metode seperti metode momen, metode kuadrat terkecil, dan metode kemungkinan maksimum (maximum likelihood). Dari ketiga metode tersebut, dalam menaksir peluang transisi yang akan digunakan disini adalah metode kemungkinan maksimum karena metode tersebut merupakan metode yang paling banyak digunakan dibandingkan kedua metode lainnya. Gagasan dari metode ini ialah bahwa penaksir parameter yang wajar yang berdasarkan informasi sampel adalah nilai parameter yang menghasilkan peluang terbesar untuk mendapatkan sampel tersebut.
3.2
Estimasi Kemungkinan Maksimum Untuk P
Sebelumnya perhatikan cara membangun matriks yang menyatakan banyaknya transisi antar keadaan sebagai berikut. Secara umum definisikan matriks banyaknya transisi M sebagai berikut, ⎡ m11 m12 … m1n ⎤ ⎢m m22 … m2 n ⎥⎥ 21 ⎢ M= ⎢ ⎥ ⎢ ⎥ ⎣⎢ mn1 mn 2 … mnn ⎥⎦
dengan mij = M (Vt = j Vt −1 = i ), i, j = 1, 2,
(3.2.1)
,n
Disini, mij menyatakan banyaknya transisi dari keadaan i ke j. Dari Tabel 1, transisi dari keadaan 1 ke 1 hanya terjadi satu kali yaitu dari hari ke-3 ke hari ke4. Dengan demikian, m11 = 1. Dengan cara yang sama, diperoleh matriks M dari hasil observasi pada Tabel 1 sebagai berikut, ⎡ m11 M = ⎢⎢ m21 ⎢⎣ m31
m12 m22 m32
m13 ⎤ ⎡1 2 2 ⎤ m23 ⎥⎥ = ⎢⎢ 4 4 1 ⎥⎥ m33 ⎥⎦ ⎢⎣1 2 2 ⎥⎦
(c 3.1)
29
Kemudian dari matriks tersebut dicari matriks peluang transisi untuk rantai Markovnya. Peluang transisi tersebut ditaksir dengan menggunakan metode kemungkinan maksimum. Dudewicz dan Mishra [4] mendefinisikan fungsi kemungkinan (likelihood) sebagai berikut, Definisi (Fungsi Kemungkinan) Andaikan V1, V2,…,Vt adalah t buah peubah acak dengan fungsi distribusi F (v1 , v2 ,..., vt θ ) dengan θ ∈Θ merupakan parameter yang tidak diketahui, maka fungsi kemungkinan ialah,
L(θ ) = f (v1 , v2 ,..., vt θ ) Dengan mengambil V1, V2,…,Vt sampel acak yang berdistribusi identik dan saling bebas, maka diperoleh fungsi kemungkinan sebagai berikut, L(θ ) = f (v1 , v2 ,..., vt θ ) = f (v1 θ ) f (v2 θ )
f (vt θ )
Setiap θ = θ (V1 ,V2 ,...,Vt ) ∈ Θ dimana L(θ ) = sup{L(θ ) : θ ∈ Θ} , disebut penaksir kemungkinan maksimum dari θ . Penaksiran θ dengan menggunakan metode kemungkinan maksimum telah banyak dibahas dalam buku-buku statistika. Namun, secara umum dalam pembahasannya fungsi distribusi F (v θ ) maupun fungsi kepadatan peluang
f (v θ ) dari peubah acak V telah diketahui sebelumnya. Dudewicz dan Mishra [4] memberikan sebuah contoh menarik tentang bagaimana melakukan penaksiran kemungkinan maksimum dengan hanya mengetahui nilai-nilai peluang dari V untuk Ө tertentu. Contoh (Penaksiran Kemungkinan Maksimum) Diketahui tabel peluang dari
peubah acak V = {0,1,2,3,4} terhadap parameter yang berbeda-beda yaitu θ1, θ 2 , θ 3 sebagai berikut,
30
V
0
1
2
3
4
p(v I Ө1)
0,00
0,05
0,05
0,08
0,10
p(v I Ө2)
0,05
0,05
0,08
0,10
0,00
p(v I Ө3)
0,9
0,08
0,02
0,00
0,00
Di sini θ (0) = θ3 , θ (1) = θ3 , θ (2) = θ 2 , θ (3) = θ1 , θ (4) = θ1 . Andaikan kita adakan sedikit
perubahan
bahwa
kita
tahu
θ ∈ Θ = {θ1 ,θ 2 }, maka
penaksiran
kemungkinan maksimum tidak lagi tunggal karena meskipun θ (0) = θ 2 ,
θ (2) = θ 2 , θ (3) = θ1 , θ (4) = θ1 nilai θ (1) memiliki dua kemungkinan yaitu θ1 dan θ 2 . Hasil dari contoh diatas konsisten dengan definisi kita semula yaitu taksiran kemungkinan maksimum dari θ adalah θ yang memaksimumkan fungsi kemungkinan L( θ ). Selain itu, hasil taksiran tersebut tidak harus tunggal. Secara umum, untuk rantai Markov berorde 1 dengan n keadaan fungsi kemungkinannya adalah, L ( p ) = P(Vtn = in , Vtn−1 = in −1 ,
, Vt1 = i1 , Vt0 = i0 )
= P (Vtn = 1n Vtn−1 = in −1 ) × = pi0 × pi0i1 × = pi0 ×
∏
× P (Vt1 = i1 Vt0 = i0 ) × P(Vt0 = i0 )
× pin−1in m
i , j =1,2, , n
pij ij
(3.2.2)
Untuk mencari pˆ kl pertama persamaan (3.2.2) ditulis, n
n
mkw L( p ) = pi0 × ∏ pij ij × ∏ pkw , i , j =1 i≠k
m
k = 1, 2,
,n
(3.2.3)
w =1
Karena pada langkah selanjutnya akan dicari turunan pertama persamaan (3.2.3) terhadap Pkl, maka untuk mempermudah fungsi kemungkinannya dilog-kan sehingga pij dan pkw terlepas dari pangkat mij dan mkw. Hal ini dapat dilakukan karena fungsi log mempertahankan sifat kemonotonan. Dengan demikian, persamaan (3.2.3) menjadi,
31
⎛ ⎞ n n mij mkw ⎟ ⎜ log L( p ) = log pi0 × ∏ pij × ∏ pkw ⎜ ⎟ i , j =1 w =1 ⎜ ⎟ i≠k ⎝ ⎠ n
n
mkw = log pi0 + ∑ log pij ij + ∑ log pkw m
i , j =1 i≠k
l =1
n
n
i , j =1 i≠k
l =1
= log pi0 + ∑ mij log pij + ∑ mkw log pkw dimana k = 1, 2,
(3.2.4)
,n
Dari (2.2.1), maka dalam persamaan (3.2.4) bagian
n
∑m w=1
kw
log pkw
dapat
dituliskan sebagai berikut, n
n
w =1
w =1 w≠ l
∑ mkw log pkw = mkl log pkl + ∑ mkw log pkw n
n
w =1 w≠ l
s =1 s≠w
= mkl log pkl + ∑ mkw log(1 − pkl − ∑ pks ),
k , l = 1, 2,
,n
sehingga persamaan (3.2.4) menjadi,
⎛ ⎞ n n ⎜ log L( p) = log pi0 + ∑ mij log pij + mkl log pkl + ∑ mkw log(1 − pkl − ∑ pks ) ⎟ ⎜⎜ ⎟⎟ i , j =1 w=1 s =1 i ≠k w≠l s≠w ⎝ ⎠ dimana k , l = 1,2, , n n
n
n
w =1 w≠l
s =1 s≠w
(3.2.5)
Bagian selain mkl log pkl + ∑ mkw log(1 − pkl − ∑ pks ) tidak mengandung unsur pˆ kl sehingga turunan pertama (3.2.5) terhadap pˆ kl adalah,
32
⎛ ⎛ ⎞⎞ n n n ⎜ log p + ⎜ mij log pij + mkl log pkl + ∑ mkw log(1 − pkl − ∑ pks ) ⎟ ⎟⎟ ∑ i0 ⎜ ⎜⎜ ⎟⎟ i , j =1 w =1 s =1 ⎜ ⎟ i k w l s w ≠ ≠ ≠ ⎝ ⎠⎠ ⎝ ⎛ ⎞ n n ∂ ⎜ mkl log pkl + ∑ mkw log(1 − pkl − ∑ pks ) ⎟ = ⎟⎟ ∂Pkl ⎜⎜ w =1 s =1 w≠l s≠w ⎝ ⎠ n m mkw = kl − ∑ n pkl w=1 (1 − p − pks ) w≠l ∑ kl (3.2.6) s =1
∂ ∂ log L( p) = ∂Pkl ∂Pkl
s≠w
Nilai taksiran kemungkinan maksimum dari pkl adalah pˆ kl yang menjadi solusi dari
∂ log L( p) = 0 ∂Pkl ⇒
⇒
mkl pkl
n
−∑
w=1 w≠l
(3.2.7)
mkw
n
(1 − pkl − ∑ pks )
=0
s =1 s≠w
⎛ ⎞ ⎛ n ⎞ n n n ⎜ ⎟ mkl × ⎜ ∏ (1 − pkl − ∑ pks ) ⎟ − mkw × ⎜ pkl × ∏ (1 − pkl − pkw − ∑ pks ) ⎟ ⎜ w=1 ⎟ ⎜ ⎟ s =1 s =1 q =1 ⎜ w≠ l ⎟ ≠ s w s≠w q ≠l ⎜ ⎟ ⎝ ⎠ s≠q ⎝ ⎠ =0 n n pkl × ∏ (1 − pkl − ∑ pks ) s =1 s≠w
w =1 w≠l
⇒
⎛ ⎞ ⎛ n ⎞ n ⎜ ⎟ n n mkl × ⎜⎜ ∏ (1 − pkl − ∑ pks ) ⎟⎟ − mkw × ⎜ pkl × ∏ (1 − pkl − pkw − ∑ pks ) ⎟ = 0 ⎜ ⎟ s =1 s =1 q =1 ⎜ ww=≠1l ⎟ s≠w s≠w q ≠l ⎜ ⎟ ⎝ ⎠ s≠q ⎝ ⎠
(3.2.8)
Dengan menyelesaikan persamaan (3.2.8) untuk k, l = 1,2,...,n beserta (2.2.1) maka diperoleh penaksir kemungkinan maksimum untuk pkl yaitu
pˆ kl =
m kl n
∑
j =1
(3.2.9)
m kj
Dengan demikian, matriks peluang transisi taksirannya adalah
33
⎡ pˆ11 ⎢ˆ p Pˆ = ⎢ 21 ⎢ ⎢ ⎣⎢ pˆ n1
pˆ12 … pˆ1n ⎤ pˆ 22 … pˆ 2 n ⎥⎥ = pˆ n 2 …
⎥ ⎥ pˆ nn ⎦⎥
⎡ m11 ⎢ n ⎢ ∑ m1 j ⎢ j =1 ⎢ ⎢ m21 ⎢ n ⎢ ∑ m2 j ⎢ j =1 ⎢ ⎢ ⎢ mn1 ⎢ n ⎢ ∑ mnj ⎢⎣ j =1
m12
n
∑ m1 j
…
j =1
m22
n
m2 j ∑ j =1 mn 2
n
mnj ∑ j =1
m1n ⎤
n
⎥
∑ m1 j ⎥
⎥ m2 n ⎥⎥ n ⎥ m2 j ⎥ ∑ j =1 ⎥ ⎥ ⎥ mnn ⎥ n ⎥ mnj ⎥ ∑ j =1 ⎥⎦ j =1
…
…
(3.2.10)
Jadi, untuk barisan observasi pada tabel 1 matriks peluang transisinya adalah ⎡ pˆ11 Pˆ = ⎢⎢ pˆ 21 ⎢⎣ pˆ 31
pˆ12 pˆ 22 pˆ 32
⎡1 ⎢ pˆ13 ⎤ ⎢ 5 4 pˆ 23 ⎥⎥ = ⎢ ⎢9 pˆ 33 ⎥⎦ ⎢ ⎢1 ⎣⎢ 5
2 2⎤ 5 5 ⎥⎥ 4 1⎥ 9 9⎥ ⎥ 2 2⎥ 5 5 ⎦⎥
(c 3.2)
Secara umum matriks M untuk Ō(r) adalah, ⎡ m00 ⎢m ⎢ 00 ⎢ ⎢ ⎢ m0 n ⎢ m10 ⎢ ⎢ m10 M=⎢ ⎢ ⎢ m1n ⎢ ⎢ ⎢ mn 0 ⎢ ⎢ mn 0 ⎢ ⎢ ⎢⎣ mnn
00
m00
01
m00
10
m00
11
m00
n0
m0 n
n1
m0 n
00
m10
01
m10
10
m10
11
m10
n0
m1n
n1
m1n
00
mn 0
01
mn 0
10
mn 0
11
mn 0
n0
mnn
n1
mnn
⎤ ⎥ 1n ⎥ ⎥ ⎥ nn ⎥ ⎥ 0n ⎥ 1n ⎥ ⎥ ⎥ nn ⎥ ⎥ ⎥ 0n ⎥ ⎥ nn ⎥ ⎥ ⎥ nn ⎥ ⎦ 0n
(3.2.11)
34
Disini mij…kl menyatakan banyaknya transisi dari keadaan i ke l dalam r langkah melalui keadaan j pada satu langkah pertama dan seterusnya hingga melalui k di langkah ke r-1 sebelum ke l . Secara matematis dapat ditulis, mij
kl
= M (Vt = l Vt −1 = k ,
, Vt − ( r −1) = j , Vt − r = i ) , i, j , k , l = 1, 2,
,n
(3.2.12)
Sementara itu matriks peluang transisi taksirannya adalah, ⎡ pˆ 00 ⎢ pˆ ⎢ 00 ⎢ ⎢ ⎢ pˆ 0 n ⎢ pˆ 10 ⎢ ⎢ pˆ 10 Pˆ = ⎢ ⎢ ⎢ pˆ 1n ⎢ ⎢ ⎢ pˆ n 0 ⎢ˆ ⎢ pn 0 ⎢ ⎢ ⎣⎢ pˆ nn
pˆ 00 pˆ 00
00 10
pˆ 00 pˆ 00
01 11
00
pˆ 0 n pˆ 10
01
pˆ 0 n pˆ 10
10
pˆ 10
11
pˆ 10
n0
pˆ 1n
n1
pˆ 1n
n0
n1
10
pˆ n 0 pˆ n 0
11
pˆ n 0 pˆ n 0
n0
pˆ nn
n1
pˆ nn
00
01
⎤ ⎥ 1n ⎥ ⎥ ⎥ nn ⎥ ⎥ 0n ⎥ 1n ⎥ ⎥ ⎥ nn ⎥ ⎥ ⎥ 0n ⎥ ⎥ nn ⎥ ⎥ ⎥ ⎥ nn ⎦ 0n
(3.2.13)
dimana
pij
kl
= P(Vtn = l Vt −1 = k , ,Vt −( r −1) = j,Vt −r = i)
(3.2.14)
dan
pˆ ij
kl
=
mij
kl
(3.2.15)
n
∑ mij h =1
kh
Berdasarkan (3.2.11) dan (3.2.13), untuk tabel 1 maka matriks banyak transisi M dan matriks peluang transisi taksiran Pˆ Ō(2) adalah,
35
⎡ m111 ⎢m ⎢ 121 ⎢ m131 ⎢ ⎢ m211 M = ⎢ m221 ⎢ ⎢ m231 ⎢m ⎢ 311 ⎢ m321 ⎢m ⎢⎣ 331
m112 m122 m132 m212 m222 m232 m312 m322 m332
m113 ⎤ ⎡ 0 m123 ⎥⎥ ⎢ 1 ⎢ m133 ⎥ ⎢ 1 ⎥ ⎢ m213 ⎥ ⎢ 1 m223 ⎥ = ⎢ 2 ⎥ ⎢ m233 ⎥ ⎢ 0 m313 ⎥ ⎢ 0 ⎥ ⎢ m323 ⎥ ⎢ 1 m333 ⎥⎥⎦ ⎢⎢⎣ 0
0 1⎤ 0 1 ⎥⎥ 1 0⎥ ⎥ 2 0⎥ 2 0⎥ ⎥ 0 1⎥ 0 1⎥ ⎥ 1 0⎥ 1 1 ⎥⎥⎦
(c 3.3)
dan
⎡ pˆ111 ⎢ pˆ ⎢ 121 ⎢ pˆ131 ⎢ ⎢ pˆ 211 Pˆ = ⎢ pˆ 221 ⎢ ⎢ pˆ 231 ⎢ pˆ ⎢ 311 ⎢ pˆ 321 ⎢ pˆ ⎣⎢ 331
pˆ112 pˆ122 pˆ132 pˆ 212 pˆ 222 pˆ 232 pˆ 312 pˆ 322 pˆ 332
⎡0 ⎢1 ⎢ ⎢2 pˆ113 ⎤ ⎢ 1 pˆ123 ⎥⎥ ⎢ ⎢2 pˆ133 ⎥ ⎢ 1 ⎥ pˆ 213 ⎥ ⎢ 3 ⎢ pˆ 223 ⎥ = ⎢ 1 ⎥ pˆ 233 ⎥ ⎢ 2 ⎢ pˆ 313 ⎥ ⎢ 0 ⎥ ⎢ pˆ 323 ⎥ ⎢ 0 pˆ 333 ⎥⎦⎥ ⎢ 1 ⎢2 ⎢ ⎢0 ⎣⎢
0 0 1 2 2 3 1 2 0 0 1 2 1 2
1⎤ 1 ⎥⎥ 2⎥ ⎥ 0⎥ ⎥ ⎥ 0⎥ ⎥ ⎥ 0⎥ ⎥ 1⎥ 1⎥ ⎥ 0⎥ ⎥ ⎥ 1⎥ 2 ⎦⎥
(c 3.4)
Menentukan berapa orde dari barisan pengamatan pada tabel 1 sama dengan menentukan matriks peluang transisi mana yang akan digunakan, yaitu apakah (c 3.2) atau (c 3.4). Disisi lain, bisa saja ke-20 hasil pengamatan tersebut ternyata saling bebas. Dalam hal ini,
pij
kl
= P(Vt = l Vt −1 = k , ,Vt −( r −1) = j,Vt −r = i) = P(Vt = l ) = pl
Untuk kasus dimana ke-20 hasil pengamatan saling bebas, matriks banyak transisi dan matriks peluang transisi taksirannya adalah,
M = [ m1 m2
m3 ] = [ 6 9 5]
(c 3.5)
dan Pˆ = [ pˆ1
pˆ 2
⎡6 pˆ 3 ] = ⎢ ⎣ 20
9 20
5⎤ 20 ⎥⎦
(c 3.6)
36
Dengan adanya (c3.6), alternatif matriks peluang transisi untuk barisan hasil perngamatan pada tabel 1 bertambah. Sebelum memeriksa matriks peluang transisi –beserta orde yang berpadanan– mana yang sesuai, ada baiknya periksa terlebih dahulu mana matriks peluang transisi yang tidak dapat digunakan untuk data yang dimiliki. Perhatikan jika barisan pengamatan pada tabel 1 diperlakukan sebagai Ō(3). Matriks M-nya adalah, ⎡ m1111 ⎢m M = ⎢ 1121 ⎢ ⎢ ⎢⎣ m3331
m1112 m1122 m3332
m1113 ⎤ ⎡ 0 0 0 ⎤ m1123 ⎥⎥ ⎢ 0 0 0 ⎥ ⎥ =⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ m3333 ⎥⎦ ⎣ 0 1 0 ⎦
(c 3.7)
Baris pertama dan kedua dari M bernilai 0. Oleh karena itu, matriks peluang transisi taksirannya adalah, ⎡ pˆ1111 ⎢ pˆ Pˆ = ⎢ 1121 ⎢ ⎢ ⎢⎣ pˆ 3331
pˆ1112 pˆ1122 pˆ 3332
⎡0 ˆp1113 ⎤ ⎢ 0 ⎢ pˆ1123 ⎥⎥ ⎢ 0 = ⎥ ⎢0 ⎥ ⎢ pˆ 3333 ⎥⎦ ⎢ ⎢⎣ 0
0 0 0 0 1
0⎤ 0⎥ ⎥ 0⎥ 0⎥ ⎥ ⎥ 0 ⎥⎦
(c 3.8)
Perhatikan bahwa ada baris dimana jumlah seluruh elemennya tidak sama dengan 1, yaitu baris 1 dan 2. Oleh karena itu, berdasarkan (2.2.1) dapat disimpulkan bahwa matriks pada (c.3.8) bukan merupakan matriks peluang transisi suatu rantai Markov. Jadi, data pengamatan pada tabel 1 tidak dapat dimodelkan sebagai Ō(r) dengan r > 2. Satu hal yang dapat disimpulkan dari ilustrasi diatas ialah, penentuan orde pada Ō(r) juga dipengaruhi oleh data barisan pengamatan yang dimiliki. Semakin tinggi orde, perulangan yang terjadi semakin sedikit. Hal ini terlihat dari matriks banyaknya transisi dimana pada contoh diatas, semakin tinggi orde yang digunakan, semakin kecil elemen-elemen matriks banyak transisinya. Oleh karena itu, semakin tinggi orde yang akan digunakan, semakin banyak pula data yang
37
diperlukan. Dalam praktek, pada kenyataannya orde tersebut secara umum tidak lebih dari lima. Dalam aplikasi, penentuan orde suatu rantai Markov seringkali bersifat subjektif tergantung pemakai. Meskipun besarnya orde dapat meningkatkan kecocokan model terhadap kenyataan sebenarnya, demi kemudahan perhitungan –juga saat membangun model– orde diharapkan serendah mungkin (parsimoni). Meski demikian, beberapa batasan dapat dipertimbangkan misalnya, (i)
Kemudahan dalam memperoleh data.
(ii)
Jika indeks parameter berkaitan dengan waktu, perhatikan indeks datanya. Jika data yang diambil merupakan data bulanan, maka indeks parameternya juga dalam bulan.
Masalah berikutnya adalah bagaimana menentukan orde rantai Markov yang tepat. Dalam tugas akhir ini, penentuan orde tersebut dilakukan bukan dengan menaksir, melainkan dengan cara memeriksa atau ’membandingkan’ orde mana yang paling sesuai dari beberapa alternatif pilihan orde yang memungkinkan. Sebagai contoh, untuk tabel 1 akan dibandingkan antara Ō(0), Ō(1), dan Ō(2).
3.3
Studi Deskriptif Ō(r)
Dalam melakukan pengujian terhadap orde rantai Markov, perlu dipahami terlebih dahulu bagaimana pengaruh orde pada suatu rantai Markov. Dari hal tersebut, dapat diperoleh karakteristik yang membedakan orde yang satu dengan lainnya. Selanjutnya, karakteristik tersebut dapat dijadikan dasar dalam pengujian orde rantai Markov. Sebagai contoh, berikut ini adalah plot dari Ō(0) dan Ō(1) dengan tiga keadaan dan t=100,
38
3
Vt 2
1 0
5
10
15
20
25
30
35
40
45
50
t
55
60
65
70
75
80
85
90
95 100
60
65
70
75
80
85
90
95 100
Gambar. 1 Plot Ō(0)
3
Vt 2
1 0
5
10
15
20
25
30
35
40
45
50
t
55
Gambar. 2 Plot Ō(1)
Secara deskriptif, terlihat perbedaan yang menarik antar rantai Markov dengan orde yang berbeda. Dari Gambar 1, perhatikan pola transisi berbentuk 3
2
(c 3.9)
1 1
2
3
dan 3
(c 3.10)
2
1 1
1.5
2
2.5
3
yang muncul secara berulang. Sementara itu, perhatikan pula pada Gambar 2 pola transisi yang berbentuk, 3
(c 3.11)
2
1 1
2
3
4
5
dan
39
3
(c 3.12)
2
1 1 1.5 2 2.5 3 3.5 4 4.5 5
yang juga muncul beberapa kali. Pola-pola tersebut merupakan pola-pola transisi yang dominan dari kedua plot rantai Markov di atas. Dari ilustrasi tersebut dapat dibuat suatu dugaan bahwa semakin tinggi orde rantai Markov, maka semakin ’lebar’ pola transisinya. Hal ini sangat masuk akal karena kemunculan suatu keadaan bergantung secara langsung kepada r-keadaan sebelumnya dimana besar r menyatakan orde dari rantai Markov yang terkait. Secara matematis, pola-pola transisi tersebut dapat dinyatakan sebagai elemen dari matriks banyak transisi M. Pola transisi (c 3.9) dapat ditulis sebagai m131, (c 3.10) m313, (c 4.11) m13331, dan (c.3.12) m13232. Dengan demikian, bentuk plot suatu rantai Markov dapat digambarkan oleh matriks banyak transisinya.
3.4
Uji Kesesuaian Suatu Rantai Markov
Seperti yang telah diutarakan sebelumnya, orde rantai Markov ditentukan dengan cara memilih orde yang paling sesuai dari beberapa alternatif pilihan orde yang memungkinkan. Kemudian, yang perlu diperhatikan adalah sangat diharapkan bahwa orde bernilai serendah mungkin. Dari sini diperoleh 2 gambaran mengenai kriteria orde yang tepat yaitu, (i)
Semakin besar orde, semakin cocok model yang dimiliki dengan kenyataan sebenarnya.
(ii)
Orde rantai Markov yang rendah lebih diharapkan daripada orde yang lebih tinggi karena akan membuat model yang dihasilkan lebih sederhana.
Dari kedua poin diatas dapat dikatakan bahwa lebih baik memilih orde yang lebih rendah bila gambaran mengenai rantai makov yang direpresentasikannya tidak berbeda terlalu signifikan jika dibandingkan dengan orde yang lebih tinggi. Dari paragraf diatas, diperoleh langkah awal dalam melakukan uji kesesuaian. Pertama, asumsikan terlebih dahulu bahwa rantai Markov tersebut Ō(r) dimana r menyatakan orde maksimum yang dikehendaki oleh peneliti atau orde yang masih dapat ditoleransi oleh data. Untuk data pada tabel 1, r = 2 karena untuk r = 3,
40
tidak terdapat matriks peluang transisinya. Selanjutnya, bandingkan dengan apa yang terjadi jika rantai Markov tersebut Ō(r-1). Dengan demikian, hipotesisnya adalah, H0 : penggunaan Ō(r) tidak berbeda dengan Ō(r-1). H1 : penggunaan Ō(r) berbeda dengan Ō(r-1). Berdasarkan definisi orde rantai Markov, hipotesis tersebut dapat ditulis sebagai, H 0 : P (Vt = l Vt −1 = k ,
, Vt −( r −1) = j , Vt − r = i) = P(Vt = l Vt −1 = k ,
, Vt −( r −1) = j )
H 0 : pij ...kl = p j ...kl
(3.3.1)
dan H1 : P(Vt = l Vt −1 = k ,
,Vt −( r −1) = j ,Vt − r = i ) ≠ P(Vt = l Vt −1 = k ,
, Vt −( r −1) = j )
H1 : pij ...kl ≠ p j ...kl
(3.3.2)
Berdasarkan Ō(r), transisi dari keadaan i ke l dalam r langkah melalui keadaan j pada satu langkah pertama dan seterusnya hingga melalui k di langkah ke r-1 sebelum ke l terjadi sebanyak mij…kl kali dimana mij…kl merupakan elemen dari matriks banyaknya transisi M untuk Ō(r). Sementara itu, mengacu kepada (3.2.15), berdasarkan Ō(r-1) transisi tersebut terjadi sebanyak, mˆ ij
dimana pj
kl
kl
= mij
k
× pj
mij
= p(Vt = l Vt −2 = k , ,Vt −(r −1) = j ) , dan mij
k
kl
= m(Vt −1 = k n
k
= ∑ mij l =1
kl
,Vt −( r −1) = j,Vt −r = i) ,
.
Pada paragraf sebelumnya, telah dibahas bahwa dalam melakukan pemilihan antara Ō(r) dengan Ō(r-1), bila gambaran mengenai rantai makov yang direpresentasikannya tidak berbeda terlalu signifikan, orde yang lebih rendah dipilih. Sementara itu, gambaran mengenai rantai Markov dapat direpresentasikan oleh banyaknya peluang transisi m. Dengan demikian, hipotesis H0 : penggunaan
41
Ō(r) tidak berbeda dengan Ō(r-1), akan didukung bila selisih banyaknya transisi antara Ō(r) dan Ō(r-1) ’kecil’ yaitu, n
n
∑∑ i =1 j =1
n
n
∑∑ mij k =1 l =1
kl
n
− mˆ ij
kl
n
= ∑∑ i =1 j =1
n
n
∑∑ m
ij kl
k =1 l =1
− (mij
k
pj
kl
) ≤c
(3.4.1)
untuk suatu konstanta c. Masalah berikutnya adalah menentukan seberapa besar c sehingga selisih pada (3.4.1) dapat dianggap ’kecil’. Sebelumnya, perhatikan sebuah Percobaan Bernoulli dengan dua kemungkinan keluaran, yaitu A1 dan A2 dengan peluangnya masing-masing adalah P(A1) = p1 dan P(A2) = p2 = 1-p1. Suatu sampel acak dari m kali percobaan diamati dimana m1 dan m2 = m - m1 masing-masing menyatakan banyaknya keluaran dari jenis A1 dan A2. Berdasarkan teorema limit pusat untuk m → ∞, (m1 − mp1 ) ∼ N (0,1) mp1 (1 − p1 )
Z= Tulis Y = Z2, maka
G ( y ) = P(Y ≤ y ) = P( Z 2 ≤ y ) = P(− y ≤ Z ≤ 0 ⎧⎪ =⎨ ⎪⎩ F ( y ) − F (− y )
y)
y<0 y≥0
(3.4.2)
Dari, (3.4.2) untuk y ≥ 0 fungsi kepadatan peluang Y, g ( y) =
(
∂G ( y ) ∂ F ( y ) − F (− y ) = ∂y ∂y
)
=
1 2 y
(f(
y ) + f (− y )
)
y − 1 ⎛ 1 − 2y 1 − 2y ⎞ 1 ⎛ 1 − 2y 1 − 2y ⎞ 1 2 = e + e ⎟= e ⎜ ⎜ e + e ⎟= 2 2 y ⎝ 2π 2π 2π y ⎝ 2 2π y ⎠ ⎠
(3.4.3)
Dari definisi fungsi beta dimana 1
B ( p, q ) = ∫ x p −1 (1 − x) q −1 dx = 0
Γ ( p )Γ ( q ) Γ( p + q)
untuk p = q dan w = 2x-1 maka, 1
Γ 2 ( p) = 22− 2 p ∫ (1 − y 2 ) p −1 dy Γ(2 p ) 0
sehingga untuk p = 0.5 berlaku
42
1
dy ⎛1⎞ ⎛1⎞ Γ ⎜ ⎟ = Γ2 ⎜ ⎟ = 2 × ∫ 2 ⎝2⎠ ⎝2⎠ 0 1− y ⎛ = 2 × ⎜ arc(sin y ) ⎜ ⎝
⎞ π ⎟⎟ = 2 × 2 0⎠
1
(3.4.4)
= π Berdasarkan (3.4.4) maka fungsi kepadatan peluang Y pada (3.4.3) dapat ditulis, y 1 y − − − 1 1 2 2 2 g ( y) = e = y e 2π y 2 π
y<0
0 ⎧ ⎪ 1 y − − ⎪ = ⎨ 1 y 2e 2 ⎪ 2Γ ⎛ 1 ⎞ ⎪⎩ ⎜⎝ 2 ⎟⎠
y≥0
yang merupakan fungsi kepadatan peluang distribusi khi kuadrat dengan derajat kebebasan 1. Jadi, Y=
(m1 − mp1 ) 2 ∼ χ 2 (1) mp1 (1 − p1 )
(3.4.5)
Dengan sedikit manipulasi aljabar, (3.4.5) dapat diubah menjadi, (m1 − mp1 ) 2 mp1 (1 − p1 )
=
(m1 − mp1 ) 2 (m1 − mp1 ) 2 + mp1 m(1 − p1 )
(m1 − mp1 ) 2 ( ( m − m1 ) − m (1 − p1 ) ) = + mp1 m(1 − p1 ) (m − mp1 ) 2 ( m2 − mp2 ) = 1 + mp1 mp2
2
2
(mi − mpi ) 2 ∼ χ 2 (1) =∑ mpi i =1 2
(3.4.6)
Dengan melihat (3.4.6), perhatikan bahwa selisih banyaknya transisi pada (3.4.1) dapat dituliskan sebagai, n
n
n
n
∑∑ ∑∑ i =1 j =1
k =1 l =1
(m
ij kl
− (mij mij
k pj
k
pj
kl
))
2
≤ c∗
(3.4.7)
kl
dengan c* untuk suatu konstanta real sembarang. Untuk memutuskan seberapa besar c* sehingga selisih tersebut dapat dikatakan ’kecil’, perlu dketahui terlebih
43
dahulu distribusi dari (3.4.7). Dengan mengacu kepada (3.4.6) tentu wajar jika ada dugaan bahwa (3.4.7) juga berdistribusi khi kuadrat dengan derajat kebebasan tertentu, sebut saja u. Tentunya hal ini perlu dibuktikan secara matematis. Akan tetapi, pembuktian disini tidak akan disajikan secara mendetail karena diperlukan dasar-dasar aljabar yang kuat dalam penurunannya. Pertama, misalkan di setiap observasi dari n observasi yang menghasilkan n titik hasil pengamatan dalam sampel, terdapat peluang pi bahwa hasil yang diperoleh berasal dari himpunan Si. Untuk setiap himpunan dari bilangan bulat non-negatif v1, ..., vn dimana
r
∑v i =1
i
= n , peluang dimana dari n observasi terdapat tepat vi hasil
observasi yang berasal dari Si untuk i = 1, ..., r, diberikan oleh, n! v1 ! vr !
p1v1
prvr
yang merupakan bentuk umum dari (p1 + ... + pr)n. Jadi, distribusi dari r grup frekuensi v1, ..., vn adalah perumuman dari distribusi binomial yang dikenal sebagai distribusi multinomial. Fungsi karakteristik gabungannya adalah, r ⎛ ⎞ E ⎜ exp ∑ it j x j ⎟ = ( p1eit1 + j =1 ⎝ ⎠
Tulis xi =
ϕ (t1 ,
+ pr eitr )
n
vi − npi , maka fungsi karakteristik gabungan dari xi adalah, npi ⎛ r r ⎛ v j − np j ⎛ ⎞ , tr ) = E ⎜ exp ∑ it j x j ⎟ = E ⎜ exp ∑ it j ⎜ ⎜ np j ⎜ j =1 j =1 ⎝ ⎠ ⎝ ⎝ ⎛ r ⎛ r t = E ⎜ exp ⎜ ∑ i j v j − ∑ it j ⎜ j =1 np j ⎜ j =1 ⎝ ⎝ ⎛ i = ⎜ p1e ⎜ ⎝
t1 np1
+
+ pr e
i
tr npr
n
⎞⎞ ⎟⎟ ⎟⎟ ⎠⎠ r ⎛ ∑i ⎞⎞ np j ⎟ ⎟ = E ⎜ e j=1 ⎜⎜ ⎟⎟ ⎠⎠ ⎝
tj np j
⎞ − ∑ it j ⎟ e j=1 ⎟⎟ ⎠ r
vj
np j
r
⎞ −i n ∑ t j j =1 ⎟ e ⎟ ⎠
pj
(3.4.8)
Dengan uraian MacLaurin, untuk t1,..., tn konstan,
44
log ϕ (t1 ,
r ⎡ i r 1 r 2 −3 ⎤ , tr ) = −i n ∑ t j p j + n log ⎢1 + t p t j + O(n 2 ) ⎥ − ∑ ∑ j j 2n j =1 n j =1 j =1 ⎣ ⎦ 2
⎞ 1 r 1⎛ r −1 = − ∑ t j 2 + ⎜ ∑ t j p j ⎟ + O(n 2 ) 2 j =1 2 ⎝ j =1 ⎠
(3.4.9)
Berdasarkan (3.4.9) maka (3.4.8) menjadi,
ϕ (t1 ,
−
, tr ) = e
1 2
1⎛ t j2 + ⎜ 2⎜ j =1 ⎝ r
∑
2
r
∑ j =1
tj
⎞ −1 p j ⎟ +O ( n 2 ) ⎟ ⎠
(3.4.10)
Selanjuntya untuk n → ∞,
lim ϕ (t1 , , tr ) = e
⎡ r ⎛ 1 − ⎢⎢ t j 2 −⎜ ⎜ 2 j =1 ⎝ ⎢⎣
∑
r
∑t j j =1
2 ⎞ ⎤ p j ⎟ ⎥⎥ ⎟ ⎠ ⎥⎦
n →∞
=e
1 − Q ( t1 , ,tr ) 2
(3.4.11)
Bentuk kuadratik
⎛ r ⎞ , tr ) = ∑ t j − ⎜ ∑ t j p j ⎟ j =1 ⎝ j =1 ⎠ r
Q(t1 , dapat
ditulis
t t = (t1 t2
sebagai
perkalian
2
(3.4.12)
2
matriks
dan
vektor
tt Λ t
dimana
tr ) dan,
Λ r×r = I − pp t ⎛1 0 ⎜ 0 1 =⎜ ⎜ ⎜⎜ ⎝0 0
0 ⎞ ⎛ p1 ⎟ ⎜ 0 ⎟ ⎜ p1 p2 − ⎟ ⎜ ⎟ ⎜ 1 ⎟⎠ ⎜ p p ⎝ 1 r
p1 p2 p2 p2 pr
p1 pr ⎞ ⎟ p2 pr ⎟ ⎟ ⎟ pr ⎟⎠
(3.4.13)
Disini Q(t1, t2, ..., tr) adalah non-negatif dengan rank r-1 (Cramer [6] hal 109-110) dan matriks Λ memiliki r-1 nilai karakteristik bernilai 1 dan nilai karakteristik ker bernilai nol. Seiring dengan n → ∞, fungsi karakteristik gabungan x1, x2, ..., xr menuju bentuk (3.4.11) yang merupakan fungsi karakteristik gabungan distribusi normal singular dengan rank r-1 dan total massa yang terletak pada hyperplane
∑x
j
pj = 0 .
Dengan demikian, secara limit, x1, x2, ..., xr terdistribusi normal singular dengan mean nol dan matriks momen Λ seperti pada (3.3.13).
45
Selanjutnya, untuk x1, x2, ..., xr berdistribusi normal dengan mean nol dan matriks y = Cx yang menggantikan
momen Λ, maka terdapat transformasi ortogonal
variabel lama x = {x1, x2, ..., xr} dengan variabel baru y = {y1, y2, ..., yr} dimana matriks momen transformasinya, B = CΛCt adalah matriks diagonal dengan r-1 eleman diagonal bernilai 1 dan satu elemen lainnya bernilai nol. Artinya, variabel y1, y2, ..., yr-1 berdistribusi normal baku dan yr berdistribusi normal dengan mean dan variansi 0 (Cramer [6] hal 313). Dengan demikian berdasarkan (3.4.5) maka, r
∑x i =1
2
i
r
r −1
i =1
i =1
= ∑ yi 2 = ∑ yi 2 ∼ χ r −1
(3.4.14)
Dengan demikian, berdasarkan (3.4.14) statistik uji (3.4.7) berdistribusi khi kuadrat dengan derajat kebebasan nr+1 – 1 yaitu, n
n
n
n
∑∑ ∑∑ i =1 j =1
k =1 l =1
(m
ij kl
− (mij mij
k
pj
k
pj kl
kl
))
2
∼ χ nr +1 −1
(3.4.15)
Selanjutnya, pada bab berikutnya akan dipaparkan langkah-langkah dalam menggunakan statistik uji (3.4.15) dalam menguji orde barisan basa nukleotida dari spesies Homo Sapiens.
46
