ml hands-on

handson

Pada post kali ini, saya ingin membagi sebuah contoh hands-on machine learning, yang kebetulan saya juga pakai untuk mengisi materi workshop Machine Learning dan Data Science di IT Del, Sumatera Utara beberapa bulan yang lalu bersama teman kantor saya.

Hands-on ini ditulis menggunakan bahasa python dengan library tambahan antara lain numpy, matplotlib, scikit-learn, dan jupyter notebook. Library bisa diinstal manual lewat pip atau anaconda.

Data yang akan dipakai adalah data gambar digit hasil scan-an pemilu yang sebelumnya pernah dijelaskan di post ini: https://mitbal.wordpress.com/2014/10/10/pemilu-presiden-indonesia-2014-pendekatan-pembelajaran-mesin/

Untuk data bisa diambil di sini: https://www.dropbox.com/s/0nwm318mggmi3ww/11c_2000i.zip?dl=0. Silakan unduh dan ekstrak di lokal komputer masing-masing.

Untuk link ke notebook nya sendiri bisa dilihat di sini: https://github.com/mitbal/itdel/blob/master/handson.ipynb

Hands-on ini diharapkan dapat memberikan gambaran mengenai tipikal workflow machine learning. Tentu saja akan ada perbedaan untuk problem yang berbeda, misalkan data preparation yang tergantung format datanya, fase feature engineering yang problem-dependent, begitu juga dengan skema evaluasinya.

Kalau viewer-nya tidak bisa dibuka, bisa coba unduh langsung source code-nya di, nanti di komputer lokal Anda cukup panggil jupyter notebook handson.ipynb maka notebook akan terbuka di browser Anda masing-masing (asal library nya sudah terinstal semua).

Selamat bereksperimen.

 

Iklan

Whatsapp Chat Analyzer

Jadi ceritanya mahasiswa asal Indonesia yang tinggal di Stockholm menggunakan whatsapp sebagai media komunikasi untuk saling berbagi informasi hingga keluh kesah. Karena jumlah pesertanya yang memang tidak terlalu banyak (20+++), saya rasa, dibanding grup Facebook atau mailing list, group chat whatsapp memang lebih cocok untuk komunikasi yang membutuhkan respons yang relatif cepat.

Berawal dari rasa ingin tahu akan esensi dari pembicaraan di grup whatsapp tersebut, saya pun penasaran dan ingin menganalisis kontennya. Bak gayung bersambut, ternyata whatsapp menyediakan fitur untuk mengarsip pembicaraan dan mengirimkannya lewat email. Caranya cukup tahan percakapan yang ingin diarsip. Kemudian pilih ‘Email conversation’ > ‘Without media’. Gambar di bawah menunjukkan langkah-langkah tersebut. Nama dan tulisan disamarkan untuk menjaga privasi.

archive conversationwithout media

Hasilnya adalah sebuah file text. Setiap pesan chat ditulis pada satu baris dengan format “Tanggal – User: Pesan” (tanpa tanda kutip). Kode python dibawah digunakan untuk membaca file tersebut.

def read(filename):
    """ Read text file containing whatsapp chat and return the list of list of time, author, and its text
    :param filename: the filename of the chat text file
    :return: chat 2d list
    """
    chat = []
    with open(filename, 'r') as f:
        for line in f:
            lines = line.split(' - ')  # Divide between date and the rest
            if len(lines) > 1:
                lines2 = lines[1].split(': ')  # Divide between user and text
                if len(lines2) > 1:
                    speaker = lines2[0]
                    text = lines2[1]
                else:
                    speaker = ''
                    text = lines2[0]
                timestamp = lines[0]
            else:
                timestamp = ''
                speaker = ''
                text = lines[0]
            chat += [[timestamp, speaker, text]]
    return chat

Sip, setelah kita punya datanya, kita bisa geledah dan obrak-abrik isinya untuk dianalisis lebih lanjut.

Statistika dasar

Pertama-pertama. kita ingin tahu ada berapa banyak baris dalam percakapan.

date

Hmmm, ternyata ada sekitar 137 pesan per harinya. Lumayan aktif juga.

Kemudian kita ingin tahu siapa sih yang paling sering bunyi dan siapa yang paling senyap. Kita hitung total berapa banyak pesan yang ditulis per orangnya.

user frequency
Nama dan nomor telepon disamarkan untuk menjaga privasi

Wah, ternyata ada perbedaan yang cukup jomplang antara yang paling sering mengirim pesan dengan yang paling tidak sering mengirim pesan. Bahkan yang nomor 1 jumlah pesannya sekitar 2 kali dari yang nomor 2. Distribusi ini sepertinya mengikuti Zipf’s Law.

Frekuensi kata

Fitur utama dalam analisis dokumen biasanya adalah frekuensi kata. Model ini biasanya dikenal dengan nama Bag-of-Words. Kalau misalkan fitur ini dirasa kurang cukup ekspresif, biasanya ditambah lagi menjadi n-grams model, dimana kata dihitung kemunculannya bersama kata lain.

Dalam analisis dokumen juga biasanya terdapat kata-kata yang frekuensinya tinggi namun tidak memberi makna lebih kepada tulisan, seperti kata hubung dan teman-temannya. Kumpulan kata-kata tersebut sering disebut dengan istilah stopwords. Stopwords ini biasanya dihilangkan sebelum analisis agar proses lebih ringkas dan akurat. Saya memakai daftar kata stopwords dari https://sites.google.com/site/kevinbouge/stopwords-lists yang ternyata mengambil dari thesis “A Study of Stemming Effects on Information Retrieval in Bahasa Indonesia“.

Sekarang kita coba hitung frekuensi kata dalam chat minus stopwords. Histogram hasilnya bisa dilihat di bawah.

word frequency

Hmm, menarik. Ternyata kata-kata yang paling banyak keluar adalah kata sambung yang salah ketik seperti ‘yg‘ dan ‘d‘.  Kata-kata lain yang banyak muncul adalah nama panggilan ke penghuni grup yang lain. Selain itu, ternyata grup ini cukup suka tertawa, terlihat dari 3 kemunculan ekspresi tawa yaitu ‘hahaha‘, ‘wkwk‘, dan ‘haha‘. Selain itu, ternyata grup ini rajin olah raga juga, dilihat dari kata ‘badminton‘ yang muncul lebih dari 200 kali.

Emoji

