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.

Iklan

[ML] Apa itu Machine Learning

Bermula dari pembicaraan dengan teman-teman saya beberapa waktu yang lalu, motivasi dari post ini adalah untuk memberikan penjelasan mengenai apa itu Machine Learning.

Machine Learning — dialihbahasakan bebas menjadi pembelajaran mesin, atau disingkat menjadi ML — adalah ilmu cabang dari kecerdasan buatan yang mempelajari bagaimana caranya belajar dari data. Istilah Machine Learning sendiri cukup membingungkan atau misleading karena hampir tidak berhubungan dengan mesin apapun (kecuali diimplementasikan di robot). Mesin disini merujuk kepada algoritma atau program yang berjalan di komputer. Istilah lain yang biasa dipakai adalah Data Mining, Pattern Recognition, atau Knowledge Discovery.

Contoh aplikasinya bermacam-macam, sebagai contoh:

1. Bidang kedokteran: bagaimana mendeteksi penyakit seseorang dari gejala-gejala yang ada, atau deteksi apakah seseorang mengidap penyakit jantung dari rekaman elektrokardiogram, dan mencari tahu gen yang terlibat pada penyakit kanker.

2. Bidang computer vision: menemukan dan memberi label muka orang pada foto (seperti di facebook) atau face recognition, pengenalan tulisan tangan menjadi teks pada komputer atau handwriting recognition.

3. Bidang teks atau information retrieval: menerjemahkan bahasa menggunakan mesin atau machine translation, mengubah suara menjadi teks atau speech recognition, atau memisahkan email antara yang spam dan yang non-spam.

Untuk bisa lebih mengenal mengenai maka kita harus tahu terlebih dahulu bahasa atau istilah-istilah yang sering dipakai dalam bidang tersebut.

Untuk bisa mengaplikasikan teknik-teknik machine learning maka hal yang pertama harus ada adalah data. Tanpa data maka algoritma machine learning tidak dapat bekerja. Data yang dimiliki biasanya dibagi menjadi 2, yaitu data training dan data testing. Data training digunakan untuk melatih algoritma, sedangkan data testing dipakai untuk mengetahui performa algoritma yang sudah dilatih sebelumnya ketika menemukan data baru yang belum pernah dilihat sebelumnya. Ini biasanya disebut dengan generalisasi. Hasil dari pelatihan tersebut bisa disebut dengan model.

Dari model tersebut kita bisa melakukan prediksi yang biasanya dibedakan menjadi dua tergantung tipe keluarannya. Jika hasil prediksinya bersifat diskrit, maka proses itu disebut klasifikasi. Contohnya klasifikasi. Jika keluarannya bersifat kontinyu, maka itu disebut dengan regresi. Contohnya prediksi cuaca besok berdasarkan data cuaca 3 bulan terakhir. Performa dari sebuah algoritma machine learning itu sendiri bisa dihitung dari akurasiĀ  atau seberapa banyak prediksi yang dibuat yang sesuai dengan kenyataan.

Skema machine learning secara garis besar
Skema machine learning secara garis besar

Lebih lanjut lagi, algoritma pembelajaran machine learning bisa dibagi menjadi 2, supervised learning dan unsupervised learning. Pada supervised learning, data yang dimiliki dilengkapi dengan label/kelas yang menunjukkan kepada kelompok mana data tersebut seharusnya berada. Misalkan untuk kasus deteksi spam, maka kita harus tahu apakah email ini termasuk spam atau tidak. Pada kasus unsupervised learning, data yang kita miliki tidak memiliki label sehingga yang ingin kita ketahui adalah terdiri dari berapa kelompok kira-kira data yang kita punya berdasarkan kemiripannya. Algoritma yang dipakai disini biasanya adalah algoritma clustering.

Pada post berikutnya saya akan membahas algoritma-algoritma populer yang sering dipakai dalam Machine Learning.

Salam,