ESTIMASI MODEL TERBAIK BANYAKNYA GEMPA BUMI MENGGUNAKAN POISSON HIDDEN MARKOV MODELS DAN ALGORITMA EM (Studi Kasus Banyaknya Gempa Bumi di Wilayah Sumatra Barat sampai Nusa Tenggara Timur)
ROSYID SURYANDARU
PROGRAM STUDI MATEMATIKA FAKULTAS SAINS DAN TEKNOLOGI UNIVERSITAS ISLAM NEGERI SYARIF HIDAYATULLAH JAKARTA 2015 M/1436 H
ESTIMASI MODEL TERBAIK BANYAKNYA GEMPA BUMI MENGGUNAKAN POISSON HIDDEN MARKOV MODELS DAN ALGORITMA EM (Studi Kasus Banyaknya Gempa Bumi di Wilayah Sumatra Barat sampai Nusa Tenggara Timur)
Skripsi Sebagai Salah Satu Syarat Untuk Memperoleh Gelar Sarjana Sains Fakultas Sains dan Teknologi Universitas Islam Negeri Syarif Hidayatullah Jakarta
Oleh: Rosyid Suryandaru 108094000024
PROGRAM STUDI MATEMATIKA FAKULTAS SAINS DAN TEKNOLOGI UNIVERSITAS ISLAM NEGERI SYARIF HIDAYATULLAH JAKARTA 2015/1436 H
i
PENGESAHAN UJIAN Skripsi berjudul “Estimasi Model Terbaik Banyaknya Gempa Menggunakan Poisson Hidden Markov Models dan Algoritma EM” yang ditulis oleh Rosyid Suryandaru, NIM 108094000024 telah diuji dan dinyatakan lulus dalam sidang Munaqosah Fakultas Sains dan Teknologi Universitas Islam Negeri Syarif Hidayatullah Jakarta pada hari Kamis, 8 Januari 2015. Skripsi ini telah diterima untuk memenuhi salah satu persyaratan dalam memperoleh gelar Sarjana strata satu (S1) Program Matematika. Menyetujui: Penguji 1,
Penguji 2,
Dr. Nina Fitriyati, M.Kom NIP. 19760414 200604 2 001
Irma Fauziah, M.Sc NIP. 19800703 201101 2 005
Pembimbing 1,
Pembimbing 2,
Suma’inna, M.Si NIP. 19791208 200701 2 015
Mahmudi, M.Si
Mengetahui, Dekan Fakultas Sains dan Teknologi,
Ketua Program Studi Matematika,
Dr. Agus Salim, M.Si NIP. 19720816 199903 1 003
Yanne Irene, M.Si NIP. 19741231 200501 2 018
ii
PERNYATAAN
DENGAN INI SAYA MENYATAKAN BAHWA SKRIPSI INI BENARBENAR HASIL KARYA SENDIRI YANG BELUM PERNAH DIAJUKAN SEBAGAI SKRIPSI ATAU KARYA ILMIAH PADA PERGURUAN TINGGI ATAU LEMBAGA MANAPUN.
Jakarta, 8 Januari 2015
Rosyid Suryandaru 108094000024
iii
PERSEMBAHAN Alhamdulillahirobbil’alamin, segala puji bagi Allah, Tuhan Semesta Alam Skripsi ini saya persembahkan untuk Ibuku tercinta Sri Kus Indrati, Kakakkakakku, Pilon dan seluru keluarga besarku, dan juga Keluarga besar Prodi Matematika, serta semua pihak yang terlibat di dalamnya. Semoga selalu diridhoi Allah SWT, selalu dalam lindungan-Nya, serta selalu dibukakan pintu rahmat, kasih sayang, dan hidayah-Nya Amin
MOTTO
A good student never will forget a good teacher :’)
Lebih baik salah karena melakukan sesuatu dari pada tidak pernah salah karena tidak pernah melakukan sesuatu
iv
ABSTRAK Rosyid Suryandaru, Estimasi Model Terbaik Banyaknya Gempa Bumi Menggunakan Algoritma EM dan Poisson Hidden Markov Models. Di bawah bimbingan Suma’inna, M.Si dan Mahmudi, M.Si.
Penelitian dilakukan dengan dilatarbelakangi oleh wilayah Indonesia yang merupakan salah satu wilayah yang rawan akan gempa bumi. Secara umum, banyaknya peristiwa gempa bumi memiliki karakteristik distribusi Poisson. Akan tetapi terjadi overdispersi relatif terhadap distribusi poisson. Salah satu cara dalam mengatasi overdispersi khususnya pada distribusi poisson adalah dengan menggunakan Poisson Hidden Markov Models (PHMMs), yang kemudian mengaplikasikan algoritma EM (Expectation Maximisation Algorithm) pada masing-masing model untuk memdapatkan parameter estimasinya. Dari modelmodel yang diperoleh, selanjutnya dipilih model terbaik berdasarkan nilai AIC (Akaike Information Criterion) terkecil. Data yang digunakan adalah data sekunder, yaitu data banyaknya gempa bumi di wilayah Sumatra Barat sampai Nusa Tenggara Timur tahun 2008-2013 dengan magnitudo ≥5 Skala Richter yang terjadi pada kedalaman gempa bumi dangkal yaitu ≤60 Km. Dari penelitia yang dilakukan didapat 3 model PHMM, yaitu model dengan keadaan tersembunyi, yaitu 𝑚 = (2,3,4). Berdasarkan nilai AIC, model dengan 3 keadaan tersembunyi merupakan model estimasi banyaknya gempa bumi terbaik dengan nilai AIC sebesar 537,429. Hasil estimasi parameter dari model terbaik PHMM, yaitu model dengan 3 keadaan tersembunyi, nilai estimasi parameter ratarata banyaknya gempa bumi yang terjadi sebanyak 2,8446121 ≈ 3 peristiwa dalam kurun waktu 15 hari.
Kata kunci: Overdispersi, Mixture Distribution, Rantai Markov, Peluang forward backward, PHMMs, algoritma EM, AIC.
v
ABSTRACT
Rosyid Suryandaru, The Best Model Estimation for Many Earthquake Happen with Using EM Algorithm in Poisson Hidden Markov Models. Under guidance Suma’inna, M.Si and Mahmudi, M.Si.
This research was conducted with Indonesian zone that one of zone have often earthquake happen. In other hand, many earthquake happen have a characteristic poisson of distribution. But often be relative overdisversion. One of method to solve overdisversion on poisson distribution specifically is with using Poisson Hidden Markov Models (PHMMs) that so applied EM algorithm (Expectation Maxsimisation Algorithm) on every models to generated estimation value. From every models resulted, next selecting the best model obtained from the least AIC value (Akaike Information Criterion). In this thesis, using secondary data, that is data of many earthquake happen in Sumatra Barat till Nusa Tenggara Timur on 2008-2013 with magnitude ≥ 5 SR that happen in earthquake low deepen is ≤ 60 km. From this research exist 3 PHMM models, that is hidden model happen, with m= (2,3,4). Obtained AIC value, 3 of hidden model happen is the best model estimation of many earthquake happen with weight of AIC value is 537,429. Estimation parameter result from the best model of PHMM is 3 of hidden happen model, average of estimation parameter value in many earthquake happen is 2,8446121≈3 moment of 15 day.
Keywords: Overdisversion, Mixture Distribution, Markov Chain, Forward Backward Probability, PHMMs, EM Algirithm, AIC.
vi
KATA PENGANTAR بسم هللا اار حمن اار حيم
Assalamu’alaikum Warahmatullahi Wabarakatuh Alhamdulillah, Segala puji bagi Allah, Tuhan Semesta Alam, yang senantiasa melimpahkan rahmat dan nikmat-Nya kepada kita semua, tak terkecuali pada penulis, hingga penulis dapat menyelesaikan skripsi “Estimasi Model Terbaik Banyaknya Gempa Bumi Menggunakan Poisson Hidden Markov Models dan Algoritma EM” dengan studi kasus “Banyaknya Gempa Bumi di Wilayah Sumatra Barat sampai Nusa Tenggara Timur”. Shalawat serta salam senantiasa tercurah kepada Nabi Muhammad SAW, manusia luar biasa karena kecerdasannya, kemuliaan akhlaqnya, keluhuran budi pekertinya, dan sunnah-sunnah Rasulullah yang tetap subur. Dalam penyusunan skripsi ini, penulis banyak mendapat dorongan, semangat, dan bimbingan serta kritikan dan sara dari berbagai pihak. Oleh karena itu, pada kesempatan ini penulis ingin mengucapkan terima kasih kepada: 1. Orang tuaku (Bapak Surais, Ibu Sulasmi dan Ibu Sumariyati), Kakakku (Mas Wahid) dan adikku (Zulfikar) tercinta serta seluruh keluarga besar penulis yang selalu memberikan kasih sayang dan selalu mendoakan penulis sehingga dapat menyelesaikan skripsi ini, 2. Dr. Agus Salim, M.Si, selaku Dekan Fakultas Sains dan Teknologi Universitas Islam Negeri Syarif Hidayatullah Jakarta,
vii
3. Ibu Yanne Irene, M.Si, selaku Ketua Program Studi Matematika Fakultas Sains dan Teknologi UIN Syarif Hidayatullah Jakarta, 4. Ibu Suma’inna, M.Si, selaku Pembimbing I dan Mahmudi, M.Si, selaku Pembimbing II. Mohon maaf atas semua kesalahan penulis selama ini, serta terima kasih banyak atas waktu dan kesabarannya dalam membimbing penulis dalam menyelesaikan penelitian ini, 5. Seluruh dosen dan karyawan Program Studi Matematika, yang telah mengajarkan banyak hal yang sangat bermanfaat bagi penulis. 6. Saudara-saudaraku, Mas Karno, Wawan, Degleng, Bono, Mas Adi, dan temanteman KremboL, 7. Shilvia ♥, yang selalu sabar dan tidak pernah berhenti memberi semangat penulis. 8. Sahabat-sahabat terbaikku selama mengenyam pendidikan di UIN Jakarta khususnya Math’08, Putra, Deny, Munir, Tedy, Ucok Lah, Azmi, Rijak, Nijaruddin, Sopyan Lampung, Ustad Luqman, dll, 9. Teman-teman Himatika UIN Jakarta. 10. Teman-teman Kosan Antalaklai (Mas Kiki dan Mbak Desi, Ka Ucal, bos Rengki, Furqon, Kochar, Dolay, Ega, Angga, Mukmin, Cahyadi, Hanip, Mali, Yopy, Azom, Zul, Onay, Awan) Bu Kosan Antalaklai, serta semua pihak yang telah membantu penulis.
viii
Penulis menyadari bahwa masih banyak kelemahan dan kekurangan yang terdapat pada skripsi ini. Atas dasar itulah penulis memohon maaf yang sebesarbesarnya kepada semua pihak jika terdapat kesalahan yang kurang berkenan. Namun, saran dan kritik selalu penulis harapkan demi perbaikan pada penelitian selanjutnya. Akhir kata, harapan yang besar bahwa skripsi ini dapat bemanfaat dan memberikan kontribusi yang berarti, baik bagi penulis khususnya dan bagi pembaca umumnya. Wassalamu’alaikum Warahmatullahi Wabarakatuh
Jakarta, 8 Januari 2015
Penulis
ix
DAFTAR ISI
HALAMAN JUDUL ........................................................................................
i
PENGESAHAN UJIAN ...................................................................................
ii
PERNYATAAN ................................................................................................ iii PERSEMBAHAN DAN MOTTO ................................................................... iv ABSTRAK ........................................................................................................
v
ABSTRACT ...................................................................................................... vi KATA PENGANTAR ...................................................................................... vii DAFTAR ISI .....................................................................................................
x
DAFTAR TABEL ............................................................................................ xii DAFTAR GAMBAR ........................................................................................ xiii DAFTAR LAMPIRAN .................................................................................... xiv BAB I PENDAHULUAN 1.1 Latar Belakang ..............................................................................
1
1.2 Rumusan Masalah .........................................................................
3
1.3 Batasan Masalah ...........................................................................
3
1.4 Tujuan Penelitian ..........................................................................
4
1.5 Manfaat Penelitian ........................................................................
4
BAB II LANDASAN TEORI 2.1 Peluang dan Variabel Acak ...........................................................
5
2.2 Distribusi Poisson .........................................................................
6
2.3 Independent Mixture Models ........................................................
9
2.4 Rantai Markov ............................................................................... 11 2.5 PHMMs (Poisson Hidden Markov Models) .................................. 17 2.6 Penskalaan Komputasi Likelihood ................................................ 22
x
2.7 Estimasi Algoritma EM (Expectaton Maximisation Algorithm) pada HMM .................................................................................... 23 2.8 Pemilihan Model Berdasarkan AIC ............................................... 27 BAB III METODE PENELITIAN 3.1 Sumber Data .................................................................................. 29 3.2 Tahap Persiapan Data ................................................................... 30 3.3 Tahap Pemodelan ........................................................................... 30 3.4 Pemilihan Model Terbaik .............................................................. 34 3.5 Alur Peneliitian .............................................................................. 35 BAB IV PEMBAHASAN 4.1 Deskripsi Data ............................................................................... 36 4.2 Pengecekan Overdispersi Data ..................................................... 38 4.3 Pemodelan Banyaknya Gempa Bumi dengan PHMMs ................ 39 BAB V KESIMPULAN DAN SARAN 5.1 Kesimpulan ................................................................................... 47 5.2 Saran ............................................................................................. 47 DAFTAR PUSTAKA ........................................................................................ 48 LAMPIRAN ....................................................................................................... 50
xi
DAFTAR TABEL Tabel 2.1 Peluang Keadaan Cuaca ..................................................................... 15 Tabel 3.1 Data Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur Tahun 2008-2013 .................................................... 29 Tabel 4.1 Data Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 SK pada Kedalaman ≤60 Km Tahun 2008-2013 ...................................................................... 36 Tabel 4.2 Banyaknya Peristiwa Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Periode 15 Hari Tahun 2008-2013 ......................................................................................... 38 Tabel 4.3
Statistik Deskriptif Banyaknya Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur Periode 15 Hari Tahun 2008-2013............................................................................... 39
Tabel 4.4 Hasil Perhitunggan Parameter 𝛌 dan Parameter 𝛅 pada 2 Keadaan Tersembunyi ...................................................................................... 41 Tabel 4.5 Peluang pada 2 Keadaan Tersembunyi .............................................. 42 Tabel 4.6 Parameter Input 𝛌, 𝛅, dan 𝚪 pada masing-masing PHMM .............. 43 Tabel 4.7 Parameter Estimasi Algoritma EM pada masing-masing PHMM ..... 45
xii
DAFTAR GAMBAR Gambar 2.1
Matriks Peluang Transisi Satu Langkah berukuran 𝑚 × 𝑚 ......... 14
Gambar 2.2
Graf Dasar HMM ......................................................................... 18
Gambar 3.1
Alur Penelitian ............................................................................. 35
Gambar 4.1
Sebaran Data Gempa Bumi dengan Magnitudo ≥5 SK pada Kedalaman ≤60 Km Tahun 2008-2013 ........................................ 37
Gambar 4.2
Peta Sebaran Gempa Bumi dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Tahun 2008-2013 ................................ 37
Gambar 4.3
Sebaran Banyaknya Peristiwa Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Periode 15 Hari Tahun 2008-2013 .......................................................................... 39
Gambar 4.4a Hasil Pengelompokan Data Banyaknya Gempa Bumi pada 2 Keadaan Tersembunyi .................................................................. 43 Gambar 4.4b Hasil Pengelompokan Data Banyaknya Gempa Bumi pada 3 Keadaan Tersembunyi .................................................................. 44 Gambar 4.4c Hasil Pengelompokan Data Banyaknya Gempa Bumi pada 4 Keadaan Tersembunyi .................................................................. 44
xiii
DAFTAR LAMPIRAN Lampiran 1. Data Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 SK pada Kedalaman ≤60 Km Tahun 2008-2013 ................................................................................ 50 Lampiran 2. Peta Sebaran Gempa Bumi .................................................................. 51 Lampiran 3. Banyaknya Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 SK pada Kedalaman ≤60 Km Periode 15 Hari Tahun 2008-2013......................................... 52 Lampiran 4. Kode Komputasi Software R untuk Peluang Forward dan Peluang Backward pada PHMM ........................................................................ 53 Lampiran 5. Kode Komputasi Software R untuk Estimasi EM pada PHMM .......... 54 Lampiran 6. Parameter Estimasi Algoritma EM untuk 𝑚 = 2................................. 55 Lampiran 7. Parameter Estimasi Algoritma EM untuk 𝑚 = 3................................. 56 Lampiran 8. Parameter Estimasi Algoritma EM untuk 𝑚 = 4................................. 57
xiv
BAB I PENDAHULUAN
1.1 Latar Belakang Indonesia merupakan negara kepulauan yang secara astronomis terletak di 60LU-110LS dan antara 950BT-1410 BT dan secara geografis teletak diantara dua benua, yaitu Benua Asia dan Benua Australia, serta diantara dua samudera, yaitu Samudera Hindia dan Samudera Pasifik. Dan secara geologis wilayah Indonesia dilalui oleh dua jalur pegunungan muda dunia, yaitu Pegunungan Mediterania di sebelah barat dan Pegunungan Sirkum Pasifik di sebelah timur, atau berada di kawasan dengan sebutan "Pacific Ring of Fire", yaitu sebuah zona dimana banyak memiliki gunung berapi yang aktif sehingga sering terjadi gempa bumi vulkanik. Selain gempa bumi vulkanik, Indonesia juga rawan akan gempa bumi tektonik, dikarenakan letak Indonesia berada pada pertemuan empat lempeng tektonik yaitu lempeng Benua Asia, lempeng Benua Australia, lempeng Samudera Hindia, dan lempeng Samudera Pasifik[10]. Sejumlah peristiwa bencana gempa bumi dengan magnitudo besar telah tejadi 10 tahun terakhir di sejumlah wilayah di Indonesia. Peristiwa gempa bumi tersebut antara lain gempa bumi di Aceh pada tanggal 26 Desember 2004 yang mengakibatkan terjadinya bencana tsunami, gempa bumi Nias pada tanggal 28 Maret 2005, gempa bumi di Yogyakarta pada tanggal 27 Mei 2006, di Pangandaran 17 Juli
1
2006, Tasikmalaya 2 September 2009 dan gempa bumi Padang Pariaman pada tanggal 30 September 2009. Tidak kurang dari 177.701 orang meninggal dunia akibat peristiwa tersebut[11]. Secara umum, karakteristik peristiwa banyaknya gempa bumi dalam periode tertentu didekati oleh distribusi Poisson, dengan sifat khas yang dimiliki yaitu ratarata dan variansinya memiliki nilai yang sama. Akan tetapi, terkadang dalam pengaplikasian distribusi Poisson, khususnya pada kasus kegempaan, terjadi overdispersi relatif dikarenakan sifat tersebut tidak sepenuhnya terpenuhi, sehingga mengakibatkan ketidaktepatan distribusi sebagai model. Terjadinya overdispersi pada distribusi tertentu diduga disebabkan adanya beberapa pengelompokan data pada distribusi tersebut, dengan masing-masing kelompok memiliki parameter yang berbeda-beda. Salah satu metode untuk mengatasi masalah pengamatan overdispersi adalah dengan mengunakan Mixture Model (model campuran)[5]. Hidden Markov Models (HMMs) adalah perluasan dari rantai Markov dimana distribusi hasil pengamatannya bergantung pada keadaan proses Markov yang tidak teramati (unobserved). Distribusi marginal dari HMMs merupakan distribusi campuran. Pada kasus kegempaan, HMMs lebih dikenal dengan istilah PHMMs (Poisson Hidden Markov Models). Dalam PHMMs, setiap observasinya dihasilkan oleh salah satu dari 𝑚 keadaan tersembunyi (hidden state)[5]. PHMMs diaplikasikan untuk mengidentifikasi pola-pola barisan keadaan yang tidak teramati yang mendasari barisan observasi.
2
Pada penelitian ini, peneliti tertarik untuk membangun sebuah model banyaknya gempa bumi menggunakan metode PHMMs serta mencari nilai estimator parameternya dengan menggunakan Alogorima EM (Expectation Maximization Algorithm). Wilayah penelitian yang diteliti adalah wilayah Propinsi Sumatra Barat hingga Propinsi Nusa Tenggara Timur Indonesia.
1.2 Rumusan Masalah Berdasarkan latar belakang yang telah dikemukakan, rumusan masalah dalam penulisan ini antara lain: 1. Bagaimana membangun model untuk mengestimasi banyaknya gempa bumi menggunakan PHMMs? 2. Apakah model terbaik untuk mengestimasi banyaknya gempa bumi di wilayah Propinsi Sumatra Barat hingga Propinsi Nusa Tenggara Timur menggunakan metode PHMMs dengan Algoritma EM? 3. Bagaimana
hasil
estimasi
parameter
banyaknya
gempa
bumi
dengan
menggunakan metode PHMMs?
1.3 Batasan Masalah Ruang lingkup pembahasan dalam skripsi ini hanya membahas banyaknya peristiwa gempa bumi yang terjadi di wilayah Propinsi Sumatra Barat hingga
3
Propinsi Nusa Tenggara Timur dengan besar magnitudo gempa bumi ≥5 Skala Richter pada kedalaman gempa dangkal yaitu ≤60 Km.
1.4 Tujuan Penulisan Tujuan penelitian ini adalah: 1. Membangun model untuk mengestimasi banyaknya gempa bumi menggunakan PHMMs, 2. Mencari model banyaknya gempa bumi terbaik untuk wilayah Propinsi Sumatra Barat hingga Propinsi Nusa Tenggara Timur menggunakan metode PHMMs dengan estimasi algoritma EM, 3. Mencari nilai estimasi parameter banyaknya gempa bumi dengan menggunakan metode PHMMs.
1.5 Manfaat Penelitian Manfaat penelitian ini diantaranya adalah: 1. Dapat diketahui nilai estimasi rata-rata banyaknya gempa bumi di wilayah Propinsi Sumatra Barat hingga Propinsi Nusa Tenggara Timur. 2. Sebagai tambahan pengetahuan tentang PHMMs dan pengaplikasian Algoritma EM untuk PHMMs. 3. Sebagai bahan referensi bagi peneliti lain yang ingin mengaplikasikan PHMMs dan Algoritma EM untuk bidang yang lain.
4
BAB II LANDASAN TEORI
2.1 Peluang dan Variabel Acak A. Teori Peluang Misalkan 𝑆 adalah ruang sampel dari suatu eksperimen acak dan 𝑨 adalah kumpulan semua titik sampel yang bisa dibentuk dari 𝑆. Peluang pada 𝑆 adalah kumpulan semua fungsi 𝑃 dengan domain 𝐴 ke daerah hasil [0,1] yang memenuhi sifat sebagai berikut[1]: 1. 0 ≤ 𝑃(𝐴) ≤ 1, ∀𝐴 ∈ 𝑨, 2. 𝑃 ∅ = 0, 3. 𝑃 𝑆 = 1 Peluang terjadinya 𝐵 bila diketahui bahwa suatu kejadian lain 𝐴 telah terjadi disebut dengan peluang bersyarat dan dilambangkan dengan 𝑃(𝐵|𝐴). Lambang 𝑃(𝐵|𝐴) dibaca “peluang terjadinya 𝐵 bila 𝐴 terjadi” atau lebih singkatnya lagi “peluang 𝐵, bila 𝐴 diketahui”. Jika 𝐴 dan 𝐵 adalah 𝑆, maka peluang bersyarat 𝐵 bila 𝐴 diketahui didefinisikan sebagai[1]: 𝑃 𝐵𝐴 =
𝑃(𝐴∩𝐵) 𝑃(𝐴)
, dengan 𝑃 𝐴 > 0
(2.1)
Jika kejadian 𝐴 dan 𝐵 keduanya dapat terjadi sekaligus, maka[1]: 𝑃 𝐴 ∩ 𝐵 = 𝑃 𝐴 .𝑃 𝐵 𝐴
5
(2.2)
Berdasarkan kaidah penggandaan umum, jika dalam suatu percobaan kejadiankejadian 𝐴1 , 𝐴2 , 𝐴3 , … dapat terjadi, maka[1]: 𝑃 𝐴1 ∩ 𝐴2 ∩ 𝐴3 ∩ … ∩ 𝐴𝑘 = 𝑃 𝐴 𝑃 𝐴2 𝐴1 𝑃 𝐴3 𝐴1 ∩ 𝐴2 … 𝑃 𝐴𝑘 𝐴1 ∩ 𝐴2 ∩ … ∩ 𝐴𝑘−1
(2.3)
B. Variabel Acak Misal 𝐸 adalah sebuah percobaan dangan ruang sampel 𝑆. Peubah acak atau variabel acak didefinisikan sebagai sebuah fungsi 𝑋 yang memetakan setiap anggota 𝑠 ∈ 𝑆 dengan sebuah bilangan real 𝑋(𝑠)[2]. Jika dalam ruang sampel mengandung titik sampel yang berhingga banyaknya atau suatu deretan anggota yang banyaknya sama dengan banyaknya bilangan bulat (terhitung), maka variabel acak yang didefinisikan pada ruang sampel tersebut adalah variabel acak diskrit[2].
2.2 Distribusi Poisson Percobaan yang menghasilkan nilai-nilai bagi suatu peubah acak 𝑋, dimana banyaknya hasil percobaan terjadi selama selang waktu tertentu atau pada daerah tertentu, disebut dengan percobaan Poisson[1]. Selang waktu yang dimaksud misalkan: semenit, sehari, seminggu, sebulan, setahun, dan sebagainya. Sedangkan daerah tertentu yang dimaksud misalkan: suatu luasan, volume, bagian bahan dan sebagainya. Dalam contoh kasus misalkan: a. Banyaknya dering telepon per jam di kantor,
6
b. banyaknya kecelakaan mobil per hari di jalan tol, c. banyaknya nasabah yang mengantri di bank dalam waktu satu jam, d. banyaknya tikus sawah per hektar, e. banyaknya bakteri dalam satu tetes air, f. banyaknya kesalahan ketik per halaman dan sebagainya. Karakteristik percobaan Poisson adalah sebagai berikut[1]: 1. Banyaknya hasil percobaan yang terjadi dalam selang waktu atau daerah tertentu, tidak bergantung pada banyaknya hasil percobaan yang terjadi pada selang waktu atau daerah lain yang terpisah. 2. Peluang terjadinya suatu hasil percobaan selama selang waktu yang singkat sekali atau dalam suatu daerah yang kecil, sebanding dengan panjang selang waktu tersebut atau besarnya daerah tersebut, dan tidak bergantung pada banyaknya hasil percobaan yang terjadi di luar selang waktu atau daerah tersebut. 3. Peluang bahwa lebih dari satu hasil percobaan akan terjadi dalam selang waktu yang singkat atau daerah yang kecil dapat diabaikan. Bilangan 𝑋 yang menyatakan banyaknya hasil percobaan dalam suatu percobaan Poisson disebut peubah acak Poisson, dan sebaran peluangnya disebut dengan distribusi Poisson[1]. Distribusi Poisson pertama kali ditemukan oleh S.D. Poisson (1781–1840), seorang ahli matematika berkebangsaan Perancis. Distribusi Poisson
7
disebut juga distribusi peristiwa yang jarang terjadi, dimana melibatkan jumlah percobaan (𝑛) yang besar dengan peluang sukses (𝑝) sangat kecil[3]. Variabel acak 𝑋 dikatakan berdistribusi Poisson dengan parameter 𝜆 > 0, dinotasikan 𝑋 ~ Poisson(𝜆), jika 𝑋 mempunyai fungsi massa peluang (PMF)[4]: 𝑒 −𝜆 𝜆 𝑥
𝑝 𝑥 =
, 𝑥 = 0, 1, 2, …
𝑥!
0 ,
(2.4)
𝑥 lainnya.
Persamaan (2.4) memenuhi syarat-syarat PMF dikarenakan[4]: (i) 𝑝 𝑥 ≥ 0 karena 𝜆 > 0 (ii) Total peluangnya adalah 1 dikarenakan: ∞
𝑝 𝑥 = 𝑥
𝑥=0
𝑒 −𝜆 𝜆𝑥 = 𝑒 −𝜆 𝑥!
∞
𝑥=0
𝜆𝑥 = 𝑒 −𝜆 𝑒 𝜆 = 1 𝑥!
Variable acak berdistribusi Poisson umumnya digunakan untuk menyatakan banyaknya peristiwa yang terjadi dalam satu satuan waktu, misalnya banyaknya kecelakaan mobil per hari di jalan tol selama bulan Oktober, banyaknya mobil yang lewat selama 10 menit di suatu ruas jalan, banyaknya antrian nasabah bank dalam waktu satu jam, dan sebagainya. Distribusi Poisson mempunyai Fungsi Pembangkit Momen (MGF)[4]: ∞
𝑀 𝑡 =𝐸 𝑒
𝑡𝑥
∞ 𝑡𝑥
=
𝑒 𝑡𝑥
𝑒 𝑝 𝑥 = 𝑥=0
𝑥=0 ∞
=𝑒
−𝜆 𝑥=0 𝑡
𝜆𝑥 𝑒 −𝜆 𝑥!
(𝜆𝑒 𝑡 )𝑥 𝑥!
= 𝑒 −𝜆 𝑒 𝜆𝑒 = 𝑒 𝜆(𝑒
𝑡 −1)
(2.5)
8
untuk setiap 𝑡 ∈ ℝ. Karena 𝑀′ 𝑡 = 𝑒 𝜆(𝑒
𝑡 −1)
𝜆𝑒 𝑡
dan 𝑀′′ 𝑡 = 𝑒 𝜆(𝑒
𝑡 −1)
𝜆𝑒 𝑡 + 𝑒 𝜆(𝑒
𝑡 −1)
𝜆𝑒 𝑡
2
maka 𝜇 = 𝑀′ 0 = 𝜆
(2.6)
dan 𝜎 2 = 𝑀′′ 0 − 𝜇 = 𝜆 + 𝜆2 − 𝜆2 = 𝜆
(2.7)
Hasil tersebut memperlihatkan bahwa pada distribusi Poisson mempunyai ciri khas bahwa mean (𝜇) dan variansinya (𝜎 2 ) memiliki nilai yang sama dengan parameternya (𝜆).
2.3 Independent Mixture Models Pada distribusi Poisson dengan fungsi massa peluang 𝑝 𝑥 = 𝑒 −𝜆 𝜆𝑥 /𝑥! memiliki sifat khas bahwa nilai variansi sama dengan nilai mean (persamaan 2.6 dan 2.7). Akan tetapi ketika nilai variansi sampel (𝑠 2 ) lebih besar dari rata-rata sampelnya (𝑥) maka akan mengakibatkan overdispersi relatif pada distribusi Poisson dan ketidaktepatan distribusi sebagai model. Salah satu metode dalam mengatasi pengamatan overdispersi yaitu dengan mengunakan Mixture Model (model campuran)[5].
9
Mixture Models dirancang untuk mengakomodasi heterogenitas yang tidak teramati (unobserved) pada populasi, dimana populasi dimungkinkan terdiri dari kelompok-kelompok yang tidak teramati[5]. Sebagai contoh misalkan 𝑋 adalah distribusi banyaknya bungkus rokok yang dibeli oleh pelanggan di supermarket. Para pelanggan tersebut dapat dibagi menjadi beberapa kelompok, misalkan dibagi menjadi dua kelompok yaitu perokok pasif dan perokok aktif.
Jika banyaknya
bungkus rokok yang dibeli oleh masing-masing kelompok pelanggan berdistribusi Poisson, maka distribusi 𝑋 belum tentu berdistribusi Poisson dan memungkinkan terjadinya overdispersi relatif pada distribusi tersebut. Anggaplah masing-masing kelompok berdistibusi Poisson memiliki rata-rata 𝜆1 dengan peluang 𝛿1 dan 𝜆2 dengan peluang 𝛿2 = 1−𝛿1 , dimana dalam menentukan nilai rata-rata tersebut didasarkan pada beberapa mekanisme acak lain yang disebut “proses parameter”. Jika proses parameter dimana serangkaian variabel acaknya saling bebas (independent) dan banyaknya juga saling bebas maka proses tersebut disebut “mixture independen” (independent mixture)[5]. Secara umum, distribusi campuran yang saling bebas (independent mixture distribution) terdiri dari sejumlah bilangan/komponen yang berhingga (𝑚 komponen) dan distribusi pencampuran yang dipilih dari komponen tersebut. Misalkan 𝛿1 , … , 𝛿𝑚 adalah peluang pada masing-masing 𝑚 komponen dan misalkan 𝑝 1 , … , 𝑝 𝑚 adalah peluang fungsi densitas masing-masing komponen tersebut. Misalkan 𝑋
10
menyatakan variabel acak yang memiliki disribusi campuran. Peluang fungsi densitas 𝑋 didefinisikan sebagai[5]: 𝑚
𝑝 𝑥 =
𝛿𝑖 𝑝𝑖 𝑥
(2.8)
𝑖=1
Sehingga, untuk kasus diskrit, peluang fungsi densitas 𝑋 adalah[5]: 𝑚
Pr 𝑋 = 𝑥 =
Pr 𝑋 = 𝑥 𝐶 = 𝑖 Pr 𝐶 = 𝑖
(2.9)
𝑖=1
Misalkan 𝑌𝑖 adalah komponen variabel acak dengan peluang fungsi densitas 𝑝 𝑥 . Nilai harapan distribusi campuran X didefinisikan sebagai[5]: 𝑚
E 𝑋 =
𝑚
Pr 𝐶 = 𝑖 E 𝑋 𝐶 = 𝑖 = 𝑖=1
𝛿𝑖 E 𝑌𝑖
(2.10)
𝑖=1
Secara umum, nilai harapan untuk distribusi campuran X dengan 𝑘 moment adalah[5]: 𝑚
𝛿𝑖 E 𝑌𝑖 𝑘 , 𝑘 ∈ ℕ
E 𝑋𝑘 =
(2.11)
𝑖=1
dimana nilai variansi campuran untuk kasus dua komponen adalah: Var 𝑋 = 𝛿1 Var 𝑌1 + 𝛿2 Var 𝑌2 + 𝛿1 𝛿2 (E 𝑌1 − E 𝑌2 )2
(2.12)
2.4 Rantai Markov A. Pengantar Stokastik Proses stokastik didefinisikan sebagai proses menyusun dan mengindeks sekumpulan variabel acak { 𝐶𝑡 , 𝑡 ≥ 0}, dengan indeks 𝑡 berada pada sekumpulan 11
𝑇[6]. Sehingga 𝑇 dinamakan ruang parameter atau ruang indeks ∀ 𝑡 ∈ 𝑇, dengan 𝑡 merupakan parameter bilangan bulat tak negatif yang merepresentasikan karakteristik terukur yang kita perhatikan pada waktu 𝑡. Himpunan variabel acak 𝐶𝑡 pada proses stokastik menggambarkan keadaaan (state) sistem pada waktu 𝑡. State merupakan posisi atau keadaan yang akan ditentukan klasifikasinya. Misalkan {𝐶𝑡 , 𝑡 = 0,1,2, … , 𝑇} adalah barisan peubah acak dengan ruang keadaan himpunan bilangan berhingga (finite). Jika 𝐶𝑡 = 𝑖, 𝑖 ∈ ℤ+, maka kita katakan bahwa proses tersebut pada waktu 𝑡 berada pada keadaan 𝑖.
B. Definisi Proses Markov Pada tahun 1906, A.A. Markov seorang ahli matematika dari Rusia yang merupakan murid Chebysev mengemukakan teori ketergantungan variabel acak dari proses acak yang dikenal dengan proses Markov[7]. Proses Markov adalah proses stokastik yang mempunyai sifat bahwa jika nilai 𝐶𝑡 telah diketahui, maka 𝐶𝑠 di mana 𝑠 > 𝑡 tidak dipengaruhi oleh 𝐶𝑢 di mana 𝑢 < 𝑡. Dengan kata lain, Proses Markov merupakan fenomena dimana kejadian masa datang hanya dipengaruhi oleh masa sekarang dan tidak dipengaruhi oleh masa lalu.
C. Rantai Markov Diskrit Analisis rantai Markov (Markov Chain Analysis) merupakan suatu teknik peluang yang menganalisis pergerakan peluang dari suatu keadaan ke keadaan
12
lainnya. Proses stokastik {𝐶𝑡 ∶ 𝑡 ∈ ℕ} dikatakan sebuah rantai Markov dengan waktu diskrit jika untuk setiap 𝑡 ∈ ℕ berlaku[5]: Pr 𝐶𝑡+1 𝐶𝑡 , 𝐶𝑡−1 , … , 𝐶1 ) = Pr 𝐶𝑡+1 𝐶𝑡 )
(2.13)
Definisi tersebut dapat diinterpretasikan bahwa untuk suatu rantai Markov, sebaran bersyarat dari sembarang keadaan yang akan datang {𝐶𝑡+1 } dengan syarat keadaan pada masa lalu adalah {𝐶𝑡−1 , … , 𝐶1 } dan keadaan sekarang adalah {𝐶𝑡 }, adalah bebas terhadap semua keadaan yang lalu dan hanya bergantung dari keadaan sekarang (kita definisikan 𝐂 (𝑡) sebagai kejadian {𝐶1 , 𝐶2 , … , 𝐶𝑡 }). Persamaan (2.13) dapat ditulis sebagai berikut[5]: Pr 𝐶𝑡+1 𝑪(𝑡) ) = Pr 𝐶𝑡+1 𝐶𝑡 )
(2.14)
Kuantitas terpenting yang terkait dengan rantai Markov adalah peluang bersyarat yang disebut dengan peluang transisi (transition probabilities) dimana pada kasus rantai Markov homogen persamaannya adalah sebagai berikut: 𝛾𝑖𝑗 = Pr 𝐶𝑛+1 = 𝑗 𝐶𝑛 = 𝑖)
(2.15)
untuk semua 𝑖, 𝑗, 𝑛 = {1,2, … }. Nilai 𝛾𝑖𝑗 merupakan elemen dari 𝚪(1), yaitu elemen dari matriks peluang transisi satu langkah (one-step transition probabilities matrix). Kita definisikan 𝚪 1 = 𝚪. Nilai 𝛾𝑖𝑗 menyatakan peluang bahwa jika proses rantai Markov berada pada keadaan 𝑖, maka keadaan berikutnya ke keadaan 𝑗 setelah satu langkah dan memenuhi syarat-syarat berikut[6]: a. 𝛾𝑖𝑗 ≥ 0, untuk semua 𝑖, 𝑗 = {1, 2, 3, ..., 𝑚} b.
𝛾𝑖𝑗 = 1, untuk semua 𝑖 = {1, 2, 3, ..., 𝑚} 13
dimana 𝑚 menyatakan banyaknya keadaan pada rantai Markov dengan jumlah masing-masing baris pada matriks 𝚪 sama dengan 1. 1 𝛾11 𝛾21 ⋮ 𝛾𝑚1
2 𝛾12 𝛾22 ⋮ 𝛾𝑚2
⋯ ⋯ ⋯ ⋱ ⋯
𝑚 𝛾1𝑚 𝛾2𝑚 ⋮ 𝛾𝑚𝑚
1 2 𝚪 1 = 𝚪 = ⋮ 𝑚 Gambar 2.1 Matriks Peluang Transisi Satu Langkah berukuran m × m Pada peluang transisi 𝑡 langkah, yaitu peluang bahwa suatau proses yang mula-mula berada pada keadaan 𝑖 akan berpindah pada keadaan 𝑗 setelah 𝑡 langkah didefinisikan sebagai[5]: 𝛾𝑖𝑗 (𝑡) = Pr 𝐶𝑠+𝑡 = 𝑗 𝐶𝑠 = 𝑖)
(2.16)
untuk semua 𝑖, 𝑗, 𝑠, 𝑡 ∈ ℕ. Nilai 𝛾𝑖𝑗 (𝑡) merupakan elemen matriks 𝚪(𝑡) (t-step transition probabilities matrices). Berdasarkan persamaan Chapman Kolmogorov matriks peluang transisi 𝑡 langkah dapat ditentukan dengan mengalikan matriks 𝚪 sebanyak 𝑡 kali dengan persamaan sebagai berikut[5]: 𝚪 𝑡 + 𝑢 = 𝚪 𝑡 . 𝚪 𝑢 , untuk semua 𝑡, 𝑢 ∈ ℕ
(2.17)
Persamaan Chapman Kolmogorov mengimplikasikan bahwa 𝚪 𝑡 = 𝚪 1 𝒕 , untuk semua 𝑡 ∈ ℕ Sejauh ini, semua peluang yang ditentukan merupakan peluang bersyarat. Misalnya 𝚪(𝑡) adalah peluang bahwa pada waktu ke 𝑡 proses berada pada keadaan 𝑗 dengan syarat keadaan mula-mula adalah 𝑖. Pada sebaran tidak bersyarat rantai Markov, nilai peluang sebaranya ditentukan oleh sebaran peluang pada keadaan awal
14
(initial state). Peluang tidak bersyarat (unconditional probabilities) Pr(𝐶𝑡 = 𝑗) dari rantai Markov berada pada keadaan tertentu pada waktu ke 𝑡 didefinisakan sebagai vektor baris sebagai berikut[5]: 𝐮 𝑡 = ( Pr 𝐶𝑡 = 1 , … , Pr 𝐶𝑡 = 𝑚 ), 𝑡 ∈ ℕ
(2.18)
dimana 𝐮 1 merupakan peluang pada keadaaan awal (initial distribution) rantai Markov. Sehingga, untuk mengetahui peluang keadaan pada waktu 𝑡 + 1 diperoleh persamaan[5]: 𝐮 𝑡+1 =𝐮 𝑡 𝚪
(2.19)
Sebagai contoh, misalkan peluang besok akan terjadi hujan bergantung pada keadaan cuaca hari ini, serta tidak bergantung pada keadaan cuaca sebelumnya (sifat rantai Markov). Misalkan diberikan peluang transisi satu langkah rantai Markov dengan dua keadaan, dimana keadaan jika cuaca hujan = 1 dan jika cuaca cerah = 2 dengan peluang transisinya ditunjukkan pada tabel berikut: Tabel 2.1 Peluang Keadaan Cuaca 𝑡+1
Hari ke 𝑡 Hujan Cerah
Hujan 0,9 0,6
Cerah 0,1 0,4
Tabel di atas menjelaskan bahwa peluang cuaca besok akan hujan sebesar 0,9 dengan syarat jika cuaca hari ini hujan dan sebesar 0,1 jika cuaca hari ini adalah cerah. Matriks peluang transisi 𝚪 untuk memodelkan kondisi di atas adalah:
𝚪 =
1 2
1 2 0,9 0,1 0,6 0,4
15
Jika cuaca hari ini (𝑡 = 1) adalah cerah, maka peluang tidak bersyarat cuaca hari ini adalah (lihat Persamaan (2.18)): 𝐮 1 = Pr 𝐶1 = 1 , Pr 𝐶1 = 2
= (0 , 1)
dan peluang cuaca untuk 𝑡 = {2, 3, … }, yaitu peluang cuaca pada hari-hari berikutnya adalah (lihat Persamaan (2.19)): 𝐮 2 = Pr 𝐶2 = 1 , Pr 𝐶2 = 2 𝐮 3 = Pr 𝐶3 = 1 , Pr 𝐶3 = 2
= 𝐮 1 . 𝚪 = (0,6 , 0,4)
= 𝐮 2 . 𝚪 = (0,78 , 0,22), dan seterusnya.
Jadi, peluang tidak bersyarat jika cuaca hari ini cerah maka peluang bahwa cuaca 2 hari yang akan datang akan turun hujan sebesar 0,78 dan cuaca akan cerah sebesar 0,22. Pada rantai Markov, setelah proses Markov berjalan dalam jangka yang panjang (𝑡 → ∞), peluang yang dihasilkan akan bernilai tetap, dengan kata lain akan menuju stasioner. Hal ini berarti, proses Markov memiliki peluang yang konvergen ke suatu nilai yang sama, dimana proses tersebut tidak lagi bergantung pada keadaan awal (𝒖 1 ). Rantai Markov dengan matriks peluang transisi 𝚪, dikatakan memiliki distribusi yang stasioner 𝜹 (vektor baris dengan elemen tak negatif) jika[5]: 𝑚
𝛅𝚪 = 𝛅 dan
𝛿𝑖 = 1
(2.20)
𝑖=1
Dalam contoh kasus cuaca sebelumnya, perhitungan 𝛅 yang stasioner pada rantai Markov dengan matriks peluang transisi keadaan cuaca 𝚪 diperoleh dengan menggunakan Sistem Persamaan Linear (SPL) sebagai berikut:
16
𝛿1 = 0,9𝛿1 + 0,6𝛿2 dan 𝛿1 + 𝛿2 = 1 sehingga diperoleh 𝛿1 =
6 7
dan 𝛿2 =
1 7
atau dapat ditulis 𝛅 =
jangka panjang, peluang cuaca cerah akan konvergen ke
1 7
1 7
(6 , 1). Artinya, dalam
dan peluang cuaca hujan
6
konvergen ke . 7
2.5 Poisson Hidden Markov Models (PHMMs) Poisson Hidden Markov Models (PHMMs) merupakan perluasan dari HMMs (Hidden Markov Models), dimana setiap observasinya dihasilkan oleh salah satu dari 𝑚 keadaan tersembunyi (hidden state) yang berdistribusi Poisson. Dalam memahami PHMMs, berikut dijelaskan mengenai dasar teori Hidden Markov Model (HMM), distribusi marginal dan moment HMM, serta likelihood HMM. 1. Hidden Markov Model (HMM) Hidden Markov Model merupakan sebuah proses stokastik yang terdiri dari dua bagian. Bagaian pertama adalah bagian yang tidak teramati (unobserved) {𝐶𝑡 , 𝑡 ∈ ℕ} yang memenuhi sifat Markov, yang disebut juga sebagai “proses parameter”
(parameter
process).
Bagian
kedua
adalah
bagian
yang
teramati/terobservasi {𝑋𝑡 , 𝑡 ∈ ℕ}, yang disebut juga sebagai “proses keadaan bergantung” (state-dependent process) dimana distribusi 𝑋𝑡 hanya bergantung pada kondisi saat 𝐶𝑡 dan bukan bergantung pada keadaan terobservasi sebelumnya 𝑋𝑡−1 . Representasi HMMs ditunjukkan pada Gambar 2.7[5].
17
Gambar 2.2 Graf Dasar HMM Hidden Markov Model {𝑋𝑡 , 𝑡 ∈ ℕ} merupakan distribusi campuran yang bergantung, dengan 𝐗 (𝑡) dan 𝐂 (𝑡) mewakili kejadian masa lalu dari waktu 1 hingga waktu 𝑡, yang disimpulkan model sederhana tersebut dengan persamaan[5]: Pr 𝐶𝑡 𝐂 (𝑡−1) ) = Pr 𝐶𝑡 𝐶𝑡−1 ), 𝑡 = 2, 3, …
(2.21)
Pr 𝑋𝑡 𝐗 (𝑡−1) , 𝐂 (𝑡) ) = Pr 𝑋𝑡 𝐶𝑡 ), 𝑡 ∈ ℕ
(2.22)
Jika rantai Markov {𝐶𝑡 } mempunyai 𝑚 keadaaan tersembunyi, kita katakan {𝑋𝑡 } adalah HMM dengan 𝑚 keadaaan. Pada kasus pengamatan diskrit, jika rantai Markov pada waktu 𝑡 berada pada keadaan 𝑖 maka fungsi masa peluang (pmf) 𝑋𝑡 didefinisikan sebagai[5]: 𝑝𝑖 𝑥 = Pr 𝑋𝑡 = 𝑥 𝐶𝑡 = 𝑖), 𝑖 = 1, 2, … , 𝑚
(2.23)
Distribusi 𝑝𝑖 dengan 𝑚 keadaan tersembunyi dapat dikatakan sebagai distribusi keadaan bergantung (state-dependent distributions). 2. Distribusi Marginal dan Moment HMM Didefinisikan 𝑢𝑖 𝑡 = Pr(𝐶𝑡 = 𝑖) untuk 𝑡 = 1, … , 𝑇. Distribusi univariat observasi 𝑋𝑡 untuk nilai diskrit sebagai berikut[5]:
18
𝑚
Pr 𝑋𝑡 = 𝑥 =
Pr(𝐶𝑡 = 𝑖)Pr 𝑋𝑡 = 𝑥 𝐶𝑡 = 𝑖) 𝑖=1 𝑚
=
𝑢𝑖 𝑡 𝑝𝑖 𝑥
(2.24)
𝑖=1
Persamaan (2.21) dapat ditulis kedalam bentuk matriks sebagai berikut: Pr 𝑋𝑡 = 𝑥 = 𝑢1 𝑡 , … , 𝑢𝑚 𝑡
𝑝1 𝑥 ⋮ 0
⋱
0 ⋮
𝑝𝑚 𝑥
1 ⋮ 1
= 𝐮 𝑡 𝐏 𝑥 𝟏′ dimana 𝐏 𝑥 didefinisikan sebagai matrik diagonal utama dengan elemen diagonal utama ke 𝑖 adalah 𝑝𝑖 𝑥 . Berdasarkan Persamaan (2.19) bahwa 𝐮 𝑡 = 𝐮 1 𝚪 𝑡−1 diperoleh: Pr 𝑋𝑡 = 𝑥 = 𝑢 1 𝚪 𝑡−1 𝐏 𝑥 𝟏′
(2.25)
Persamaan (2.25) berlaku jika rantai Markov adalah homogen dan tidak harus stasioner. Jika kita asumsikan rantai Markov stasioner dengan distribusi stasioner 𝛅, dimana 𝛅𝚪 𝑡−1 = 𝜹 maka[5]: Pr 𝑋𝑡 = 𝑥 = 𝛅𝐏 𝑥 𝟏′
(2.26)
Pada kasus distribusi bivariat, dimana terdapat empat variabel acak 𝑋𝑡 , 𝑋𝑡+𝑘 , 𝐶𝑡 , 𝐶𝑡+𝑘 dengan 𝑘 ∈ +ℤ, distribusi marginalnya adalah[5]: Pr 𝑋𝑡 = 𝑣, 𝑋𝑡+𝑘 = 𝑤 𝑚
𝑚
=
Pr 𝑋𝑡 = 𝑣, 𝑋𝑡+𝑘 = 𝑤, 𝐶𝑡 = 𝑖, 𝐶𝑡+𝑘 = 𝑗 𝑖=1 𝑗 =1
19
𝑚
𝑚
=
Pr (𝐶𝑡 = 𝑖) 𝑝𝑖 𝑣 Pr(𝐶𝑡+𝑘 = 𝑗|𝐶𝑡 = 𝑖) 𝑝𝑗 (𝑤) 𝑖=1 𝑗 =1 𝑚
𝑚
=
𝑢𝑖 𝑡 𝑝𝑖 𝑣 𝛾𝑖𝑗 (𝑘) 𝑝𝑗 (𝑤)
(2.27)
𝑖=1 𝑗 =1
Persamaan (2.27) dapat ditulis kedalam bentuk matriks sebagai berikut: Pr 𝑋𝑡 = 𝑣, 𝑋𝑡+𝑘 = 𝑤 = 𝐮 𝑡 𝐏 𝑣 𝚪 𝑘 𝐏 𝑤 𝟏′
(2.28)
Jika rantai Markov stasioner maka: Pr 𝑋𝑡 = 𝑣, 𝑋𝑡+𝑘 = 𝑤 = 𝛅𝐏 𝑣 𝚪 𝑘 𝐏 𝑤 𝟏′
(2.29)
Pada kasus stasioner, nilai harapan keadaan bergantung (state-dependent) E 𝑔(𝑋𝑡 ) dan E 𝑔(𝑋𝑡 , 𝑋𝑡+𝑘 ) untuk setiap fungsi 𝑔 adalah sebagai berikut: 𝑚
E 𝑔(𝑋𝑡 ) =
𝛿𝑖 E 𝑔(𝑋𝑡 ) 𝐶𝑡 = 𝑖
(2.30)
𝑖=1
dan 𝑚
E 𝑔(𝑋𝑡 , 𝑋𝑡+𝑘 ) =
E 𝑔(𝑋𝑡 , 𝑋𝑡+𝑘 |𝐶𝑡 = 𝑖, 𝐶𝑡+𝑘 = 𝑗) 𝛿𝑖 𝛾𝑖𝑗 (𝑘) 𝑖,𝑗 =1 𝑚
=
E 𝑔1 𝑋𝑡 𝐶𝑡 = 𝑖
E 𝑔2 (𝑋𝑡+𝑘 )|𝐶𝑡+𝑘 = 𝑗 𝛿𝑖 𝛾𝑖𝑗 𝑘
(2.31)
𝑖,𝑗 =1
dimana 𝑔 𝑋𝑡 , 𝑋𝑡+𝑘 = 𝑔1 (𝑋𝑡 )𝑔2 (𝑋𝑡+𝑘 ) dan 𝛾𝑖𝑗 𝑘 = (𝚪 𝑘 )𝑖𝑗 , untuk setiap 𝑘 ∈ ℕ. Persamaan (2.30) dan Persamaan (2.31) berguna untuk mendapatkan nilai kovarian dan korelasi pada HMMs. Sebagai contoh, misalkan jika terdapat dua keadaan HMM dimana rantai Markov berdistribusi Poisson dan stasioner maka:
20
E 𝑋𝑡 = 𝛿1 𝜆1 + 𝛿2 𝜆2 ; Var 𝑋𝑡 = E 𝑋𝑡 + 𝛿1 𝛿2 𝜆2 − 𝜆1
2
≥ E 𝑋𝑡 ;
Cov 𝑋𝑡 , 𝑋𝑡+𝑘 = 𝛿1 𝛿2 𝜆2 − 𝜆1 2 (1 − 𝛾12 − 𝛾21 )𝑘 , untuk setiap 𝑘 ∈ ℕ. 3. Likelihood HMMs Misalkan diketahui barisan observasi 𝑥1 , 𝑥2 , … , 𝑥𝑇 yang dihasilkan oleh model. Peluang likelihood 𝐿𝑇 pada barisan obsevasi tersebut dimana diberikan 𝑚 keadaaan HMM yang memiliki distrbusi inisial/awal 𝛅 dan matriks peluang transisi 𝚪 pada rantai Markov, dan peluang keadaan bergantung (state-dependent probability) 𝑝𝑖 yang merupakan elemen dari matriks 𝐏 𝑥 sebagai berikut[5]: 𝐿𝑇 = Pr 𝐗
𝑻
=𝐱
𝑻
= 𝛅𝐏 𝑥1 𝚪𝐏 𝑥2 𝚪𝐏 𝑥3 … 𝚪𝐏(𝑥𝑇 )𝟏′
(2.32)
Jika 𝛅, yaitu distribusi 𝐶1 adalah distribusi stasioner pada rantai Markov, maka: 𝐿𝑇 = Pr 𝐗
𝑻
=𝐱
𝑻
= 𝛅𝚪𝐏 𝑥1 𝚪𝐏 𝑥2 𝚪𝐏 𝑥3 … 𝚪𝐏(𝑥𝑇 )𝟏′
Misalkan kita definisikan vektor 𝜶𝑡 , untuk 𝑡 = 1, 2, … , 𝑇 (lihat Persamaan (2.32)) dengan 𝑡
𝜶𝑡 = 𝛅𝐏 𝑥1 𝚪𝐏 𝑥2 … 𝚪𝐏 𝑥𝑡 = 𝛅𝐏 𝑥1
𝚪𝐏(𝑥𝑠 ) 𝑠=2
Sehingga kita peroleh: 𝜶1 = 𝛅𝐏 𝑥1 𝜶𝑡 = 𝜶𝑡−1 𝚪𝐏 𝑥𝑡 untuk 𝑡 = 2, 3, … , 𝑇 𝐿𝑇 = 𝜶 𝑇 𝟏′
21
dengan 𝜶0 = 𝛅 untuk kasus rantai Markov stasioner. Elemen dari 𝜶𝑡 biasanya dikenal dengan peluang forward dimana akan dijelaskan pada sub bab berikutnya[5].
2.6 Penskalaan Komputasi Likelihood (Scaling the Likelihood Computation) Pada kasus distribusi state-dependent diskrit, hasil peluang elemen 𝜶𝑡 akan semakin kecil nilainya saat 𝑡 meningkat, dan pada akhirnya dibulatkan menjadi 0, atau biasa disebut underflow. Salah satu cara dalam mengatasi masalah underflow digunakan penskalaan likelihood, yaitu menghitung nilai log dari likelihood, atau bisa disebut dengan penskalaan vektor peluang forward 𝜶𝑡 . Didefinisikan sebuah vektor, untuk 𝑡 = 0,1, … , 𝑇 adalah[5]: 𝜙𝑡 = 𝜶𝑡 /𝑤𝑡 dimana 𝑤𝑡 =
𝒊 𝛼𝑡 (𝑖)
= 𝜶𝑡 𝟏′.
Pertama kita perhatikan akibat langsung dari definisi 𝝓𝑡 dan 𝑤𝑡 : 𝑤0 = 𝜶0 𝟏′ = 𝛅𝟏′ = 1; 𝜙𝑡 = 𝛅; 𝑤𝑡 𝜙𝑡 = 𝑤𝑡−1 𝜙𝑡−1 𝐁𝑡 ;
(2.33)
𝐿 𝑇 = 𝜶 𝑇 𝟏 ′ = 𝑤𝑇 𝜙 𝑇 𝟏 ′ = 𝑤𝑇 dimana kita definisikan 𝐁𝑡 adalah sebuah matriks dengan 𝐁𝑡 = 𝚪𝐏 𝑥𝑡 . Karena 𝐿 𝑇 = 𝑤𝑇 =
𝑇 𝑡=1(𝑤𝑡 /𝑤𝑡−1 ),
berdasarkan Persamaan (2.33) diperoleh: 𝑤𝑡 = 𝑤𝑡−1 (𝝓𝑡−1 𝐁𝑡 𝟏′)
dan sehingga
22
𝑇
log 𝐿𝑇 = 𝑡=1
𝑤𝑡 log = 𝑤𝑡−1
𝑇
log 𝝓𝑡−1 𝐁𝑡 𝟏′
2.34
𝑡=1
2.7 Estimasi Algoritma EM (Expectation Maximisation Algorithm) pada HMM Salah satu metode umum yang digunakan untuk mengestimasi parameterparameter pada HMMs adalah dengan menggunakan metode Algoritma EM (Estimation Maximization Algorithm). Dalam konteks HMMs, algoritma EM dikenal sebagai algoritma Baum-Welch, dimana rantai Markov pada HMM adalah homogen dan tidak diharuskan stasioner. Parameter HMM yang diestimasi dengan Algoritma EM adalah distribusi keadaan bergantung 𝑝𝑖 , matriks peluang transisi 𝚪, dan distribusi inisial 𝛅. Dalam penerapannya, Algoritma EM membutuhkan perangakat yaitu peluang forward dan peluang backward, dimana kedua peluang tersebut dapat digunakan untuk prediksi state[5]. 1. Peluang Forward Peluang forward 𝜶𝑡 untuk 𝑡 = 1, 2, … , 𝑇 didefinisikan sebagai vektor baris[5]: 𝑡
𝜶𝑡 = 𝛅𝐏 𝑥1 𝚪𝐏 𝑥2 … 𝚪𝐏 𝑥𝑡 = 𝛅𝐏 𝑥1
𝚪𝐏(𝑥𝑠 )
(2.35)
𝑠=2
dengan 𝜹 adalah distribusi inisial rantai Markov. Berdasarkan definisi peluang forward di atas, untuk 𝑡 = 1, 2, … , 𝑇 −1, dapat ditulis 𝜶𝑡+1 = 𝜶𝑡 𝚪𝐏 𝑥𝑡+1 , atau dalam bentuk skalar: 𝑁
𝛼𝑡+1 𝑗 =
𝛼𝑡 𝑖 𝛾𝑖𝑗 𝑝𝑗 𝑥𝑡+1 , 𝑖=1
23
Artinya, 𝛼𝑡 𝑗 dimana 𝑗 adalah komponen 𝜶𝒕 adalah peluang bersama Pr(𝑋1 = 𝑥1 , 𝑋2 = 𝑥2 , … , 𝑋𝑡 = 𝑥𝑡 , 𝐶𝑡 = 𝑗) Dalil Peluang Forward[5]: Untuk 𝑡 = 1, 2, … , 𝑇 dan 𝑗 = 1, 2, … , 𝑚, 𝛼𝑡 𝑗 = Pr 𝐗
𝑡
= 𝐱 𝑡 , 𝐶𝑡 = 𝑗
2. Peluang Backward Peluang backward 𝜷𝑡 untuk 𝑡 = 1, 2, … , 𝑇 didefinisikan sebagai vektor baris[5]: 𝑡
𝜷′𝑡 = 𝚪𝐏 𝑥𝑡+1 𝚪𝐏 𝑥𝑡+2 … 𝚪𝐏 𝑥𝑇 𝟏′ =
𝚪𝐏 𝑥𝑠
𝟏′
(2.36)
𝑠=𝑡+1
dimana untuk 𝑡 = 𝑇, 𝜷 𝑇 = 1. Berdasarkan definisi peluang backward di atas, untuk 𝑡 = 1, 2, … , 𝑇 −1, dapat ditulis 𝜷′𝑡 = 𝚪𝐏 𝑥𝑡+1 𝜷′𝑡+1 . Dalil Peluang Backward[5]: Untuk 𝑡 = 1, 2, … , 𝑇 − 1 dan 𝑖 = 1, 2, … , 𝑚, 𝛽𝑡 𝑖 = Pr 𝑋𝑡+1 = 𝑥𝑡+1 , 𝑋𝑡+2 = 𝑥𝑡+2 , … , 𝑋𝑇 = 𝑥𝑇 , 𝐶𝑡 = 𝑖 , dengan ketentuan bahwa Pr 𝐶𝑡 = 𝑖 > 0. Dalam notasi lebih sederhana dapat ditulis: 𝑇 𝛽𝑡 𝑖 = Pr 𝐗 𝑇𝑡+1 = 𝐱 𝑡+1 𝐶𝑡 = 𝑖 ,
dimana 𝐗 𝑏𝑎 merupakan vektor (𝑋𝑎 , 𝑋𝑎+1 , … , 𝑋𝑏 ). Dalil di atas mengidentifikasi 𝛽𝑡 (𝑖) sebagai peluang bersyarat, yaitu peluang obeservasi 𝑥𝑡+1 , … , 𝑥𝑇 dimana diberikan rantai Markov berada pada keadaan 𝑖 pada waktu 𝑡.
24
3. Peluang Forward dan Peluang Backward Gabungan antara peluang forward dan peluang backward 𝛼𝑡 dan 𝛽𝑡 dapat diterapkan untuk menghitung peluang Pr 𝐗
𝑇
=𝐱
𝑇
, 𝐶𝑡 = 𝑖 , dimana gabungan
peluang tersebut dibutuhkan dalam pengaplikasian algoritma EM pada HMMs. Dalil Peluang Forward dan Peluang Backward[5]: Untuk 𝑡 = 1, 2, … , 𝑇 dan 𝑖 = 1, 2, … , 𝑚, 𝛼𝑡 𝑖 𝛽𝑡 𝑖 = Pr 𝐗 dan akibatnya 𝜶𝒕 𝜷′𝒕 = Pr 𝐗
𝑇
=𝐱
𝑇
𝑇
=𝐱
𝑇
, 𝐶𝑡 = 𝑖 ,
= 𝐿𝑇 , untuk setiap 𝑡.
Dan dalam pengaplikasian algoritma EM pada HMMs juga dibutuhkan dua sifat berikut[5]: Dalil Peluang Forward dan Peluang Backward[5]: Petama, untuk 𝑡 = 1, 2, … , 𝑇, Pr 𝐶𝑡 = 𝑗 𝐗
𝑇
=𝐱
𝑇
= 𝛼𝑡 𝑗 𝛽𝑡 𝑗 / 𝐿𝑇 ;
(2.37)
dan yang kedua, untuk 𝑡 = 2, … , 𝑇 Pr 𝐶𝑡−1 = 𝑗, 𝐶𝑡 = 𝑘 𝐗
𝑇
=𝐱
𝑇
= 𝛼𝑡−1 𝑗 𝛾𝑗𝑘 𝑝𝑘 𝑥𝑡 𝛽𝑡 𝑘 / 𝐿𝑇 (2.38)
Pada HMM, yaitu dimana barisan keadaan rantai Markov tidak teramati, dimungkinkan terdapat data hilang (missing value) pada barisan tersebut, yang berakibat data tidak lengkap (incomplete data). Algoritma EM merupakan sebuah metode iteratif yang juga berfungsi untuk menghitung estimasi maksimum likelihood (maximum likelihood estimation) untuk data tidak lengkap, sehingga diperoleh data
25
lengkap log-likelihood[5]. Dalam setiap iterasi Algoritma EM terdapat 2 tahap, yaitu tahap Ekspektasi atau tahap E (E-step) dan tahap Maksimisasi atau tahap M (M-step). Pada kasus HMM, barisan keadaan 𝑐1 , 𝑐2 , … , 𝑐𝑇 rantai Markov dengan variabel acak nol-satu didefinisikan sebagai[5]: 𝑢𝑗 𝑡 = 1 jika dan hanya jika 𝑐𝑡 = 𝑗, (𝑡 = 1, 2, … , 𝑇) dan 𝑣𝑗𝑘 𝑡 = 1 jika dan hanya jika 𝑐𝑡−1 = 𝑗 dan 𝑐𝑡 = 𝑘 (𝑡 = 2, 3, … , 𝑇). Data lengkap log-likelihood (CDLL) HMM, yaitu dimana terdapat barisan observasi 𝑥1 , 𝑥2 , … , 𝑥𝑇 serta data hilang 𝑐1 , 𝑐2 , … , 𝑐𝑇 , adalah[5]: log Pr (𝐗
𝑻
,𝐜
𝑻
)
𝑚
𝑚
𝑚
𝑢𝑗 1 log 𝛿𝑗 +
= 𝑗 =1
𝑚
𝑇
𝑣𝑗𝑘 𝑡 𝑗 =1 𝑘=1
𝑇
log 𝛾𝑗𝑘 +
𝑡=2
𝑢𝑗 𝑡 log 𝑝𝑗 𝑥𝑡
(2.39)
𝑗 =1 𝑡=1
= Bentuk 1 + Bentuk 2 + Bentuk 3 dimana 𝜹 adalah distribusi inisial rantai Markov (distribusi 𝐶1 ) yang tidak diharuskan stasioner. Proses atau 2 tahapan Algoritma EM pada HMM adalah sebagai berikut: 1. Tahap E (E-step) Mengganti semua nilai 𝑣𝑗𝑘 dan 𝑢𝑗 𝑡 diberikan observasi 𝐱
𝑇
dengan ekpektasi bersyaratnya jika
[5]:
𝑢𝑗 𝑡 = Pr 𝐶𝑡 = 𝑗| 𝐱
𝑇
= 𝛼𝑡 𝑗 𝛽𝑡 𝑗 /𝐿𝑇 ;
(2.40)
dan 𝑣𝑗𝑘 𝑡 = Pr 𝐶𝑡−1 = 𝑗, 𝐶𝑡 = 𝑘| 𝐱
𝑇
= 𝛼𝑡−1 𝑗 𝛾𝑗𝑘 𝑝𝑘 𝑥𝑡 𝛽𝑡 𝑘 /𝐿𝑇 .
(2.41)
26
Sebagai catatan bahwa dalam perhitungan pada E-step, diperlukan peluang forward dan peluang backward pada HMM (Persamaan (2.37) dan Persamaan (2.38), dimana untuk peluang forward tidak mengasumsikan rantai Markov{𝐶𝑡 } stasioner. Akan tetapi pada peluang backward tidak terpengaruh dengan stasioneritas rantai Markov{𝐶𝑡 }[5]. 2. Tahap M (M-step) Setelah mengganti nilai 𝑣𝑗𝑘 dan 𝑢𝑗 𝑡
dengan 𝑢𝑗 𝑡
dan 𝑣𝑗𝑘 𝑡 , langkah
berikutnya adalah memaksimalkan CDLL (Persamaan (2.39)), yang berkenaan dengan tiga set parameter, yaitu distribusi inisial 𝜹, matriks peluang transisi 𝚪, dan parameter distribusi state-dependen (dalam kasus Poisson-HMM adalah 𝜆1 , 𝜆2 , … , 𝜆3 )[5]. Kedua langkah algoritma EM di atas diulang hingga mencapai kekonvergenan pada masing-masing parameter. 2.8 Pemilihan Model berdasarkan AIC (Akaike Information Criterion) Akaike Information Criterion (AIC) diperkenalkan pertama kali oleh Akaike (1974) untuk mengidentifikasikan model dari suatu dataset. Metode ini merupakan salah satu dari metode yang menerapkan pendekatan Maximum Likelihood[9]. Persamaan AIC dalam melakukan pemilihan model adalah sebagai berikut[5]: 𝐀𝐈𝐂 = −2 log 𝐿 + 2𝑝,
(2.42)
dimana log 𝐿 adalah log-likelihood pada model dan 𝑝 adalah banyaknya parameter bebas pada model. Misalakan, jika model diberikan 2 keadaan tersembunyi 𝑚 = 2,
27
maka model tersebut memiliki parameter bebas sebanyak 𝑝 = 𝑚2 + 𝑚 − 1 = 5. Hal tersebut dikarenakan parameter 𝚪 memiliki parameter bebas sebanyak 𝑚2 , parameter 𝛅 sebanyak 𝑚 − 1, dan parameter 𝛌 sebanyak 𝑚 parameter bebas.
28
BAB III METODOLOGI PENELITIAN
3.1 Sumber Data Data yang digunakan dalam penelitian ini adalah data gempa bumi Wilayah II dari mulai tanggal 4 Maret 2008 hingga 17 Desember 2013. Data yang diambil berbentuk data sekunder yang diambil dari Balai Besar Meteorologi dan Geofisika Wilayah II Ciputat dengan cakupan wilayah dari Propinsi Sumatra Barat (Sumbar) hingga Propinsi Nusa Tenggara Timur (NTT). Dari data tersebut menjelaskan waktu, sistem koordinat, kedalaman serta besarnya kekuatan pada gempa bumi yang terjadi. Berikut adalah data gempa bumi wilayah dari Propinsi Sumatra Barat hingga Propinsi Nusa Tenggara Timur: Tabel 3.1 Data Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur Tahun 2008-2013 Tanggal/Bulan/Tahun
Jam
Menit
Detik
Lintang
Bujur
Kedalaman (Km)
Magnitudo (SR)
04/03/2008
09
45
41,0
-3,26
102,06
10,00
4,2
16/03/2008
23
05
28,0
-2,95
100,80
57,00
5,4
19/03/2008
11
18
21,0
-3,09
102,15
25,00
4,4
20/03/2008
08
26
53,9
-8,33
104,25
33,00
4,4
....
....
....
....
....
....
....
....
13/12/2013
05
29
41
-6,74
102,66
30,18
4,9
13/12/2013
15
22
12
-6,82
102,59
33,05
4,3
16/12/2013
09
03
11
-6,80
102,58
34,31
5,0
29
3.2 Tahap Persiapan Data Berdasarkan tujuan penelitian ini, yaitu mencari model estimasi banyaknya gempa bumi terbaik menggunakan metode estimasi algoritma EM (Estimation Maximization algorithm) pada PHMM (Poisson Hidden Markov Models), maka data yang perlu dipersiapkan adalah data banyaknya gempa bumi. Peneliti menyaring data pada Tabel 3.1 berdasarkan kriteria kekuatan gempa bumi, yaitu gempa bumi dengan magnitudo ≥5 Skala Richter yang terjadi kedalaman gempa bumi dangkal yaitu ≤60 Km dari permukaan terjadinya gempa bumi tersebut[12]. Data hasil penyaringan, selanjutnya dihitung banyaknya gempa bumi yang terjadi dalam kurun waktu atau periode tertentu dengan kriteria yang sudah dijelaskan di atas. Dari proses pengambilan data banyaknya gempa bumi di atas, data memiliki karakteristik distribusi Poisson. Sehingga, langkah berikutnya adalah melakukan pengecekan overdispersi data terhadap distribusi Poisson, yaitu dengan membandingkan nilai rata-rata dan variansi dari data banyaknya gempa bumi. Jika nilai variansi lebih besar dibandingkan dengan nilai rata-ratanya, maka data terjadi overdispersi terhadap distribusi Poisson.
3.3 Tahap Pemodelan A. Penentuan Parameter Input Penentuan Parameter Input merupakan tahapan dalam mencari nilai parameter awal untuk masing-masing model, yaitu mencari parameter rata-rata banyaknya gempa bumi 𝜆𝑖 = 𝜆1 , … , 𝜆𝑚 dimana untuk setiap 𝜆𝑖 memliki kriteria
30
distribusi Poisson dengan peluang awal kejadiannya 𝛅 dan matriks peluang transisi 𝚪. Sebagai contoh, misalkan diketahui barisan variabel acak banyaknya peristiwa gempa bumi {𝑋𝑡 , 𝑡 = 1, … , 50 } hasil dari suatu observasi dengan periode tertentu, dimana 𝑋𝑡 memiliki karakteristik distribusi Poisson sebagai berikut: 𝑋𝑡 = {7,7,12,8,7,7,4,6,8,8,7,8,7,4,1,3,8,9,10,11,10,2,6,12,10,16,7,1, 5,6,4,11,12,7,15,10,10,2,5,5,10,13,16,11,5,6,7,14,11,14} Karena s 2 ≈ 14,01 > 𝑥 ≈ 8,1 maka terjadi overdispersi terhadap distribusi Poisson. Sehingga, diduga terdapat pengelompokan data yang tidak teramati (unobseved) pada variabel observasi, misalkan pada periode tertentu, rata-rata banyaknya peristiwa gempa bumi terjadi dengan intensitas tinggi dan sedikit. Misalkan diberikan 2 kelompok keadaan tersembunyi (𝑚 = 2) dengan kelompok tersembunyi 1 panjang interval kelas banyaknya gempa bumi dari 1 sampai 6 dan interval kelas sisanya masuk dikelompok tersembunyi 2. Berdasarkan pembagian kelompok di atas diperoleh tiga parameter sebagai berikut: 1. 𝜆1 = (4+6+4+1+3+2+6+1+5+6+4+2+5+5+5+6)/16 = 4,0625 𝜆2 = (7+7+12+8+7+7+8+8+7+8+7+8+9+10+11+10+12+10+16+7+11+ 12+7+15+10+10+10+13+16+11+7+14+11+14)/34 = 10 2. 𝛿1 =
𝛿2 =
16 50 34 50
= 0,32
= 0,68
31
3. 𝚪 =
𝜆1 𝜆2
𝜆1
𝜆2
9/15 6/33
6/15 27/33
𝜆1 =
𝜆1 𝜆2
𝜆2
0,6 0,4 0,18 0,82
Dari hasil perhitungan di atas, dipeoleh nilai rata-rata kejadian gempa bumi pada kelompok 1 sebesar 4,0625 dan pada kelompok 2 sebesar 10 kejadian, atau dapat dinotasikan dalam bentuk vektor baris 𝛌 = (4,0625 , 10), dengan peluang inisial/awal kejadiannya 𝛅 =(0,32 , 0,68) dan mariks peluang transisi kelompok keadaan tersembunyi 𝚪. Berlaku juga untuk model dengan 𝑚 = {3,4,5} dalam mendapatkan ketiga parameter tersebut. Ketiga perameter, yaitu 𝛌, 𝛅, dan 𝚪 pada masing-masing model kemudian dimasukkan kedalam software R versi 2.12.0 sebagai parameter input untuk mendapatkan model estimasi parameter tersebut.
B. Estimasi Parameter PHMMs dengan Algoritma EM Salah satu metode umum yang digunakan dalam mengestimasi HMMs adalah dengan menggunakan metode Algoritma EM (Estimation Maximization algorithm). Perangakat yang diperlukan sebelum menerapkan Algoritma EM adalah menghitung nilai peluang forward dan peluang backward yang ditunjukkan pada Persamaan (2.35) dan Persamaan (2.36). Masing-masing hasil dari peluang forward dan peluang backward tersebut kemudian dilakukan penskalaan likelihood (Scaling the Likelihood Computation) untuk mengatasi
32
kasus underflow yang ditunjukkan pada Persamaan (2.34). Berikut adalah algoritma dalam penerapan penskalaan likelihood pada Persamaan (2.34): set 𝜙0 ← 𝛅 dan 𝑙 ← 0 untuk 𝑡 = 1, 2, ... , 𝑇 𝐯 ← 𝜙𝑡−1 𝚪𝐏 𝑥𝑡 𝑢 ← 𝐯𝟏′ 𝑙 ← 𝑙 + log 𝑢 𝜙𝑡 = 𝐯/𝑢 return l Sebagai catatan bahwa 𝚪 dan 𝐏 𝑥𝑡 adalah matriks berukuran 𝑚 × 𝑚, 𝐯 dan 𝑢 adalah vektor dengan panjang 𝑚, dan 𝑙 adalah skalar log-likelihood yang diakumulasikan. Algoritma EM terdiri dari dua langkah yaitu Tahap E (E step) dan Tahap M (M step). Tahap E merupakan langkah dalam menghitung ekspektasi bersyarat dari data yang hilang (missing data) berdasarkan Persamaan (2.40) dan Persamaan (2.41). Hasil dari E step kemudian mengganti data pengamatan yang hilang sehingga diperoleh data lengkap log-likelihood (complete data log-likelihood). Dari data lengkap log-likelihood kemudian dimaksimalkan pada Tahap M untuk masing-masing bentuk (lihat Persamaan (2.39)). Berikut adalah solusi dalam penerapan pada Tahap M: Bentuk 1 Set → 𝛿𝑗 = 𝑢𝑗 1 /
𝑚 𝑗 =1 𝑢𝑗
1 = 𝑢𝑗 1
Bentuk 2 Set → 𝛾𝑗𝑘 = 𝑓𝑗𝑘 /
𝑚 𝑘=1 𝑓𝑗𝑘 ,
dimana 𝑓𝑗𝑘 =
𝑇 𝑡=2 𝑣𝑗𝑘
𝑡 .
33
Bentuk 3 untuk Poisson-HMM, 𝑝𝑗 𝑥 = 𝑒 −λ 𝑗 . λ𝑗 𝑥 /𝑥! Set → 0 =
𝑢𝑗 𝑡 (−1 + 𝑥𝑡 /λ𝑗 ); 𝑡
sehingga: 𝑇
Set → λ𝑗 =
𝑇
𝑢𝑗 𝑡 𝑥𝑡 𝑡=1
𝑢𝑗 𝑡 . 𝑡=1
3.4 Pemilihan Model Terbaik Kriteria dalam pemilihan model estimasi terbaik pada penelitian ini adalah berdasarkan nilai AIC (Akaike Information Criterion). Kriteria AIC didefinisikan pada Persamaan (2.42), yaitu: AIC = −2 𝑙𝑜𝑔 𝐿 + 2𝑝, dimana 𝑙𝑜𝑔 𝐿 adalah nilai log-likelihood masing-masing model dan 𝑝 adalah jumlah parameter pada model tersebut. Pada penelitian ini akan dicari 3 model estimasi dengan menggunakan metode Algoritma EM, yaitu model dengan keadaan tersembunyi 𝑚 = (2, 3, 4). Model estimasi terbaik berdasarkan kriteria AIC adalah model yang memiliki nilai AIC terkecil.
34
3.5 Alur Penelitian
Mulai
Data Tahap Persiapan Data
Penyaringan Data
Pengecekan Overdispersi
Tahap Pemodelan
Penentuan input Parameter PHMMs
Menghitung Banyaknya Gempa
Menghitung Peluang Forward dan Backward
Penskalaan Likelihood
Estimasi Parameter Algoritma EM
Kesimpulan (Model Terbaik)
Selesai Gambar 3.1 Alur Penelitian
35
BAB IV HASIL DAN PEMBAHASAN
4.1 Deskripsi Data Data dalam penelitian ini adalah data gempa bumi di sepanjang Propinsi Sumatra Barat (Sumbar) sampai Propinsi Nusa Tenggara Timur (NTT) dari tanggal 4 Maret 2008 hingga tahun 17 Desember 2013 yang kemudian disaring berdasarkan besarnya magnitudo yaitu ≥5 Skala Richter yang terjadi pada kedalaman gempa bumi dangkal yaitu ≤60 Km dari permukaan terjadinya gempa bumi tersebut[12]. Kriteria penyaringan data tersebut didasarkan pada besarnya resiko dan bahaya yang ditimbulkan oleh gempa bumi tersebut, khususnya bagi masyarakat di sepanjang wilayah yang diamati. Berikut adalah data gempa bumi dari hasil penyaring berdasarkan magnitudo dan kedalaman yang telah ditentukan di atas: Tabel 4.1 Data Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 SK pada Kedalaman ≤60 Km Tahun 2008-2013 Tanggal/Bulan/Tahun
Lintang
Bujur
Kedalaman (Km)
Magnitudo (SR)
16/3/2008
-2,95
100,80
57
5,4
31/3/2008
-3,00
100,89
20
5,2
02/4/2008
-4,26
102,64
31
6,1
03/4/2008
-4,19
102,22
20
5,0
...
...
...
...
...
28/11/2013
-0,26
98,562
51,39
5,1
10/12/2013
-57,27
101,99
11,62
5,3
16/12/2013
-6,799
102,57
34.31
5,0
36
dan berikut adalah sebaran data gempa bumi beserta peta sebarannya berdasarkan Tabel 4.1:
Gambar 4.1 Sebaran Data Gempa Bumi dengan Magnitudo ≥5 SK pada Kedalaman ≤60 Km Tahun 2008-2013
Gambar 4.2 Peta Sebaran Gempa Bumi dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Tahun 2008-2013 Pada Gambar 4.1 terlihat bahwa terjadinya peristiwa gempa bumi dengan magnitudo berskala richter besar, lebih sering diikuti dengan semakin rapatnya
37
frekuensi gempa bumi yang terjadi dalam periode yang berdekatan. Serta pada Gambar 4.2 terlihat bahwa, peristiwa gempa bumi pada wilayah yang diteliti paling sering terjadi di wilayah barat Pulau Sumatra dan begitu pula di bagian selatan Pulau Jawa. Hal ini disebabkan kedua wilayah tersebut dilewati oleh lempengan tektonik Samudera Hindia yang mengakibatkan berpotensi besar sering terjadinya gempa bumi.
4.2 Pengecekan Overdispersi Data Pengecekan overdispersi dilakukan pada data banyaknya gempa bumi dalam periode 15 hari. Cara yang dilakukan adalah menghitung banyaknya kejadian gempa bumi pada Tabel 4.1. dalam periode 15 hari, kemudian membandingakan nilai rata-rata dan variansi dari data tersebut. Jika nilai variansi banyaknya gempa bumi dalam periode 15 hari lebih besar dari pada nilai rataratanya maka terjadi overdispersi terhadap data tersebut. Diperoleh data banyaknya gempa bumi dalam periode 15 hari sebanyak 141 data, yang disajikan pada Tabel 4.2 dan untuk sebaran datanya ditunjukkan pada Gambar 4.2. Tabel 4.2 Banyaknya Peristiwa Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Periode 15 Hari Tahun 2008-2013 15 Hari Ke1 2 3
Banyaknya Gempa ≥5 SK 3 2 1
...
...
139 140 141
0 1 2
38
Gambar 4.3 Sebaran Banyaknya Peristiwa Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Periode 15 Hari Tahun 2008-2013 Dari proses pengambilan data banyaknya gempa bumi pada Tabel 4.2, karakteristik pengambilan datanya mengikuti pengambilan data yang berdistribusi Poisson. Tahap selanjutnya adalah tahap pengecekan overdispersi data banyaknya gempa bumi terhadap distribusi Poisson. Berikut adalah statistik deskriptif berdasarkan data pada Tabel 4.2: Tabel 4.3 Statistik Deskriptif Banyaknya Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur Periode 15 Hari Tahun 2008-2013 𝑁
𝑥
s2
141
2,035461
4,148733
Maksimum Minimum 14
0
Dari Tabel 4.3 terilihat bahwa nilai s2 ≈ 4,15 > 𝑥 ≈ 2,036. Sehingga, data banyaknya gempa bumi pada Tabel 4.2 terjadi overdispersi terhadap distribusi Poisson. Oleh karena itu, untuk analisis selanjutnya dapat digunakan PHMMs. 4.3 Pemodelan Banyaknya Gempa Bumi dengan PHMMs (Poisson Hidden Markov Models) Pada penelitian ini akan dicari model estimasi banyaknya gempa bumi terbaik dari tiga model. Tiga model tersebut adalah model dengan keadaan
39
tersembunyi, yaitu 𝑚 = (2,3,4). Penulis hanya memodelkan hingga keadaan tersembunyi 𝑚 = 4 dikarenakan data tidak memadai untuk keadaan tersembunyi 𝑚 ≥ 5. Metode yang digunakan adalah metode PHMMs (Poisson Hidden Markov Models) dengan estimasi Algoritma EM (Expectation Maximisation Algorithm). Berikut ini adalah langkah-langkah yang dilakukan beserta hasil pengolahannya untuk mendapatkan ketiga model estimasi banyaknya gempa bumi tersebut serta penentuan model yang terbaik. A. Penentuan Parameter Input Pada tahap ini adalah menghitung nilai 𝛌 = 𝜆1 , … , 𝜆𝑚 , yaitu parameter rata-rata banyaknya gempa bumi per-15 hari, dengan peluang awal kejadiannya 𝛅 = (𝛿1 , … , 𝛿𝑚 ) dan matriks peluang transisi keadaan tersembunyi 𝚪 berukuran 𝑚 × 𝑚. Langkah awal dalam mencari parameter-parameter tersebut adalah membuat tabel distribusi frekuensi dari data banyaknya gempa bumi (Tabel 4.2), dengan jumlah kelas pada tabel distribusi frekuensi tersebut ditentukan oleh banyaknya keadaan tersembunyi yang diberikan. Pada peneliti ini, peneliti mengambil nilai range 0 sampai 11 banyaknya gempa bumi untuk pembagian interval masing-masing kelas secara seragam. Sebagai contoh model dengan 𝑚 =2, misalkan jika diberikan 2 keadaan tersembunyi dengan rata-rata banyaknya gempa bumi 𝛌 = 𝜆1 , 𝜆2 , maka terdapat 2 kelas dengan panjang interval masingmasing kelasnya:
𝑐=
range banyaknya kelas
=
12 2
=6
40
Berdasarkan nilai 𝑐 di atas, ruang sampel banyaknya gempa bumi yang ada pada keadaan tersembunyi 1 adalah {0,1,2,3,4,5} dan sisanya masuk pada keadaan tersembunyi 2. Langkah selanjutnya adalah memasukkan data pada Tabel 4.2 ke masing-masing kelompok keadaan tersembunyi berdasarkan ruang sampelnya, sehingga diperoleh nilai parameter 𝛌 = λ1 , λ2 . Parameter 𝛅 = (𝛿1 , 𝛿2 ), yaitu peluang awal pada keadaan tersembunyi 1 dan keadaan tersembunyi 2, diperoleh dengan menghitung jumlah frekuensi pada masing-masing kelompok keadaan tersembunyi dan kemudian frekuensi dari masing-masing kelompok tersebut dibagi dengan frekuensi keseluruhan keadaan tersembunyi. Hasil perhitungan parameter 𝛌 dan 𝛅 untuk kasus 2 keadaan tersembunyi disajikan pada Tabel 4.4 berikut:
1 1 1 1 1 2
3
-
2
-
1
-
4
-
2
-
-
6
1
3
2
2
3
1
4
4
5
2
6
6 ......
KEADAAN TERSEMBUNYI
......
BANYAKNYA GEMPA ≥ 5 SK
139
0
140
1
141
2
𝑓 𝛌 𝛅
1 1 1 141
......
BANYAKNYA GEMPA KEADAAN TERSEMBUNYI 2
15 HARI KE-
......
BANYAKNYA GEMPA KEADAAN TERSEMBUNYI 1
......
Tabel 4.4 Hasil Perhitunggan Parameter 𝝀 dan Parameter 𝜹 pada 2 Keadaan Tersembunyi
0
-
1
-
2
-
133 1,699248 133/141 = 0,943262
8 7,625 133/8 = 0,056738
Berdasarkan Tabel 4.4, diperoleh nilai 𝛌 = (1,699248 , 7,625) dan nilai 𝛅 = (0,943262 , 0,056738). Artinya, dalam periode waktu 15 hari, keadaan 41
tersembunyi 1 memiliki nilai rata-rata banyaknya gempa sebesar 1,699248 kejadian dengan peluang awal kejadiannya sebesar 0,943262 dan 7,625 adalah rata-rata kejadian pada keadaan tersembunyi 2 dengan peluang awal sebesar 0,056738. Nilai parameter 𝚪, yaitu matriks peluang transisi keadaan tersembunyi, dimana elemen-elemen pada matriks tersebut diperoleh dengan cara menghitung frekuensi pada masing-masing kemungkinan perpindahan keadaan tersembunyi yang kemudian dibagi dengan jumlah total masing-masing baris keadaan tersembunyi. Pada kasus 2 keadaan tersembunyi, terdapat empat kemungkinan perpindahan keadaan tersembunyi yang ditunjukkan pada Tabel 4.5 berikut: Tabel 4.5 Peluang pada 2 Keadaan Tersembunyi 𝜆𝑖
1
2
1 2
125/132 7/8
7/132 1/8
sehingga diperoleh matriks peluang transisi pada 2 keadaan tersembunyi sebagai berikut:
𝚪 =
1
2
1
0,94697
0,05303
2
0,875
0,125
Berdasarkan matriks peluang transisi di atas dapat diartikan, jika periode ini berada pada keadaan tersembunyi 1, maka peluang 15 hari yang akan datang berada pada keadaan tersembunyi 1 adalah 0,946970. Jika periode ini berada pada
42
keadaan tersembunyi 1, maka peluang 15 hari yang akan datang berada pada keadaan tersembunyi 2 adalah 0,05303 dan seterusnya. Proses perhitungan parameter 𝛌, 𝛅, dan 𝚪, dimana diberikan keadaaan tersembunyi 𝑚 = 3 dan 𝑚 = 4 dilakukan dengan cara yang sama. Berikut adalah tabel lengkap hasil dari proses perhitungan parameter 𝛌, 𝛅, dan 𝚪 dengan keadaaan tersembunyi 𝑚 = (2,3,4), serta Gambar 4.4 yang merupakan hasil dari pengelompokan masing-masing keadaan tersembunyi: Tabel 4.6 Parameter Input 𝝀, 𝜹, dan 𝜞 pada masing-masing PHMM Model
𝑖
𝛌
𝛅
𝑚=2
1 2
1,699248 7,625
𝑚 =3
1 2 3
𝑚 =4
1 2 3 4
𝚪 1
2
3
4
0,943262 0,056738
0,94697 0,875
0,05303 0,125
-
-
1,336207 4,73913 11,5
0,822695 0,163121 0,014184
0,826087 0,826087 0,5
0,156522 0,173913 0,5
0,017391 0 0
-
0,924731 3,5
0,659575 0,283688
0,68478 0,625
0,26087 0,325
0,03261 0,05
0,02174 0
6,333333 11,5
0,042553 0,014184
0,66667 0,5
0,33333 0
0 0,5
0 0
Gambar 4.4a Hasil Pengelompokan Data Banyaknya Gempa Bumi pada 2 Keadaan Tersembunyi
43
Gambar 4.4b Hasil Pengelompokan Data Banyaknya Gempa Bumi pada 3 Keadaan Tersembunyi
Gambar 4.4c Hasil Pengelompokan Data Banyaknya Gempa Bumi pada 4 Keadaan Tersembunyi B. Penaksiran Parameter-parameter PHMMs dengan Algorima EM Pada bagian ini adalah menghitung nilai estimasi parameter 𝛌, 𝛅 dan 𝚪 untuk masing-masing model dengan menggunakan algoritma EM (Expectation Maximisation Algorithm). Setelah diperoleh hasil estimasi parameter, langkah selanjutnya adalah membandingkan nilai AIC (Akaike Information Criterion) dimana nilai AIC yang terkecil adalah model estimasi banyaknya gempa bumi terbaik. Perhitungan estimasi algorima EM dikerjakan dengan menggunakan software R versi 2.12.0 dengan nilai toleransi sebesar 1𝑒 − 06. Berikut adalah hasil output estimasi parameter dengan Algoritma EM pada masing-masing keadaan tersembunyi:
44
Tabel 4.7 Parameter Estimasi Algoritma EM pada masing-masing PHMM Model
𝑝
−𝑙
AIC
𝑚=2
5
268,7069
547,4138
𝑚=3
11
𝑚=4
19
257,7146
255,53
537,4292
549,060
𝑖
𝛌
𝛅
1
1,611815
2
𝚪 1
2
3
4
1
0,905652
0,094348
-
-
5,597163
0
0,786729
0,213271
-
-
1
0,852232
0
0,831224
0,145735
0,023041
-
2
2,844612
1
0,152336
0,847664
0
-
3
12,31558
0
0
1
0
-
1
0,861266
0
0,854086
0,1294405
0
0,016474
2
2,257127
1
0,060919
0,3065668
0,632514
0
3
3,802863
0
0,240358
0,7596419
0
0
4
13,75019
0
0
0
1
0
Berdasarkan Tabel 4.7, yaitu hasil estimasi dengan menggunakan Algoritma EM pada PHMMs, terlihat bahwa nilai AIC terkecil berada pada saat diberikan 3 keadaaan tersembunyi yaitu sebesar 537,429. Sehingga dapat dikatakan bahwa model dengan 3 keadaaan tersembunyi merupakan model terbaik dari pada model dengan 𝑚 = 2 dan 𝑚 = 3. Hal tersebut dikarenakan pada model 𝑚 = 3, memiliki selisih antara nilai 𝑙 dengan nilai 𝑝 yang lebih besar dari pada model dengan 𝑚 = 2 dan 𝑚 = 4, yang mengakibatkan nilai AIC kecil. Berikut adalah hasil estimasi parameter dari model terbaik PHMM, yaitu model dengan 3 keadaan tersembunyi:
𝛌 = 0,8522318, 2,844612, 12,3155828 ; 𝛅 = 0,1,0 ;
𝚪=
0,8312241 0,1457353 0,02304062 ; 0,152336 0,847664 0 0 1 0
dengan nilai ekpektasi dan varansi PHMM:
45
3
E 𝑋𝑡 =
𝛿𝑖 𝜆𝑖 𝑖=1
= 𝛿1 𝜆1 + 𝛿2 𝜆2 + 𝛿3 𝜆3 = (0×0,8522318) + (1×2,844612) + (0×12,3155828) = 2,8446121 Var 𝑋𝑡 = E 𝑋𝑡 = 2,8446121 Sehingga dapat disimpulkan bahwa dari ketiga model estimasi banyaknya gempa bumi di sepanjang Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur, model dengan 3 keadaan tersembunyi merupakan model estimasi banyaknya gempa bumi terbaik dengan nilai estimasi parameter rata-rata banyaknya gempa bumi yang terjadi sebanyak 2,8446121 ≈ 3 peristiwa dalam kurun waktu 15 hari.
46
BAB V KESIMPULAN DAN SARAN
5.1 Kesimpulan Kesimpulan yang dapat diambil berdasarkan hasil dan pembahasan sebelumnya adalah sebagai berikut: 1. Model penentuan rata-rata banyaknya gempa bumi dapat menggunakan metode PHMM (Poisson Hidden Markov Model) dengan estimasi Algoritma EM (Expectation Maximisation Algorithm). 2. Model terbaik banyaknya gempa bumi pada penelitian ini adalah model dengan 3 keadaaan tersembumyi dengan nilai estimasi parameter rata-rata banyaknya gempa adalah 2,8446121 peristiwa dalam kurun waktu 15 hari.
5.2 Saran Dalam penelitian ini, peneliti hanya melakukan analisis penentuan mendapatkan model banyaknya gempa bumi terbaik berdasarkan kriteria pengujian. Oleh sebab itu, penulis mempunyai saran untuk peneliti lain yang juga tertarik dengan materi ini: 1. Pada penelitian selanjutnya, estimasi parameter banyaknya gempa bumi terbaik dapat dilakukan dengan mempertimbangkan pembagian lokasi pengamatan berdasarkan sebaran gempa bumi yang paling sering terjadi. 2. Penelitian ini dapat dilanjutkan kembali untuk memprediksi banyaknya gempa bumi pada periode berikutnya.
47
DAFTAR PUSTAKA
[1] Walpole, Ronald E. 1993, Pengantar Statistika, Edisi ketiga, Jakarta: Gramedia Pustaka Utama [2] Herrhyanto Nar, Gantini Tuti, 2009, Pengantar Statistika Matematis, Bandung: Yrama Widya. [3] J. Supranto, 2009, Statistik Teori dan Aplikasi, Edisi ketujuh Jilid 2, Jakarta: Erlangga. [4] Nurhayati, Nunung, Teori Peluang (PAM 2231)-Unsoed. [5] Zucchini, Walter dan Iain L, MacDonald, 2009, Hidden Markov Models for Time Series, London: CRC Press. [6] Hillier, F. S., dan Lieberman, G. J. 2008, Introduction to Operation Research, 8th Edition Jilid 2, Yogyakarta: Andi. [7] R. Rabiner, Lawrence, A Tutorial on
Hidden Markov Models and Selected
Applications in Speech Recognition, IEEE, Vol. 77, No. 2, Februari, 1989. [8] Taylor Howard, Karlin Samuel, 1984, An Introduction to Stochastic Modeling, London: Academic press. [9] Agusta, Yudi, Mixture Modelling Menggunakan Prinsip Minimum Message Length, Jurnal Sistem dan Informatika Vol. 1 (Agustus 2005), 1-16, Denpasar: STIKOM Bali.
48
[10] Zonasi Gempa bumi di Indonesia, http://www.academia.edu/4517794/Zonasi_Gempa_bumi_di_Indonesia, [1/8/2014 08.57 WIB]. [11] Daftar gempa bumi di Indonesia, http://id.wikipedia.org/wiki/Daftar_gempa_bumi_di_Indonesia, [5/8/2014 08.45 WIB]. [12] Gempa Bumi, http://id.wikipedia.org/wiki/Gempa_bumi, [4/7/2014 10.45 WIB].
49
Data Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Tahun 2008-2013 Date
Magnitudo (SR)
16-03-2008 31-03-2008 02-04-2008 03-04-2008 22-04-2008 27-04-2008 27-04-2008 01-05-2008 03-05-2008 03-05-2008 18-05-2008 20-05-2008 21-05-2008 21-05-2008 21-05-2008 23-05-2008 18-06-2008 20-06-2008 24-06-2008 07-07-2008 09-07-2008 11-07-2008 12-07-2008 20-07-2008 24-07-2008 08-08-2008 08-08-2008 08-08-2008 20-08-2008 22-08-2008 26-08-2008
5,40 5,20 6,10 5,00 5,30 5,40 5,50 5,40 5,70 5,10 6,00 6,10 5,30 5,70 5,40 5,00 5,00 5,30 5,50 5,10 5,60 6,10 5,00 5,90 5,00 6,00 5,70 5,00 5,20 5,20 6,60
29-08-2008 31-08-2008 08-09-2008
5,10 5,00 5,40
⋮
⋮
22-11-2011 27-11-2011 29-11-2011 16-12-2011 19-12-2011 22-12-2011 22-12-2011 30-12-2011 05-01-2012 24-02-2012 28-03-2012 06-04-2012 06-04-2012 30-04-2012 27-05-2012 04-06-2012 19-06-2012 07-07-2012 03-09-2012 04-09-2012 04-09-2012 04-09-2012 11-09-2012 13-09-2012 13-09-2012 14-09-2012 15-09-2012 15-09-2012 18-09-2012
5,20 5,40 5,30 5,20 5,00 5,10 5,00 5,30 5,10 5,00 5,20 5,50 5,10 5,00 5,20 5,90 5,10 5,00 6,10 5,00 5,20 5,20 5,30 5,10 5,40 6,20 5,70 5,30 5,20
50
09-11-2012 10-11-2012 21-11-2012 22-11-2012 25-12-2012 30-12-2012 31-01-2013 02-02-2013 06-02-2013 13-02-2013 25-03-2013 05-04-2013 08-04-2013 16-04-2013 16-04-2013 06-05-2013 30-05-2013 03-06-2013 11-06-2013 13-06-2013 22-06-2013 06-07-2013 09-07-2013 09-07-2013 18-07-2013 24-07-2013 14-08-2013 22-08-2013 23-09-2013 11-10-2013 14-10-2013 28-11-2013 10-12-2013 16-12-2013
5,30 5,10 5,40 5,20 5,10 5,00 5,00 5,00 5,30 5,30 5,20 5,00 5,60 5,20 5,00 5,10 5,30 5,10 5,10 5,40 5,10 6,00 5,60 5,10 5,00 5,20 5,00 5,10 5,20 5,30 5,00 5,10 5,30 5,00
51
Banyaknya Peristiwa Gempa Bumi di Propinsi Sumatra Barat sampai Propinsi Nusa Tenggara Timur dengan Magnitudo ≥5 Skala Richter pada Kedalaman ≤60 Km Periode 15 Hari Tahun 2008-2013 15 HARI KE-
BANYAKNYA GEMPA ≥ 5SK
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
3 2 1 4 2 6 0 3 4 2 3 4 4 1 3 5 1 5 2 2 6 2 0 0 0 1 1 1 0 0 3
32 33 34 35 36 37 38 39 40 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
3 0 3 3 1 9 2 6 4 2 3 1 2 1 3 1 2 4 1 4 2 0 0 1 0 0 3 0 1 0 2 1 1 0
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
0 0 4 7 0 0 0 4 0 0 2 0 2 2 0 0 2 3 1 0 2 2 1 4 1 1 1 0 1 2 0 0 1 2
52
Kode Komputasi Software R untuk Peluang Forward dan Backward pada PHMM. pois.HMM.lalphabeta <- function(x,m,lambda,gamma,delta=NULL) { if(is.null(delta))delta <- solve(t(diag(m)-gamma+1),rep(1,m)) n <- length(x) lalpha <- lbeta<-matrix(NA,m,n) allprobs <- outer(x,lambda,dpois) foo <- delta*allprobs[1,] sumfoo <- sum(foo) lscale <- log(sumfoo) foo <- foo/sumfoo lalpha[,1] <- log(foo)+lscale for (i in 2:n) { foo <- foo%*%gamma*allprobs[i,] sumfoo <- sum(foo) lscale <- lscale+log(sumfoo) foo <- foo/sumfoo lalpha[,i] <- log(foo)+lscale } lbeta[,n] <- rep(0,m) foo <- rep(1/m,m) lscale <- log(m) for (i in (n-1):1) { foo <- gamma%*%(allprobs[i+1,]*foo) lbeta[,i] <- log(foo)+lscale sumfoo <- sum(foo) foo <- foo/sumfoo lscale <- lscale+log(sumfoo) } list(la=lalpha,lb=lbeta) }
53
Kode Komputasi Software R untuk Estimasi EM pada PHMM. pois.HMM.EM <- function(x,m,lambda,gamma,delta,maxiter=1000,tol=1e-6,...) { n <- length(x) lambda.next <- lambda gamma.next <- gamma delta.next <- delta for (iter in 1:maxiter) { lallprobs <- outer(x,lambda,dpois,log=TRUE) fb <- pois.HMM.lalphabeta(x,m,lambda,gamma,delta=delta) la <- fb$la lb <- fb$lb c <- max(la[,n]) llk <- c+log(sum(exp(la[,n]-c))) for (j in 1:m) { for (k in 1:m) { gamma.next[j,k] <- gamma[j,k]*sum(exp(la[j,1:(n-1)] +lallprobs[2:n,k]+lb[k,2:n]-llk)) } lambda.next[j] <- sum(exp(la[j,]+lb[j,]-llk)*x)/sum(exp(la[j,]+lb[j,]-llk)) } gamma.next <- gamma.next/apply(gamma.next,1,sum) delta.next <- exp(la[,1]+lb[,1]-llk) delta.next <- delta.next/sum(delta.next) crit <- sum(abs(lambda-lambda.next)) + sum(abs(gamma-gamma.next)) + sum(abs(delta-delta.next)) if(crit
54
MODEL 2 STATE ALGORITMA EM Iterasi
𝜸𝟏𝟐
𝜸𝟐𝟏
𝝀𝟏
𝝀𝟐
𝜹𝟏
−𝒍
1 2 30 50 100 500 750 1000
0,05224961 0,05429220 0,08365144 0,08919704 0,09333757 0,09434776 0,09434777 0,09434777
0,8657903 0,8548018 0,8061127 0,7965288 0,7887221 0,7867287 0,7867287 0,7867287
1,731595 1,726842 1,640528 1,625604 1,614518 1,611815 1,611815 1,611815
7,099168 6,928970 5,871929 5,725293 5,621742 5,597163 5,597163 5,597163
0,98671385 0,99603296 1 1 1 1 1 1
269,25270 268,97910 268,71590 268,70880 268,70700 268,70690 268,70690 268,70690
Konvergen
0,09434777
0,7867287
1,611815
5,597163
1
268,70690
55
MODEL 3 STATE ALGORITMA EM 𝜸𝟏𝟑 𝜸𝟐𝟏 𝜸𝟐𝟑
Iterasi
𝜸𝟏𝟐
1 2 30 50 100 500 750 1000
0,1599348 0,1618725 0,1519598 0,1460954 0,1457358 0,1457353 0,1457353 0,1457353
0,01206616 0,01100507 0,02309635 0,02304903 0,02304063 0,02304062 0,02304062 0,02304062
0,8145327 0,80069415 0,15926120 0,15249610 0,15233620 0,15233600 0,15233600 0,15233600
Konvergen
0,1457353
0,02304062
0,152336
𝜸𝟑𝟏
𝜸𝟑𝟐
0 0 0 0 0 0 0 0
0,1982348 0,06828712 0 0 0 0 0 0
0,8017652 0,9317129 1 1 1 1 1 1
0
0
1
𝜹𝟐
−𝒍
MODEL 3 STATE ALGORITMA EM 𝝀𝟐 𝝀𝟑 𝜹𝟏
Iterasi
𝝀𝟏
1 2 30 50 100 500 750 1000
1,4429880 1,4607280 0,8523790 0,8515037 0,8522308 0,8522318 0,8522318 0,8522318
4,2916670 4,1488420 2,8533920 2,8443633 2,8446117 2,8446121 2,8446121 2,8446121
12,165311 12,7771740 12,2878520 12,3163394 12,3155839 12,3155828 12,3155828 12,3155828
0,770955014 0,68873699 0 0 0 0 0 0
0,228785052 0,31126097 1 1 1 1 1 1
266,9192 265,0246 257,7232 257,7146 257,7146 257,7146 257,7146 257,7146
Konvergen
0,8522318
2,8446121
12,3155828
0
1
257,7146 56
MODEL 4 STATE ALGORITMA EM 𝜸𝟏𝟑 𝜸𝟏𝟒 𝜸𝟐𝟏
Iterasi
𝜸𝟏𝟐
1 2 30 50 100 500 750 1000 1500 2000
0,2845207 0,2837699 0,1548558 0,1460763 0,1396978 0,1295602 0,1294455 0,1294405 0,1294403 0,1294403
0,02411115 0,02088013 0,00000119 0 0 0 0 0 0 0
0,01433387 0,01347044 0,01747817 0,01722440 0,01699953 0,01647856 0,01647393 0,01647373 0,01647372 0,01647372
Konvergen
0,1294403
0
0,01647372
𝜸𝟐𝟑
𝜸𝟐𝟒
0,56246940 0,52169900 0,11956610 0,11024430 0,10075460 0,06141822 0,06094021 0,06091941 0,06091846 0,06091846
0,04379472 0,04394784 0,09171166 0,11694570 0,17113420 0,63201140 0,63249300 0,63251380 0,63251480 0,63251480
0 0 0 0 0 0 0 0 0 0
0,06091846
0,63251480
0
𝜸𝟒𝟐
𝜸𝟒𝟑
MODEL 4 STATE ALGORITMA EM 𝜸𝟑𝟐 𝜸𝟑𝟒 𝜸𝟒𝟏
Iterasi
𝜸𝟑𝟏
1 2 30 50 100 500 750 1000 1500 2000
0,58797890 0,52423990 0,25324450 0,27441180 0,27265320 0,23972794 0,24033186 0,24035812 0,24035932 0,24035932
0,4120211 0,4757601 0,7467555 0,7255882 0,7273468 0,7602721 0,7596681 0,7596419 0,7596407 0,7596407
0 0 0 0 0 0 0 0 0 0
0,1747956 0,0480094 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0,82520444 0,95199060 1 1 1 1 1 1 1 1
Konvergen
0,24035932
0,7596407
0
0
0
1 57
MODEL 4 STATE ALGORITMA EM Iterasi
𝝀𝟏
𝝀𝟐
𝝀𝟑
𝝀𝟒
𝜹𝟏
1 2 30 50 100 500 750 1000 1500 2000
1,0668800 1,0885540 0,7959117 0,8136253 0,8293037 0,8609924 0,8612545 0,8612659 0,8612664 0,8612664
3,2461260 3,1354070 2,5150618 2,4884172 2,4385864 2,2568877 2,2571168 2,2571267 2,2571272 2,2571272
5,3857540 5,1396950 5,3842984 5,1595332 4,7956545 3,8030547 3,8028711 3,8028632 3,8028628 3,8028628
12,5653080 0,345648072 0,620369282 0,03375290 13,2102780 0,1454705 0,82953450 0,02499426 13,8249372 0 1 0,00000001 13,8284619 0 1 0 13,8341688 0 1 0 13,7504940 0 1 0 13,7501980 0 1 0 13,7501852 0 1 0 13,7501846 0 1 0 13,7501846 0 1 0
266,9108 262,5749 255,9739 255,9359 255,8944 255,5300 255,5300 255,5300 255,5300 255,5300
Konvergen
0,8612664
2,2571272
3,8028628
13,7501846
255,5300
0
𝜹𝟐
1
𝜹𝟑
0
−𝒍
58