if a picture paints a thousand words then why can’t I paint you“. Kurang lebih itulah sepenggal lagu If dari Bread. Penggunaan emoji atau emoticon sangat kentara dalam chatting karena satu gambarnya bisa mengekspresikan hal yang sulit jika dicoba disampaikan dengan kata-kata. Penggunaan emoji ini sangat pesat hingga menimbulkan sentimen kalau nanti kita bisa kembali ke zaman Mesir kuno dengan tulisan hieroglyph-nya. Karena umumnya penggunaan emoji ini, saya pun mencoba menghitung frekuensi dari tiap-tiap emoji.

Yang agak menyulitkan mungkin adalah emoji ini direpresentasikan dalam unicode dan di-encode dalam UTF-8 sehingga harus dibuat fungsi konversinya dulu antar keduanya. Gambar untuk plot didapat dari https://github.com/github/gemoji/ yang saya dapat link-nya dari apps.timwhitlock.info/emoji/tables/unicode. Kesulitan lain adalah ternyata Python Image Library punya sedikit untuk png dengan transparency. Apa boleh buat, gambar emoji-nya pun harus saya konversi ke jpg dengan utilitas mogrify dari imagemagick. Histogram hasilnya bisa dilihat di gambar di bawah.

emoji frequency

Hmm, menarik. Bisa dilihat, emoji teratas adalah senyum sambil berkeringat yang saya baca sebagai ekspresi dari ‘hehe’. Sisanya adalah berbagai ekspresi untuk senyum dan tertawa. Tetapi emoji ke 3 yang agak aneh. Mata. Ada apa dengan mata? Emoji ke 8, ‘see no evil monkey’ juga aneh. Apakah ada hubungannya dengan emoji ke 3?

Topik

Target berikutnya adalah Topic modelling. Lazimnya, setahu saya ini memakai teknik Latent Dirichlet Allocation (LDA). Karena saya belum pernah implementasi teknik ini sebelumnya dan kekurangan waktu untuk mencoba, saya mencoba mencari apakah ada library lain siap pakai. Saya ketemu situs startup untuk klasifikasi teks http://www.monkeylearn.com/. Di situs tersebut terdapat servis untuk melakukan klasifikasi topik dokumen secara general. Untuk free account diberikan 1000 kali panggilan API. Yah, cukuplah untuk proyek kecil-kecilan ini.

monkeylearn1monkeylearn2

Karena model-nya dibangun dari korpus Bahasa Inggris, maka teks harus diterjemahkan dulu ke bahasa tersebut. Hal ini sebenarnya bisa mengurangi akurasi, akan tetapi akurasi bukan tujuan utama yang dicari. Pilihan pertama saya, Google Translate API, ternyata tidak gratis. Ini tentunya cukup memberatkan saya yang cuma mahasiswa. Saya pun memakai layanan alternatif dari Yandex, search engine buatan Rusia.

topic

Dan inilah hasilnya. Ternyata percakapan di grup whatsapp ini dikategorikan ke ‘Entertainment & Recreation’. Hmm, bisa jadi. Yang agak aneh adalah kategori kedua, ‘Anime’, yang merupakan subkategori dari Entertainment di atas karena seingat saya tidak ada diskusi tentang anime sama sekali di grup. Hmm, inilah sulitnya memakai model buatan orang. Sulit untuk dianalisis.

Sebenarnya masih ada lagi yang mau saya coba, seperti menggunakan t-SNE untuk melihat kemiripan antar user. Namun, apa daya TTM (thesis telah memanggil).

Kode lengkap dapat dilihat di https://github.com/mitbal/wca. Cukup ikuti petunjuk di README untuk menjalankan program.

Semoga bermanfaat. Salam.

Pemilu Presiden Indonesia 2014: Pendekatan Pembelajaran Mesin

Pemilihan Umum Presiden Indonesia tahun 2014 sudah berakhir. Kalau boleh saya bilang ini pemilu presiden paling ramai yang pernah saya ikuti (walaupun saya belum ikut banyak pemilu sebenarnya, baru 2 kali sejak dapat KTP 7 tahun lalu). Orang berantem setiap hari di media sosial karena beda pilihan presiden (Facebook belum terlalu beken lima tahun lalu kayaknya). Yang bikin kurang sedap adalah banyak munculnya berita/tudingan/tuduhan yang kebenarannya tidak dapat dipertanggungjawabkan. Setelah saya pikir-pikir, mungkin hal-hal seperti ini ada baiknya juga sebagai bahan pembelajaran untuk kita semua warga Indonesia agar semakin cerdas dan dewasa, salah satunya dalam masalah politik. Nanti ke depannya, pada pemilu yang berikutnya misalnya, kita bisa memilih dan menyikap dengan lebih bijak.

Kebetulan saya kemarin bertugas menjadi salah satu petugas pelaksana pemilu luar negeri (PPLN) untuk negara Stockholm dan Latvia. Pemilu di luar negeri tentunya mempunyai nuansa yang agak berbeda dibandingkan di dalam negeri. Di luar negeri, penggunaan hak pilih dapat dilakukan melalui pos atau datang langsung ke TPS.

Jumlah pemilih di seluruh Swedia mencapai sekitar 700an. Untuk TPS, di Stockholm ada 1 TPS yang didirikan di Wisma Duta KBRI Stockholm. Alhamdulillah pemilu berjalan lancar sampai akhir. Hasil penghitungan akhir dapat dilihat di http://indonesiskaambassaden.se/hasil-pemilu/.

Persiapan dan pelaksanaannya sendiri tidak terlalu susah. Yang paling sulit mungkin adalah masalah pendataan pemilih yang datanya sudah tidak valid (sudah pindah atau ganti kewarganegaraan) atau tidak terdaftar sebagai pemilih karena tidak melapor ketika datang ke KBRI Swedia. Mungkin solusinya bisa menggunakan personnummer?

tps
TPS di Stockholm. Sepi ya, enggak kayak di Hong Kong.
20140909_120934
Perlengkapan pemilu untuk yang memilih lewat pos.

Ketatnya pemilu kali ini memunculkan wacana “kawal pemilu” untuk menjaga agar tidak terjadi kecurangan. Bak gayung bersambut, KPU ternyata juga menyediakan hasil scan formulir C1 yang berisi hasil penghitungan suara di situsnya https://pemilu2014.kpu.go.id/c1.php. Bermunculanlah karyakarya orang-orang kreatif dan inisiatif Indonesia di seantero dunia yang menyediakan sarana bagi masyarakat untuk merekapitulasi suara dari form C1. Fenomena ini dijabarkan dengan apik oleh dosen saya Pak Ruli di sini: http://theconversation.com/open-election-data-mass-interaction-indonesian-public-as-watchdog-29450.

Fenomena ini sangat positif, apalagi mungkin bagi para pakar e-government yang mendukung open data. Akan tetapi, karena bidang ilmu saya adalah Machine Learning, saya melihat upaya crowdsourcing ini dalam cahaya yang berbeda. Saya melihat ini adalah kesempatan untuk mendapatkan data dengan labelnya yang nantinya bisa dipakai untuk supervised learning.

Data adalah hal yang berlimpah ruah di era internet seperti sekarang. Akan tetapi, kebanyakan data-data yang tersedia tersebut bersifat unsupervised atau tidak ada label yang menunjukkan bahwa data tersebut termasuk ke kategori apa. Pada zaman dahulu kala, biasanya data harus diklasifikasikan manual oleh expert (Kalau contohnya pengkategorian objek pada gambar, maka manusia, asal punya akal sehat, menjadi expert-nya. Kalau misalkan data medis, maka harus dokter yang melakukan proses labelling).

Di era sekarang, biasanya dataset dibuat disusun dengan menggunakan jasa crowdsourcing berbayar seperti misalnya Amazon Mechanical Turk (walau tetap harus disupervisi). Contoh datasetnya misalkan Image-Net dan Microsoft COCO. Tapi tetap yang bisa membuat dataset adalah institusi riset yang mempunyai funding atau korporasi besar yang memiliki kapital yang besar juga (Coba bayangkan berapa biaya yang dikeluarkan Microsoft untuk mengumpulkan data hingga ratusan ribu gambar ketika mendesain algoritma deteksi skeleton untuk Kinect).

Oleh karena itu, ini adalah kesempatan emas untuk mendapatkan data tulisan tangan digit angka orang Indonesia dengan jumlah yang banyak. Dalam satu form C1 terdapat 12 angka. Jika dikalikan dengan jumlah TPS yang mencapai lebih dari 400 ribu maka total kita bisa mendapatkan data training sebanyak hampir 5 juta data. Sebagai pembanding, data handwritten digit yang sering dipakai untuk benchmarking performa algoritma machine learning baru, MNIST (walaupun ini sebenarnya hanya subset dari dataset yang lebih besar lagi ukurannya), memiliki jumlah 60.000 data untuk training dan 10.000 data untuk testing.

Berangkat dari motivasi di atas, maka saya memutuskan untuk mencoba mengekstraksi gambar tulisan angka dari gambar dan mengumpulkannya menjadi sebuah dataset kemudian melakukan eksplorasi dan eksperimen klasifikasi dengan machine learning sehingga nantinya mungkin bisa dibuat program yang otomatis akan menghitung rekap suara hanya dengan meng-upload gambar C1. Ini saya kerjakan pada liburan musim panas kemarin disambung sela-sela kuliah semester ganjil sekarang.

Akuisisi Data

Ada beberapa situs yang menyediakan jasa crowdsourcing yang bisa dipilih. Saya memilih kawalpemilu.org karena sepertinya ini yang paling lengkap (kawalpemilu.org menyelesaikan scan C1 dari 478.685 TPS yang tersebar di seluruh Indonesia dengan mengandalkan tenaga 700an relawan dalam waktu kira-kira lima hari) dan sepertinya datanya dapat dengan mudah diekstrak (halamannya cuma pakai tag table jadi saya kira cukup parsing halaman HTML saja). Saya putuskan saya harus men-download dulu semua file gambar dan rekap karena saya tidak tahu sampai kapan data tersebut akan dibuka.

Ekspektasi kedua saya ternyata salah karena situs kawalpemilu adalah satu halaman aplikasi javascript. Tetapi setelah menengok isi source code-nya, ternyata hasil rekap bisa dirikues dari webservice (yang di-host di Google Apps Engine) dalam format json. Dari source code saya juga menjadi tahu kalau scan C1 juga dapat dirikues ke server dengan cukup menyediakan nama file dengan format nomor daerah sebanyak 7 digit, nomor tps di daerah tersebut sebanyak 3 digit, dan angka yang menunjukkan halaman C1 yang diminta, dari 01 sampai 04 (yang diperlukan halaman keempat).

Kurang lebih potongan kode untuk men-download semua file adalah sebagai berikut (file scrap.py)

import requests as req

def prefix_pad(s, n):
    """ Prepend string s with 0 until it has length of n
    :param s: The string to be prepended
    :param n: The final length of prepended string
    :return: The prepended string
    """
    while len(s) < n:
        s = '0' + s
    return s

def image_filename(i, j):
    """ Prepare the KPU formatted image filename
    :param i: the area code
    :param j: the station code
    :return: The formatted filename
    """
    return prefix_pad(i, 7) + prefix_pad(j, 3) + '04'

# Loop for every area in Indonesia
start_area = 1
end_area = 80000
tps = 0
for i in xrange(start_area, end_area):
    print 'Now donwloading data from area: ', i
    r = req.get('http://kawal-pemilu.appspot.com/api/tps/' + str(i))

    if r.text == 'Error':
        continue

    json = r.json()
    for j in xrange(len(json)):
        # Check to see whether there is no error in json and also there is vote count data inside
        if json[j][6] and json[j][7] == '0':
            # Download C1 scanned image
            print 'TPS ke: ', tps
            imfname = image_filename(str(i), str(j + 1))
            imreq = req.get('http://scanc1.kpu.go.id/viewp.php?f=' + imfname + '.jpg')

            # Save to file
            f = open('scrap/' + str(i) + '_' + str(j) + '.txt', 'w')
            f.write('prabowo, jokowi, suara sah, tidak sah' + '\n')
            prabowo_count = str(json[j][2])
            jokowi_count = str(json[j][3])
            valid_count = str(json[j][4])
            invalid_count = str(json[j][5])
            f.write(prabowo_count + ',' + jokowi_count + ',' + valid_count + ',' + invalid_count)
            f.close()
            f = open('scrap/' + str(i) + '_' + str(j) + '.jpg', 'wb')
            f.write(imreq.content)
            f.close()
            tps += 1

Setelah kode tersebut dieksekusi, butuh waktu hingga 3 minggu untuk men-download semua file gambar dan rekap json. Programnya sendiri harus di-restart beberapa kali. Hasil akhirnya adalah File gambar dengan jumlah cukup banyak sekitar 417 ribu dengan ukuran total 60 gigabytes yang sementara tersimpan dengan aman di harddisk eksternal saya. Setelah melakukan bersih-bersih dengan menghapus data yang korup, sekarang kita bisa melakukan hal-hal menarik terhadap data yang kita punya.

c1
Contoh beberapa file scan form C1 yang diunduh. Hampir semua gambar memiliki ukuran, resolusi, dan kondisi yang beraneka ragam seperti terbalik atau tertukan dengan halaman lain. Hal-hal seperti ini yang mempersulit proses ekstraksi nantinya.

Ektraksi Data

Menentukan Region of Interest.

Untuk mempercepat proses komputasi dalam mengekstrak angka yang ada di dalam gambar, kita potong gambar hingga hanya menjadi region of interest-nya saja. Hal yang mempersulit adalah resolusi gambar yang berbeda-beda. Jika harus ditangani satu per satu maka kode-nya bisa meledak. Setelah melakukan investigasi singkat, kira-kira ada 3 jenis ukuran resolusi, kecil, sedang, dan tinggi. Saya buatlah aturan yang mudah-mudahan mengakomodir semua gambar.

    dim = c1.shape
    y0 = 350
    y1 = y0 + 450
    if dim[0] < 1100:
        y0 = dim[0] * 300 / 1700
        y1 = y0 + dim[0] * 400 / 1700
    elif 1700 < dim[0] < 2400:
        y0 = 350
        y1 = y0 + 450
    else:
        y0 = dim[0] * 350 / 1700
        y1 = y0 + dim[0] * 450 / 1700
    x0 = dim[1] * 19 / 24

Ekstraksi Digit

Setelah kita mendapatkan area yang mengandung rekap suara, proses ekstraksi dapat dilakukan dengan mudah dan cepat. Ada 2 pendekatan yang saya gunakan. Pertama, ekstraksi menggunakan Hough transform yang mengambil inspirasi dari sini sudokugrab.blogspot.se/2009/07/how-does-it-all-work.html. Pendekatan kedua menggunakan Harris corner detection.

Extraction based Hough transform

Pertama, gambar perlu di-threshold menggunakan Otsu threshold untuk mendapatkan edge dari gambar. Tidak perlu teknik sulit-sulit macam Canny edge detection, karena gambar relatif bebas dari noise dan tidak terlalu kompleks.

Kemudian kita panggil fungsi untuk melakukan Hough transform. Dari Hough space ini kita bisa mendapat semua garis yag ada di gambar. Dari semua garis yang didapat, kita seleksi 4 garis yang paling dekat dengan sisi terluar gambar (2 vertikal, 2 horizontal). Dari pasangan garis tadi kita hitung titik perpotongan antara garis vertikal dan horizontal untuk mendapat pojok terluar dari kotak rekap.

Setelah mendapat 4 pojok dari kotak, kita lakukan transformasi Affine. Transformasi ini bisa dalam bentuk skala, rotasi, dan translasi di mana garis paralel sebelum transformasi akan tetap paralel setelah transformasi. Kita transformasi agar orientasinya lurus dan ukurannya seragam semua 400 x 150. Setelah ditransformasi cukup dipotong kotak-kotak sebanyak 12.

ekstrak digit
Tahapan proses, dimulai dari region of interest, thresholding, hough transform, dan terakhir transformasi.

 

Extraction based Harris corner

Berbeda dengan pendekatan menggunakan Hough transform, dengan corner detection kita sudah otomatis mendapatkan titik-titik yang potensial menjadi ujung dari kotak rekap suara tanpa harus melewati pencarian garis terlebih dahulu. Selanjutnya cukup mencari titik yang terdekat dengan keempat ujung gambar. Transformasi yang sama juga dilakukan setelah itu.

Akan tetapi setelah melihat hasil empiris, maka saya lebih memilih teknik Hough transform karena file gambar hasil akhir lebih rapi dibanding dengan corner detection. Hal ini ditengarai karena corner dapat terdeteksi pada posisi-posisi yang salah dan ini mengakibatkan angka setelah ditransformasi menjadi terdistorsi.

Batch process

Setelah sukses melakukan ekstraksi pada satu gambar C1, berikutnya proses ekstraksi  akan dijalankan ke semua gambar scan C1. Niatnya  awalnya ke semua gambar, tapi saya jalankan untuk 100.000 data pertama dahulu.

Hasil akhirnya ratusan ribu gambar angka berhasil terekstrak. Tetapi ternyata banyak gambar-gambar yang rusak karena salah ekstraksi (region of interest-nya salah sehingga kotak yang berisi angka tidak berhasil ditemukan atau terpotong). Oleh karena itu masih diperlukan seleksi secara manual (walaupun sebenarnya cepat untuk dilakukan, cukup Ctrl+X Ctrl+V, tapi membosankan untuk dilakukan). Hasil sementara ada 11 kelas (0 sampai 9 tambah kelas khusus X untuk digit ratusan atau puluhan yang nilai tidak sampai situ) dengan masing-masing kelas mempunyai 2000 data. Masih ada banyak lagi data yang bisa diambil, mungkin solusinya perlu di-crowdsourcing juga.

Semua proses di atas ada dalam modul extract.py di repositori.

digit.pn
Contoh subset dataset dengan 100 data per kelas. Datanya lumayan variatif kan…

Training & Testing

Setelah kita mendapatkan data per kelas yang cukup banyak kita bisa memulai proses training. Untuk proses ini saya menggunakan framework machine learning favorit saya, Weka. Weka men-support banyak algoritma serta memiliki kemampuan untuk menyimpan dan memuat model hasil pelatihan sebelumnya. Untuk eksperimen awal ini, kita coba menggunakan data yang masih raw dan classifier nearest neighbor (di weka namanya adalah IB1). Kita coba pada dataset yang kecil terlebih dahulu, dengan ukuran 500 instance per kelas.

weka load data
500 data per kelas, masing-masing data memiliki 5000 fitur.

Dengan skema 10-fold cross validation, akurasi yang didapatkan adalah sekitar 60%. Hmm, tidak terlalu buruk. Tentunya hasil akurasi di atas bisa ditingkatkan lebih lanjut, dengan menggunakan teknik-teknik ekstraksi fitur yang lebih canggih atau menggunakan classifier yang lebih powerful dan dilakukan fine tuning hingga hanya memiliki error yang sangat kecil.

weka run
Hasil klasifikasi dalam bentuk akurasi dan confusion matrix.

Terakhir, kita buat tampilan grafik yang digunakan untuk memuat dan melakukan ekstraksi dan klasifikasi digit dengan model yang telah didapatkan. GUI dibuat dengan Tkinter.

gui
Tampilan GUI awal; setelah me-load scan lembar C1 dan model; serta setelah diproses.

Kalau dilihat-lihat hasilnya ada salah prediksi untuk angka 2 dan karakter X. Yah, lumayanlah…

Penutup

Demikian kira-kira bagaimana ilmu machine learning bisa ikut-ikutan pada hajatan pemilu 2014. Kalau kawalpemilu butuh 5 hari untuk merekap semua suara, maka dengan bantuan komputer seharusnya itu bisa selesai lebih cepat lagi. Jika ada yang mempunyai ide berbeda untuk tahap pengekstrakan agar lebih efektif, bisa langsung dengan me-fork repositori dan langsung implementasi atau langsung kirim email untuk berdiskusi.

Untuk file dataset dapat di-download di http://db.orangedox.com/ontgRUzd9EjjlxZVoD/11c_2000i.zip. Kalau perlu gunakan module prepare_input.py untuk membaca file gambar dan menjadikannya ke file CSV yang siap diproses di Weka.

Rencananya akan saya update begitu saya bisa memisahkan lebih banyak data lagi per kelasnya. Sementara kode untuk semua bisa dilihat di repository github: https://github.com/mitbal/pemilu

Post-post berikutnya kita akan membahas metode ekstraksi fitur dan klasifikasi yang lebih powerful. Mudah-mudahan bisa coba sampai metode representasi fitur yang state-of-the-art macam fitur CNN.

Semoga berguna. Salam.

Baseline Wander Removal dengan Wavelet

Pada post kali ini, kita akan melakukan proses menghilangkan baseline wander (baseline wander removal) pada sinyal EKG dengan menggunakan transformasi wavelet seperti yang ditulis oleh Sargolzaei et al pada paper dengan judul “A new robust wavelet based algorithm for baseline wandering cancellation in ECG signals“.

Latar belakang
Elektrokardiogram (EKG) adalah pembacaan sinyal elektrik jantung dengan cara menempelkan elektroda ke posisi tertentu pada tubuh dan kemudian membaca nilai perbedaan potensial listrik alias tegangan yang dihasilkan. Dari pembacaan ini kita bisa mencari tahu mengenai kondisi jantung. Akan tetapi, sinyal yang baru dibaca ini tidak lepas dari noise sehingga bisa mengganggu proses diagnosis. Sumber noise ini berasal dari berbagai macam, seperti elektroda dari elektroda itu sendiri atau jika tubuh bergerak ketika dilakukan pengukuran. Salah satu kondisi noise yang sering menyerang adalah situasi yang disebut baseline wander. Baseline wander terjadi apabila sinyal EKG tidak lurus pada sumbu x, malah naik turun. Contoh bisa dilihat di gambar 1.

bwr
Contoh sinyal yang terkena baseline wander dan hasilnya setelah dihilangkan. Gambar diambil tanpa izin dari paper yang dirujuk di atas.

Baseline wander removal (BWR) adalah salah satu tahap preprocessing pada sinyal EKG untuk menghilangkan baseline drift ini. Ada beberapa teknik yang bisa dipakai, seperti melakukan filtering. Pada post ini, kita akan mengimplementasikan salah satu teknik yang memanfaatkan transformasi wavelet. Penjelasan dasar mengenai wavelet salah satunya bisa dilihat di http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html. Kelebihan dari teknik ini adalah kita tidak perlu menentukan parameter seperti ketika menggunakan high pass filter (frequency cut-off) sehingga metode kita bekerja secara non-supervised.

Sinyal EKG yang sudah bersih dari noise semacam ini kemudian bisa digunakan untuk berbagai macam kegunaan. Untuk pembahasan mengenai kegunaan sinyal EKG untuk mengenali tipe detak arrhythmia atau diagnosis tahapan tidur, bisa merujuk ke paper ini dan ini.

Algoritma
Berikut adalah tahapan dalam algoritma ini. Pertama, kita lakukan dekomposisi kepada sinyal asli dengan transformasi wavelet. Dalam hal ini kita memilih Daubechies orde 4 sebagai fungsi basisnya. Sinyal didekomposisi menjadi bagian frekuensi rendah/aproksimasi dan frekuensi tinggi/detail. Kemudian kita hitung nilai energi pada sinyal frekuensi tinggi. Kita cari kondisi dimana nilai energi pada level dekomposisi tersebut lebih rendah dari pada nilai pada level dekomposisi sebelumnya dan sesudahnya (alias lokal minima). Setelah kita temukan level tersebut, kita rekonstruksi sinyal aproksimasi dari level ini dengan membuang nilai pada sinyal frekuensi tinggi (atau dijadikan 0 semua). Sinyal hasil rekonstruksi ini kita sebut dengan baseline. Untuk menghilangkan baseline wander pada sinyal asli, maka kita kurangkan sinyal asli dengan sinyal baseline.

level
Gambar atas: plot nilai energi pada sinyal detail pada berbagai level dekomposisi. Tanda panah menunjukkan level ketika nilainya adalah lokal minima. Gambar bawah: Sinyal asli, baseline, dan sinyal asli yang sudah dikurang baseline. Lagi-lagi diambil dari paper rujukan di atas.

 

Implementasi
Pertama, kita implementasi metode konvolusi antara dua sinyal.

def conv(x, h):
    """ Perform the convolution operation between two input signals. The output signal length
    is the sum of the lenght of both input signal minus 1."""
    length = len(x) + len(h) - 1
    y = [0]*length

    for i in xrange(len(y)):
        for j in xrange(len(h)):
            if i-j >= 0 and i-j < len(x):
                y[i] += h[j] * x[i-j]

    return y

Lalu kita implementasi metode dekomposisi wavelet dan jangan lupa deklarasi koefisien basis fungsi yang dipakai. Dalam kasus ini adalah Daubechies dengan 4 koefisien. Dekomposisi dilakukan dengan cara mengkonvolusi sinyal asli dengan koefisien low pass and high pass. Sinyal keluaran dari proses ini kemudian di-downsampling menjadi setengahnya dengan cara hanya mengambil nilai pada posisi ganjil atau genap saja. Hal ini dilakukan berulang terhadap sinyal keluaran low pass (atau disebut aproksimasi) sebanyak parameter level.

c0 = (1+sqrt(3))/(4*sqrt(2))
c1 = (3+sqrt(3))/(4*sqrt(2))
c2 = (3-sqrt(3))/(4*sqrt(2))
c3 = (1-sqrt(3))/(4*sqrt(2))

def db4_dec(x, level):
""" Perform the wavelet decomposition to signal x with Daubechies order 4 basis function as many as specified level"""

    # Decomposition coefficient for low pass and high pass
    lpk = [c0, c1, c2, c3]
    hpk = [c3, -c2, c1, -c0]

    result = [[]]*(level+1)
    x_temp = x[:]
    for i in xrange(level):
        lp = conv(x_temp, lpk)
        hp = conv(x_temp, hpk)

        # Downsample both output by half
        lp_ds=[0]*(len(lp)/2)
        hp_ds=[0]*(len(hp)/2)
        for j in xrange(len(lp_ds)):
            lp_ds[j] = lp[2*j+1]
            hp_ds[j] = hp[2*j+1]

        result[level-i] = hp_ds
        x_temp = lp_ds[:]

    result[0] = lp_ds
    return result

Fungsi rekontruksi digunakan untuk mencari baseline dari sinyal asal. Rekonstruksi bekerja dengan melakukan konvolusi dengan koefisien rekonstruksi pada kedua sinyal low pass dan high pass, melakukan upsampling dengan menyelipkan 0 di setiap nilai pada sinyal dan kemudian masing-masing posisi pada sinyal saling dijumlahkan.

def db4_rec(signals, level):
    """ Perform reconstruction from a set of decomposed low pass and high pass signals as deep as specified level"""

    # Reconstruction coefficient
    lpk = [c3, c2, c1, c0]
    hpk = [-c0, c1, -c2, c3]

    cp_sig = signals[:]
    for i in xrange(level):
        lp = cp_sig[0]
        hp = cp_sig[1]

        # Verify new length
        length = 0
        if len(lp) > len(hp):
            length = 2*len(hp)
        else:
            length = 2*len(lp)

        # Upsampling by 2
        lpu = [0]*(length+1)
        hpu = [0]*(length+1)
        index = 0
        for j in xrange(length+1):
            if j%2 != 0:
                lpu[j] = lp[index]
                hpu[j] = hp[index]
                index += 1

        # Convolve with reconstruction coefficient
        lpc = conv(lpu, lpk)
        hpc = conv(hpu, hpk)

        # Truncate the convolved output by the length of filter kernel minus 1 at both end of the signal
        lpt = lpc[3:-3]
        hpt = hpc[3:-3]

        # Add both signals
        org = [0]*len(lpt)
        for j in xrange(len(org)):
            org[j] = lpt[j] + hpt[j]

        if len(cp_sig) > 2:
            cp_sig = [org]+cp_sig[2:]
        else:
            cp_sig = [org]

    return cp_sig[0]

Method calcEnergy menghitung nilai energi dari sebuah sinyal berdasarkan definisinya yaitu jumlahan kuadrat sinyal di tiap titik

def calcEnergy(x):
    """ Calculate the energy of a signal which is the sum of square of each points in the signal."""
    total = 0
    for i in x:
        total += i*i
    return total

Kemudian method bwr berikut adalah implementasi dari algoritma yang dijabarkan di atas.

def bwr(raw):
    """ Perform the baseline wander removal process against signal raw. The output of this method is signal with correct baseline
    and its baseline """
    en0 = 0
    en1 = 0
    en2 = 0
    n = 0

    curlp = raw[:]
    num_dec = 0
    last_lp = []
    while True:
        print 'Iterasi ke' + str(num_dec+1)
        print len(curlp)

        # Decompose 1 level
        [lp, hp] = db4_dec(curlp,1)

        # Shift and calculate the energy of detail/high pass coefficient
        en0 = en1
        en1 = en2
        en2 = calcEnergy(hp)
        print en2

        # Check if we are in the local minimum of energy function of high-pass signal
        if en0 > en1 and en1 < en2:
            last_lp = curlp
            break

        curlp = lp[:]
        num_dec = num_dec+1

    # Reconstruct the baseline from this level low pass signal up to the original length
    base = last_lp[:]
    for i in xrange(num_dec):
        base = db4_rec([base,[0]*len(base)], 1)

    # Correct the original signal by subtract it with its baseline
    ecg_out = [0]*len(raw)
    for i in xrange(len(raw)):
        ecg_out[i] =  raw[i] - base[i]

    return (base, ecg_out)

Contoh
Untuk contoh kita ambil data dari situs Physionet khususnya MIT-BIH Arrhythmia database. Data diambil dari salah satu record pasien dengan kode nomor 101. Sinyal EKG diambil sepanjang 1 menit. Sinyal EKG yang diambil berasal dari lead II dan V5.

Contoh kode untuk memanggil modul dan fungsi yang sudah dibuat di atas adalah sebagai berikut

import bwr
import matplotlib.pyplot as plt

# Read input csv file from physionet
f = open('samples1.csv', 'r')
lines = f.readlines()
f.close()

# Discard the first two lines because of header. Takes either column 1 or 2 from each lines (different signal lead)
signal = [0]*(len(lines)-2)
for i in xrange(len(signal)):
	signal[i] = float(lines[i+2].split(',')[1])

# Call the BWR method
(baseline, ecg_out) = bwr.bwr(signal)

plt.subplot(2,1,1)
plt.plot(signal, 'b-')
plt.plot(baseline, 'r-')

plt.subplot(2,1,2)
plt.plot(ecg_out, 'b-')
plt.show()

Berikut contoh pertama. Garis merah pada plot merupakan baseline dari sinyal EKG.

samples1

Berikut adalah contoh kedua.

samples2

Contoh berikut diambil dari lead V5.

samples3

Kode lengkap dan sampel sinyal dapat dilihat di: https://github.com/mitbal/py-bwr

Semoga berguna. Salam.

UPGMA

Pada post kali ini kita akan mengimplementasi algoritma UPGMA atau Unweighted Pair Group Method with Arithmetic Mean untuk melakukan clustering.

Algoritma
UPGMA bekerja dengan prinsip yang sederhana. Pada setiap iterasi, pilih pasangan point dengan point, atau point dengan cluster, atau cluster dengan cluster, yang memiliki jarak yang terpendek. Gabung kedua pasangan ini kedalam satu cluster. Hal ini dilakukan terus menerus hingga jumlah cluster berkurang menjadi jumlah yang diinginkan. Untuk menghitung jarak antar datapoint, kita bisa menggunakan berbagai kriteria jarak. Pada implementasi ini, kita menggunakan euclidean distance untuk menghitung jarak. Sedangkan untuk menghitung jarak antar cluster, kita hitung dengan cara merata-ratakan jarak antar semua pasang point dari cluster pertama ke cluster kedua, atau dengan kata lain:

\frac{1}{|A|.|B|} \sum_{x \in A} \sum_{y \in B} d(x,y)

Teknik UPGMA ini bisa digolongkan kepada teknik hierarchical clustering karena kita bisa melihat hirarki sebuah cluster yang dibentuk dari cluster-cluster lain yang lebih kecil.

Implementasi
Kita implementasikan UPGMA dalam bahasa Python. Pertama, kita butuh struktur data untuk menyimpan point yang terdapat pada sebuah cluster.

class Node:
    def __init__(self, p):
        self.points = p
        self.right = None
        self.left = None

Fungsi UPGMA menerima dua parameter. Parameter pertama adalah set yang berisi datapoints sebanyak n dengan masing-masing point memiliki dimensi d. Parameter kedua adalah jumlah cluster yang kita inginkan.

def upgma(points, k):
    """ Cluster based on distance matrix dist using Unweighted Pair Group Method with Arithmetic Mean algorithm up to k cluster"""

    # Initialize each cluster with one point
    nodes = []
    n = len(points)
    for i in xrange(n):
        node = Node([points[i]])
        nodes = nodes + [node]

    # Iterate until the number of clusters is k
    nc = n
    while nc > k:
        # Calculate the pairwise distance of each cluster, while searching for pair with least distance
        c1 = 0; c2 = 0; i1 = 0; i2 = 0;
        sdis = 9999999999
        for i in xrange(nc):
            for j in xrange(i+1, nc):
                dis = euclidistance(nodes[i], nodes[j])
                if dis < sdis:
                    sdis = dis
                    c1 = nodes[i]; c2 = nodes[j];
                    i1 = i; i2 = j;
        # Merge these two nodes into one new node
        node = Node(c1.points + c2.points)
        node.left = c1; node.right = c2;
        
        #Remove the previous nodes, and add the new node
        new_nodes = []
        for i in xrange(nc):
            if i != i1 and i != i2:
                new_nodes = new_nodes + [nodes[i]]
        new_nodes = new_nodes + [node]
        nodes = new_nodes[:]
        nc = nc - 1

    return nodes

Kemudian terakhir, kita definisikan fungsi jarak yang kita pakai. Selain jarak euclidean, kita bisa yang digunakan kriteria jarak lainnya seperti Manhattan distance atau Chebyshev distance.

import math

def euclidistance(c1, c2):
    """ Calculate the distance between two cluster """
    dist = .0
    n1 = len(c1.points)
    n2 = len(c2.points)
    for i in xrange(n1):
        for j in xrange(n2):
            p1 = c1.points[i]
            p2 = c2.points[j]
            dim = len(p1)
            d = 0
            for k in xrange(dim):
                d = d + (p1[k]-p2[k])**2
            d = math.sqrt(d)
            dist = dist + d
    dist = dist / (n1*n2)
    return dist

Contoh
Untuk contoh pertama kita pakai dataset sintesis yang dihasilkan oleh dua buah distribusi normal multidimensi.

import upgma
import random
import matplotlib.pyplot as plt
import math

# Example 1
datapoints = [(random.normalvariate(2.5, 1.0), random.normalvariate(1.5,1.0)) for i in xrange(100)] + \
				[(random.normalvariate(-1, 0.5), random.normalvariate(3,0.5)) for i in xrange(100)]

# Plot datapoints before clustering
plt.plot([x for (x,y) in datapoints], [y for (x,y) in datapoints], 'k^')
plt.show()

# Cluster the data
nodes = upgma.upgma(datapoints, 2)
plt.plot([x[0] for x in nodes[0].points], [x[1] for x in nodes[0].points], 'b*')
plt.plot([x[0] for x in nodes[1].points], [x[1] for x in nodes[1].points], 'ro')
plt.show()

dan berikut keluaran hasil clustering jika kita pilih 2 sebagai jumlah cluster.
nocluster1cluster1

Contoh kedua kita ambil dataset old faithful geyser (http://www.stat.cmu.edu/~larry/all-of-statistics/=data/faithful.dat).

f = open('faithful.dat.txt', 'r')
lines = f.readlines()
f.close()

datapoints = []
for line in lines:
	tokens = line.split()
	datapoints += [(float(tokens[1]), float(tokens[2]))]

Karena fitur kedua berbeda skalanya dengan fitur pertama (puluhan berbanding dengan satuan), kita lakukan normalisasi z-score pada tiap masing-masing fitur dengan cara mengurangi dengan mean distribusi dan kemudian dibagi dengan standar deviasi.

avg1 = sum([x for (x,y) in datapoints])
avg2 = sum([y for (x,y) in datapoints])
centered_datapoints = map(lambda (x,y): (x-avg1, y-avg2), datapoints)
std1 = math.sqrt(sum(map(lambda x: x*x, [x for (x,y) in centered_datapoints])))
std2 = math.sqrt(sum(map(lambda x: x*x, [y for (x,y) in centered_datapoints])))
normalized_datapoints = map(lambda (x,y): (x/std1, y/std2), centered_datapoints)

Hasil clustering-nya adalah sebagai berikut.

# Before clustering
plt.plot([x for (x,y) in normalized_datapoints], [y for (x,y) in normalized_datapoints], 'k^')
plt.show()

# Cluster the data
nodes = upgma.upgma(normalized_datapoints, 2)
plt.plot([x[0] for x in nodes[0].points], [x[1] for x in nodes[0].points], 'b*')
plt.plot([x[0] for x in nodes[1].points], [x[1] for x in nodes[1].points], 'ro')
plt.show()

nocluster2cluster2

Kode lengkap dapat dilihat di: https://github.com/mitbal/py-upgma

Semoga berguna. Salam.

Gibbs Sampling

Pada post ini, saya akan menjelaskan mengenai implementasi algoritma Gibbs sampling untuk mendeteksi pola pada deret DNA atau populer dengan istilah motif finding seperti yang dijabarkan oleh Lawrence di paper-nya pada tahun 1993 Detecting Subtle Sequence Signals: A Gibbs Sampling Strategy for Multiple Alignment.

Sekilas mengenai Gibbs sampling. Gibbs sampling adalah salah satu algoritma keluarga dari Markov Chain Monte Carlo (MCMC). Kurang lebih intinya adalah kita bisa menghitung joint probability distribution dengan cara melakukan sampling satu per satu terhadap setiap variabel dengan berdasarkan nilai variabel lainnya alias full conditional probability. Materi pengenalan Gibbs sampling yang lebih lanjut bisa dipelajari di sini dan di sini.

Motif itu sendiri adalah pola yang terdapat pada beberapa deret DNA yang memiliki tingkat kemiripan yang tinggi. Contohnya bisa dilihat pada gambar di bawah. Seperti yang bisa dilihat, terdapat kemunculan pola dengan bentuk CAGGAT pada 3 deret dengan lokasi yang berbeda-beda. Motif finding bertujuan menemukan pola tersebut berdasarkan data semua deret string yang dimiliki. Karakter yang muncul selain pada motif biasa disebut background.

motif

Nah, bagaimana caranya kita bisa menggunakan Gibbs sampling untuk melakukan deteksi motif? Ada 2 asumsi sebelum kita memulai memakan Gibbs sampling. Pertama, pada satu deret DNA hanya ada 1 buah pola dan kedua, kita mengetahui berapa panjang dari pola tersebut. Kita akan implementasi Gibbs sampling dengan menggunakan bahasa pemrograman Python.

Pertama, kita inisialisasi state posisi secara acak

pos = [randint(0, len(seq)-w+1) for seq in sequences]

Kemudian kita lakukan iterasi sampai posisi yang kita cari konvergen. Kita lakukan iterasi terhadap semua deret satu per satu. Sebutlah deret yang aktif pada iterasi ini sebagai i.

    last_pos = None
    while last_pos != pos:
        last_pos = pos[:]
        
        # We pick the sequence, well, in sequence starting from index 0
        for i in xrange(K):

Kemudian dari semua deret yang kita punya, selain deret yang sedang aktif dalam iterasi, kita hitung model probabilitas motif dan background. Model untuk background cukup dihitung dengan menghitung berapa kali karakter suatu karakter muncul. Sedangkan untuk model motif, kita harus menghitung berapa kali sebuah karakter muncul pada posisi ke j pada motif tersebut.

 seq_minus = sequences[:]; del seq_minus[i]
 pos_minus = pos[:]; del pos_minus[i]
 q, p = compute_model(seq_minus, pos_minus, alphabet, w)
...
def compute_model(sequences, pos, alphabet, w):
    """
    This method compute the probability model of background and word based on data in 
    the sequences.
    """
    q = {x:[1]*w for x in alphabet}
    p = {x: 1 for x in alphabet}
    
    # Count the number of character occurrence in the particular position of word
    for i in xrange(len(sequences)):
        start_pos = pos[i]        
        for j in xrange(w):
            c = sequences[i][start_pos+j]
            q[c][j] += 1
    # Normalize the count
    for c in alphabet:
        for j in xrange(w):
            q[c][j] = q[c][j] / float( len(sequences)+len(alphabet) )
    
    # Count the number of character occurrence in background position
    # which mean everywhere except in the word position
    for i in xrange(len(sequences)):
        for j in xrange(len(sequences[i])):
            if j < pos[i] or j > pos[i]+w:
                c = sequences[i][j]
                p[c] += 1
    # Normalize the count
    total = sum(p.values())
    for c in alphabet:
        p[c] = p[c] / float(total)
    
    return q, p

Setelah mendapatkan model untuk motif dan background, kita coba hitung probabilitas untuk semua posisi motif yang mungkin pada deret i dengan cara mengalikannya dengan dua model yang sudah dihitung sebelumnya. Semua karakter dimulai dari posisi awal hingga panjang motif dikalikan dengan nilai probabilitas kemunculan karakter tersebut pada posisi yang bersangkutan.

N = len(sequences[i])
qx = [1]*(N-w+1)
px = [1]*(N-w+1)
for j in xrange(N-w+1):
    for k in xrange(w):
        c = sequences[i][j+k]
        qx[j] = qx[j] * q[c][k]
        px[j] = px[j] * p[c]

Setelah itu, kita bisa menghitung rasio antara motif dan background untuk setiap kemungkinan posisi motif pada deret tersebut. Rasio yang paling besar menandakan bahwa posisi tersebut adalah posisi yang paling mungkin sebagai awal dari motif.

Aj = [x/y for (x,y) in zip(qx, px)]
norm_c = sum(Aj)
Aj = map(lambda x: x/norm_c, Aj)

Kita lakukan sampling posisi baru berdasarkan distribusi rasio. Rasio yang lebih besar tentunya memiliki peluang untuk terpilih lebih besar. Hal ini dilakukan untuk memberikan randomness dan menghindari local maxima (walaupun mungkin kurang berhasil).

pos[i] = sample(range(N-w+1), Aj)
...
def sample(alphabet, dist):
    """ This method produce a new discrete sample list by alphabet with probability
    distribution given by dist.
    The length of alphabet and dist must be same."""
    sampl = None
    cum_dist = np.cumsum(dist)
    r = rand()
    for i in xrange(len(dist)):
        if r < cum_dist[i]:
            sampl = alphabet[i]
            break
    
    return sampl

Hal ini diulang untuk semua deret, dan dilakukan hingga tidak ada lagi perubahan posisi untuk motif.

Kode lengkapnya bisa dilihat di https://github.com/mitbal/gibbs-sampler-motif-finding pada file gibbs.py

Sekarang untuk tes algoritma kita buat satu file python baru dan masukkan kode berikut

seqs = ['muhammadiqbal', 'iqbaltawakal', 'miqbalt']
k = 5

new_pos = gibbs.sampling(seqs, k)

words = [seqs[i][new_pos[i]:new_pos[i]+k] for i in xrange(len(seqs))]
print words

Jika sukses, seharusnya program di atas mengeluarkan output seperti berikut

['iqbal', 'iqbal', 'iqbal']

Namun terkadang Gibbs sampling bisa terjebak di local maxima karena inisialisasi posisi awal yang kurang baik atau pemilihan sampling posisi baru yang kurang beruntung.

['uhamm', 'iqbal', 'iqbal']

Untuk mengurangi risiko tersebut, kita bisa menjalankan gibbs sampling beberapa kali dan menghitung berapa kali output yang dihasilkan keluar. Dengan memilih output yang paling sering keluar kita bisa meningkatkan akurasi dari Gibbs sampling. Istilahnya adalah multiple chains.

result = {}
for i in xrange(20):
    new_pos = gibbs.sampling(seqs, k)
    #print new_pos
    tnp = tuple(new_pos)    
    if tnp in result:
        result[tnp] += 1
    else:
        result[tnp] = 1

max_vote = 0
max_pos = None
for key in result:
    #print key, result[key]
    if result[key] > max_vote:
        max_pos = list(key)
        max_vote = result[key]

words = [seqs[i][max_pos[i]:max_pos[i]+k] for i in xrange(len(seqs))]
print words

Selain cara diatas, module gibbs bisa dipanggil langsung dari command line. Pertama masukkan data mengenai deret dalam sebuah file dengan format baris pertama adalah panjang motif, dan kemudian setiap barisnya adalah satu deret. Contoh file “test.txt”

3
thequickdog
browndog
dogwood

Kemudian panggil program dari terminal dengan cara

 cat test.txt | python gibbs.py

atau

 python gibbs.py test.txt

Output-nya adalah posisi motif untuk setiap deret dalam bentuk list.

Semoga berguna. Salam.