🌑 🤷🏾 ♈️ Kami membongkar algoritma-EM menjadi batu bata kecil ☝🏿 💣 🏙️

Pada artikel ini, seperti yang mungkin sudah Anda duga, kita akan berbicara tentang perangkat EM-algoritma. Pertama-tama, artikel ini mungkin menarik bagi mereka yang secara perlahan bergabung dengan komunitas dataainists. Materi yang disajikan dalam artikel akan lebih bermanfaat bagi mereka yang baru saja mulai mengambil kursus ketiga "Mencari struktur data" dalam spesialisasi "Pembelajaran Mesin dan Analisis Data" dari MIPT dan Yandex.

Materi yang disajikan dalam artikel, dalam arti tertentu, merupakan tambahan pada minggu pertama pelatihan dalam kursus tersebut, yaitu, memungkinkan Anda untuk menjawab beberapa pertanyaan penting mengenai prinsip pengoperasian algoritma EM. Untuk pemahaman yang lebih baik tentang materi, disarankan bagi pembaca kami yang terhormat untuk dapat melakukan operasi dengan matriks (perkalian matriks, menemukan penentu matriks dan matriks terbalik), memahami dasar-dasar teori probabilitas dan matstat, dan tentu saja, setidaknya memiliki pemahaman dasar tentang algoritma pengelompokan dasar dan memahami apa pengelompokan terjadi dalam pembelajaran mesin. Meskipun, tentu saja, tanpa pengetahuan ini, Anda dapat membaca artikel, sesuatu dan ya itu akan menjadi jelas :)

Juga, menurut tradisi lama, artikel itu tidak akan berisi penelitian teoritis yang mendalam, tetapi akan diisi dengan contoh-contoh sederhana dan dapat dimengerti. Setiap contoh selanjutnya akan menjelaskan operasi algoritma EM sedikit lebih dalam dari yang sebelumnya, yang pada akhirnya membawa kita langsung ke analisis algoritma itu sendiri. Untuk setiap contoh, kode akan ditulis. Semua kode ditulis dengan python 2.7, dan untuk ini saya minta maaf sebelumnya. Ternyata sekarang saya menggunakan versi khusus ini, tetapi setelah beralih ke python 3, saya akan mencoba mengubah kode dalam artikel.

Mari kita berkenalan dengan garis besar artikel: 1) Mari kita perhatikan secara umum bagaimana algoritma EM bekerja.

2) Kami mengulangi poin utama dari teorema Bayes (rumus):

P_{x_{i}} (B_{j}) = \frac{P (B_{j}) \cdot P_{B_{j}} (x_{i})}{P (x_{i})}

$P_{x_i}(B_j) = \frac{P(B_j) \centerdot P_{B_j}(x_i)}{P(x_i)}$

3) Pertimbangkan contoh pertama pada teorema Bayesian, di mana semua elemen rumus (atau lebih tepatnya probabilitas

P (B_{j}), P_{B_{j}} (x_{i})

$P(B_j), P_{B_j}(x_i)$ ) untuk tanda sama di sebelah kanan diketahui. Dalam hal ini, kita hanya perlu memahami dengan benar angka mana yang harus diganti dan kemudian melakukan operasi dasar.

4) Perhatikan contoh kedua pada teorema Bayes. Untuk mengatasi masalah tersebut, kita perlu juga menghitung kepadatan probabilitas suatu peristiwa tertentu (

x_{i}

$x_i$ ) tunduk pada hipotesis (

B_{j}

$B_j$ ) -

ρ_{B_{j}} (x_{i})

$\rho_{B_j}(x_i)$ . Untuk perhitungan, kita akan diberikan parameter variabel acak, yaitu, ekspektasi matematis -

μ

$\mu$ dan standar deviasi -

σ

$\sigma$ . Perhitungan kepadatan probabilitas akan dilakukan sesuai dengan rumus berikut:

ρ_{B_{j}} (x_{i}) = \frac{1}{σ_{j} \cdot \sqrt{2 π}} \cdot \exp (- \frac{(x_{i} - μ_{j})^{2}}{2 σ_{j}^{2}})

$\rho_{B_j}(x_i) = \frac{1}{\sigma_j \centerdot \sqrt{2 \pi}} \centerdot \exp(- \frac{(x_i-\mu_j)^2}{2 \sigma_j^2})$

Penggunaan rumus di atas mengasumsikan bahwa variabel acak memiliki dimensi satu dimensi dan didistribusikan secara normal (jika tidak mereka mengatakan distribusi Gaussian atau Gaussian-Laplace adalah distribusi probabilitas).

5) Pertimbangkan contoh ketiga, yang berbeda dari yang sebelumnya hanya dalam hal ini mempertimbangkan kasus distribusi normal multidimensi (dalam versi dua dimensi kami) dari variabel acak. Dalam hal ini, perhitungan kepadatan dilakukan sesuai dengan rumus berikut:

ρ_{B_{j}} (x_{i}) = \frac{1}{(2 π)^{k / 2} \cdot ∣ Σ_{j} ∣^{1 / 2}} \cdot \exp (- \frac{1}{2} \cdot (x_{i} - μ_{j})^{T} Σ_{j}^{- 1} (x_{i} - μ_{j}))

$\rho_{B_j}(x_i) = \frac{1}{(2\pi)^{k/2} \centerdot \mid \Sigma_j\mid ^{1/2}} \centerdot \exp(- \frac{1}{2} \centerdot (x_i-\mu_j)^T \Sigma_j^{-1} (x_i-\mu_j))$

6) Akhirnya, kami memodifikasi contoh ketiga sedemikian rupa sehingga dengan jelas menunjukkan operasi algoritma EM.

7) Sebagai kesimpulan, kami membandingkan kualitas algoritme yang kami simpulkan dengan kualitas algoritme EM, yang dibangun di pustaka sklearn (class sklearn.mixture.GaussianMixture)

Jika Anda melihat karakter Jepang alih-alih rumus yang ditentukan dalam garis besar artikel, harap jangan takut - setelah teorema Bayes (paragraf kedua dari garis besar makalah ini), kami akan menganalisis contoh-contoh sedemikian rupa sehingga Anda mungkin akan mengerti setidaknya secara intuitif. Ayo pergi!

Algoritma EM secara umum

Pada kursus, berdasarkan artikel yang ditulis, algoritma EM disajikan sebagai salah satu metode pengelompokan. Dengan kata lain, ini adalah metode pembelajaran mesin tanpa guru, ketika kita tidak tahu jawaban yang benar sebelumnya. Dalam artikel kami, kami juga akan mempertimbangkan algoritma sebagai bagian dari salah satu metode pengelompokan variabel acak dengan distribusi normal. Namun, perlu dicatat bahwa algoritma ini memiliki aplikasi yang lebih luas. Algoritma EM (Bahasa Inggris ekspektasi-maksimisasi) adalah metode umum untuk menemukan perkiraan fungsi likelihood dalam model dengan variabel tersembunyi, yang dari campuran distribusi memungkinkan Anda untuk membangun (perkiraan) distribusi probabilitas kompleks.

Algoritma EM dalam masalah pengelompokan digunakan sebagai algoritma iteratif, yang melakukan dua langkah pada setiap iterasi:

E-step. Pada langkah-E pertama, dalam beberapa cara, misalnya, secara acak, kami memilih variabel tersembunyi, dalam kasus kami itu akan menjadi harapan matematika -

μ

$\mu$ dan standar deviasi -

σ

$\sigma$ . Menggunakan variabel yang dipilih, kami menghitung probabilitas menempatkan setiap objek ke kluster tertentu. Langkah-E berikutnya menggunakan variabel tersembunyi yang didefinisikan dalam langkah-M.

Langkah-M . Pada langkah-M, kita, sesuai dengan probabilitas menempatkan setiap objek ke kluster tertentu yang diperoleh pada langkah-E, menghitung ulang variabel tersembunyi

μ

$\mu$ dan

σ

$\sigma$

Iterasi diulang sampai terjadi konvergensi. Kita dapat mengasumsikan bahwa konvergensi terjadi ketika nilai-nilai variabel tersembunyi tidak berubah secara signifikan, misalnya, dalam konstanta yang diberikan. Dalam contoh terakhir, yang dipertimbangkan dalam artikel, kami tidak akan menetapkan konstanta dan, oleh karena itu, kami tidak akan menentukan perubahan dalam nilai variabel tersembunyi di setiap iterasi. Mari kita lakukan dengan lebih sederhana - kami akan membatasi algoritme kami hingga sejumlah langkah tetap.

Langkah-M cukup sederhana untuk dipahami, dan kita akan melihat ini dalam contoh terakhir. Kami akan fokus pada E-step dalam artikel.

Langkah-E didasarkan pada teorema Bayesian, atas dasar yang kita tentukan probabilitas setiap objek untuk ditugaskan ke satu atau lain cluster. Mari kita ingat tentang teorema ini.

Ingat formula Bayes

P_{x_{i}} (B_{j}) = \frac{P (B_{j}) \cdot P_{B_{j}} (x_{i})}{P (x_{i})}

$P_{x_i}(B_j) = \frac{P(B_j) \centerdot P_{B_j}(x_i)}{P(x_i)}$

dimana

P_{x_{i}} (B_{j})

$P_{x_i}(B_j)$ - probabilitas bahwa dengan acara tersebut

x_{i}

$x_i$ ada hipotesis

B_{j}

$B_j$

P (B_{j})

$P(B_j)$ - probabilitas hipotesis

B_{j}

$B_j$

P_{B_{j}} (x_{i})

$P_{B_j}(x_i)$ - probabilitas terjadinya suatu peristiwa

x_{i}

$x_i$ tunduk pada hipotesis

B_{j}

$B_j$

P (x_{i})

$P(x_i)$ - probabilitas terjadinya suatu peristiwa

x_{i}

$x_i$ tunduk pada pemenuhan setiap hipotesis (dihitung dengan rumus probabilitas total acara)

Asalkan acara tersebut

x_{i}

$x_i$ sudah terjadi, kemungkinan memenuhi hipotesis

B_{j}

$B_j$ dinilai kembali menggunakan rumus di atas. Kita bisa mengatakan itu

P (B_{j})

$P(B_j)$ ini adalah probabilitas a priori dari hipotesis (sebelum ujian), dan

P_{x_{i}} (B_{j})

$P_{x_i}(B_j)$ Apakah probabilitas posterior dari hipotesis yang sama dievaluasi setelah kejadian

x_{i}

$x_i$ yaitu, mengingat fakta bahwa acara tersebut

x_{i}

$x_i$ andal terjadi.

Teorema Contoh Satu Oleh Bayes

Bayangkan kami menerima bagian yang diproduksi pada dua mesin yang berbeda. Karakteristik berikut dari produk yang diterima di gudang diketahui:

Secara total, bagian - 10.000 pcs (

N = 10000

$N = 10000$ ), bagian mana yang diproduksi pada mesin pertama - 6000 pcs. (

N 1 = 6000

$N1 = 6000$ ), pada detik - 4000 pcs. (

N 2 = 4000

$N2 = 4000$ )

Proporsi produk standar (mis., Tidak cacat) yang diproduksi pada mesin pertama adalah 0,9, pangsa produk standar yang diproduksi pada mesin kedua adalah 0,8.

Kami secara acak mengekstraksi satu bagian dari bagian yang diterima, dan ternyata menjadi standar.

Kita perlu menentukan probabilitas bahwa:

1) bagian diproduksi pada mesin 1;

2) barang itu diproduksi di mesin ke-2.

Keputusan

1) Acara

x

$x$ dalam contoh kita, mengekstraksi produk standar.

2) Sekarang, mari kita tentukan hipotesisnya

B_{j}

$B_j$ . Kami memiliki dua bagian, yang berarti dua hipotesis:

Hipotesis

B_{1}

$B_1$ : produk yang tidak sengaja dilepas diproduksi di mesin No. 1.

Hipotesa

B_{2}

$B_2$ : produk yang tidak sengaja dilepas diproduksi pada mesin No. 2.

Dengan demikian, kemungkinan mengekstraksi produk yang diproduksi pada mesin pertama adalah

P (B_{1})

$P(B_1)$ membuat

N 1 / N = 6000 / 10000 = 0.6

$N1/N = 6000/10000 = 0.6$ .

Probabilitas mengekstraksi produk yang diproduksi pada mesin kedua adalah

P (B_{2})

$P(B_2)$ membuat

N 2 / N = 0.4

$N2/N = 0.4$ . Dan kami tidak peduli dengan produk standar atau yang cacat, penting dari batch mana itu.

3) Probabilitas mengekstraksi produk standar dari produk yang diproduksi pada mesin No. 1 (yaitu, berdasarkan hipotesis pertama) sesuai dengan proporsi produk standar yang diproduksi pada mesin No. 1 dan 0,9,

P_{B_{1}} (x_{i}) = 0.9

$P_{B_1}(x_i) = 0.9$ .

Probabilitas mengekstraksi produk adalah subjek produk standar untuk hipotesis No. 2 -

P_{B_{2}} (x_{i}) = 0.8

$P_{B_2}(x_i) = 0.8$

4) Tentukan kemungkinan mengekstraksi produk standar dari seluruh rangkaian produk -

P (x_{i})

$P(x_i)$ . Sesuai dengan rumus untuk probabilitas total suatu peristiwa

P (x_{i}) = P (B_{1}) \cdot P_{B_{1}} (x_{i}) + P (B_{2}) \cdot P_{B_{2}} (x_{i}) = 0.6 \cdot 0.9 + 0.4 \cdot 0.8 = 0.54 + 0.32 = 0.86

$P(x_i) = P(B_1) \centerdot P_{B_1}(x_i) + P(B_2) \centerdot P_{B_2}(x_i) = 0.6 \centerdot 0.9 + 0.4 \centerdot 0.8 = 0.54 + 0.32 = 0.86$ . Di sini kita mengerti itu

0.14

$0.14$ - ini adalah kemungkinan mengekstraksi bagian yang rusak dari semua bagian yang tiba di gudang dan total ternyata 1

(0.86 + 0.14 = 1)

$(0.86 + 0.14 = 1)$ .

5) Sekarang kita memiliki semua data untuk menyelesaikan masalah.

5.1) Menentukan probabilitas bahwa bagian standar yang diekstraksi secara acak dihasilkan pada mesin

№ 1

$№1$ :

P_{x_{i}} (B_{1}) = \frac{P (B_{1}) \cdot P_{B_{1}} (x_{i})}{P (x_{i})} = \frac{0.6 \cdot 0.9}{0.86} \approx 0.63

$P_{x_i}(B_1) = \frac{P(B_1) \centerdot P_{B_1}(x_i)}{P(x_i)} = \frac{0.6 \centerdot 0.9}{0.86} \approx 0.63$

5.2) Kami menentukan probabilitas bahwa bagian standar yang diekstraksi secara acak diproduksi pada mesin

№ 2

$№2$ :

P_{x_{i}} (B_{2}) = \frac{P (B_{2}) \cdot P_{B_{2}} (x_{i})}{P (x_{i})} = \frac{0.4 \cdot 0.8}{0.86} \approx 0.37

$P_{x_i}(B_2) = \frac{P(B_2) \centerdot P_{B_2}(x_i)}{P(x_i)} = \frac{0.4 \centerdot 0.8}{0.86} \approx 0.37$

Jadi, kami mengevaluasi kembali hipotesis

B_{1}

$B_1$ dan

B_{2}

$B_2$ . Setelah penilaian ulang, hipotesis juga membentuk kelompok peristiwa lengkap:

0.63 + 0.37 = 1

$0.63 + 0.37 = 1$ .

Nah, sekarang adalah saatnya untuk beralih ke contoh kedua.

Contoh kedua untuk teorema Bayes menggunakan parameter distribusi normal dari variabel acak: ekspektasi matematis $\mu$ dan standar deviasi $\sigma$

Mari kita sedikit memodifikasi kondisi dari contoh sebelumnya. Kami berasumsi bahwa kami tidak mengetahui proporsi produk standar yang diproduksi pada mesin No. 1 dan No. 2, tetapi anggaplah, sebaliknya, kami tahu dimensi rata-rata diameter produk dan standar deviasi diameter produk pada setiap mesin.

Mesin No. 1 menghasilkan bagian dengan diameter 64 mm dan standar deviasi 4 mm.

Mesin No. 2 menghasilkan bagian dengan diameter 52 mm dan standar deviasi 2 mm.

Suatu kondisi penting adalah bahwa seluruh rangkaian produk dijelaskan oleh distribusi normal atau distribusi Gaussian.

Kalau tidak, kondisinya sama, kami menuliskannya:

N 1 = 6000, μ_{1} = 64, σ_{1} = 4

$N1 = 6000, \mu_1 = 64, \sigma_1 = 4$

N 2 = 4000, μ_{2} = 52, σ_{2} = 2

$N2 = 4000, \mu_2 = 52, \sigma_2 = 2$

Kami berasumsi bahwa dalam proses penerimaan produk ada insiden kecil, akibatnya semua produk dicampur.

Tugas kita adalah memilah-milah detail dan untuk masing-masing menentukan probabilitas bahwa itu diproduksi pada mesin No. 1 atau mesin No. 2. Kami juga akan mengasumsikan bahwa bagian itu diproduksi pada mesin, kemungkinan yang lebih tinggi.

Solusi

Pertama-tama, kami akan menganalisis algoritma solusi. Kita dapat dengan mudah menghitung probabilitas setiap hipotesis, namun, kita sudah tahu nilainya dari contoh sebelumnya:

P (B_{1}) = 0.6

$P(B_1)=0.6$ ,

P (B_{2}) = 0.4

$P(B_2)=0.4$ . Tetapi bagaimana kita menghitung probabilitas merebut produk standar secara individual untuk masing-masing hipotesis yang tersedia? Semuanya sederhana! Kami tidak akan mempertimbangkan probabilitas, sebagai gantinya, kami akan menentukan kepadatan probabilitas item yang akan dihapus dengan nilai diameternya untuk setiap hipotesis secara terpisah.

Untuk melakukan ini, kami menggunakan rumus terkenal:

ρ_{B_{j}} (x_{i}) = \frac{1}{σ_{j} \cdot \sqrt{2 π}} \cdot \exp (- \frac{(x_{i} - μ_{j})^{2}}{2 σ_{j}^{2}})

$\rho_{B_j}(x_i) = \frac{1}{\sigma_j \centerdot \sqrt{2 \pi}} \centerdot \exp(- \frac{(x_i-\mu_j)^2}{2 \sigma_j^2})$

dimana

x_{i}

$x_i$ - nilai acak (dalam kasus kami, ukuran sebenarnya diameter produk),

μ_{j}

$\mu_j$ - Ekspektasi matematis variabel acak

j

$j$ hipotesis (dalam kasus kami, ukuran rata-rata diameter produk yang dihasilkan

j

$j$ alat mesin)

σ_{j}

$\sigma_j$ - standar deviasi variabel acak

j

$j$ hipotesis (dalam kasus kami, penyimpangan dari ukuran rata-rata diameter produk yang dihasilkan

j

$j$ alat mesin).

Dengan demikian, kami membuat substitusi pintar probabilitas untuk kepadatan probabilitas dalam pembilang rumus Bayes, kami akan melakukan penggantian serupa dalam penyebut.

Kami menyesuaikan formula khusus untuk kasus kami.

Kemungkinan bagian dengan dimensi diameter

x_{i}

$x_i$ diproduksi pada mesin No. 1 ditentukan oleh rumus:

P_{x_{i}} (B_{1}) = \frac{w_{1} \cdot ρ_{B_{1}} (x_{i})}{w_{1} \cdot ρ_{B_{1}} (x_{i}) + w_{2} \cdot ρ_{B_{1}} (x_{i})}

$P_{x_i}(B_1) = \frac{w_1 \centerdot \rho_{B_1}(x_i)}{w_1 \centerdot \rho_{B_1}(x_i) + w_2 \centerdot \rho_{B_1}(x_i)}$

Kemungkinan bagian dengan dimensi diameter

x_{i}

$x_i$ diproduksi pada mesin No. 2:

P_{x_{i}} (B_{2}) = \frac{w_{2} \cdot ρ_{B_{2}} (x_{i})}{w_{1} \cdot ρ_{B_{1}} (x_{i}) + w_{2} \cdot ρ_{B_{1}} (x_{i})}

$P_{x_i}(B_2) = \frac{w_2 \centerdot \rho_{B_2}(x_i)}{w_1 \centerdot \rho_{B_1}(x_i) + w_2 \centerdot \rho_{B_1}(x_i)}$

Perhatikan bahwa kami telah mengganti notasi untuk probabilitas hipotesis dengan

P (B_{j})

$P(B_j)$ di

w_{j}

$w_j$ . Ini karena penunjukan yang sesuai pada kursus yang ditunjuk sebelumnya. Selanjutnya kita akan menggunakan notasi seperti itu saja.

Sekarang kami 100% selesai dan siap untuk menyelesaikan masalah.

Kami melihat kodenya

Impor perpustakaan

#  , 
import numpy as np
import random
import math
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.metrics import accuracy_score

Fungsi yang Digunakan dalam Contoh Kedua

#      ,    
#      : .,  . 
def gaus_func_01(mu,sigma,x):
    return math.e**(-(x-mu)**2/(2*sigma**2)) / (sigma*(2*math.pi)**0.5)

#        
def proba_x(X, w1, w2, mu_1, mu_2, sigma_1, sigma_2):
    for i in X:
        P1_x = gaus_func_01(mu_1,sigma_1,i)
        P2_x = gaus_func_01(mu_2,sigma_2,i)
        P_x = w1*P1_x + w2*P2_x
        P_x_1 = (w1*P1_x)/P_x
        P_x_2 = (w2*P2_x)/P_x
        proba_temp = []
        proba_temp.append(P_x_1)
        proba_temp.append(P_x_2)
        proba_X.append(proba_temp)
    return proba_X

#        
def pred_x(proba_X, limit_proba):
    pred_X = []
    for x in proba_X:
        if x[0] >= limit_proba:
            pred_X.append(1)
        else:
            pred_X.append(2)
    return np.array(pred_X)

#    
def graph_01(X, pred_X, mu_1, sigma_1, mu_2, sigma_2):
    true_pred = []
    false_pred_1 = []
    false_pred_2 = []
    for i in range(X.shape[0]):
        if pred_X[i] == y[i]:
            true_pred.append([X[i], -0.025])
        else:
            if y[i] == 1:
                false_pred_1.append([X[i], -0.0075])
            else:
                false_pred_2.append([X[i], -0.015])

    false_pred_1 = np.array(false_pred_1)            
    false_pred_2 = np.array(false_pred_2)
    true_pred = np.array(true_pred)

    x_theory = np.linspace(42, 85, 20000)
    y_theory_1 = []
    for x in x_theory:
        y_theory_1.append(gaus_func_01(mu_1,sigma_1,x))
    y_theory_2 = []
    for x in x_theory:
        y_theory_2.append(gaus_func_01(mu_2,sigma_2,x))

    plt.figure(figsize=(18, 8))    
    plt.plot(
        x_theory, y_theory_1, color = 'green', lw = 2, label = 'Theoretical probability density for machine 1')
    plt.plot(
        x_theory, y_theory_2, color = 'firebrick', lw = 2, label = 'Theoretical probability density for machine 2')
    plt.hist(
        X[:N1], bins = 'auto', color='#539caf', normed = True, alpha = 0.35, label = 'machine tool products 1')
    plt.hist(
        X[N1:N], bins = 'auto', color='sandybrown', normed = True, alpha = 0.75, label = 'machine tool products 2')
    plt.plot(mu_1, 0, 'o', markersize = 11, color = 'blue', label = 'Mu 1')
    plt.plot(mu_2, 0, 'o', markersize = 11, color = 'red', label = 'Mu 2')

    plt.plot([mu_1 - sigma_1, mu_1 - sigma_1], [0,0.85*np.max(y_theory_1)],
             ':', lw = 3, color = 'blue', alpha = 0.55, label = 'Mu1 - sigma1')
    plt.plot([mu_1 + sigma_1, mu_1 + sigma_1], [0,0.85*np.max(y_theory_1)],
             ':', lw = 3, color = 'blue', alpha = 0.55, label = 'Mu1 + sigma1')
    plt.plot([mu_2 - sigma_2, mu_2 - sigma_2], [0,0.85*np.max(y_theory_2)],
             ':', lw = 3, color = 'red', alpha = 0.55, label = 'Mu2 - sigma2')
    plt.plot([mu_2 + sigma_2, mu_2 + sigma_2], [0,0.85*np.max(y_theory_2)],
             ':', lw = 3, color = 'red', alpha = 0.55, label = 'Mu2 + sigma2')

    plt.plot([mu_1 - 2 * sigma_1, mu_1 - 2 * sigma_1], [0, 0.9*0.5 * np.max(y_theory_1)],
             ':', lw = 2.5, color = 'blue', alpha = 0.35, label = 'Mu1 - 2*sigma1')
    plt.plot([mu_1 + 2 * sigma_1, mu_1 + 2 * sigma_1], [0, 0.9*0.5 * np.max(y_theory_1)],
             ':', lw = 2.5, color = 'blue', alpha = 0.35, label = 'Mu1 + 2*sigma1')
    plt.plot([mu_2 - 2 * sigma_2, mu_2 - 2 * sigma_2], [0, 0.9*0.5 * np.max(y_theory_2)],
             ':', lw = 2.5, color = 'red', alpha = 0.35, label = 'Mu2 - 2*sigma2')
    plt.plot([mu_2 + 2 * sigma_2, mu_2 + 2 * sigma_2], [0, 0.9*0.5 * np.max(y_theory_2)],
             ':', lw = 2.5, color = 'red', alpha = 0.35, label = 'Mu2 + 2*sigma2')

    plt.plot(false_pred_1[:,0], false_pred_1[:,1], 'o', markersize = 2.5, color = 'blue', alpha = 0.2, label = 'errors1')
    plt.plot(false_pred_2[:,0], false_pred_2[:,1], 'o', markersize = 2.5, color = 'red', alpha = 0.3, label = 'errors2')
    plt.plot(true_pred[:,0], true_pred[:,1], 'o', markersize = 3, color = 'green', alpha = 0.2, label = 'right answers')

    plt.xlabel('Caliber')
    plt.ylabel('Probability density')
    plt.legend()
    plt.show()

Tentukan parameternya

#    
#      №1
N1 = 6000
#      №2
N2 = 4000
#      
N = N1+N2

#    №1
mu_1 = 64.
#        №1
sigma_1 = 3.5

#    №2
mu_2 = 52
#        №2
sigma_2 = 2.

Kami akan menghitung

X = np.zeros((N))
np.random.seed(seed=42)
#    ,   №1
X[:N1] = np.random.normal(loc=mu_1, scale=sigma_1, size=N1)
#  ,   №2
X[N1:N] = np.random.normal(loc=mu_2, scale=sigma_2, size=N2)

#   
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#     ,    №1
w1 = float(N1)/N
#     ,    №2
w2 = float(N2)/N

#           
proba_X = []
proba_X = proba_x(X, w1, w2, mu_1, mu_2, sigma_1, sigma_2)

#   ,   ,        
limit_proba = 0.5

#     
pred_X = []
pred_X = pred_x(proba_X, limit_proba)

#    
print '   :', round(accuracy_score(y, pred_X),3)
print

print ' №1'
graph_01(X, pred_X, mu_1, sigma_1, mu_2, sigma_2)

Bagan # 1

Apa yang baru saja kita lakukan? Kami menghasilkan dengan cara pseudo-acak 10.000 nilai yang dijelaskan oleh distribusi normal, yang 6.000 nilai dengan rata-rata 64 mm dan standar deviasi 4 mm, nilai 4.000 dengan rata-rata 52 mm dan standar deviasi 2 mm. Selanjutnya, sesuai dengan algoritma di atas, untuk setiap bagian, probabilitas produksinya pada mesin No. 1 atau No. 2 ditentukan. Setelah itu, kami memilih hipotesis tentang produksi produk pada mesin tertentu, tergantung pada hipotesis mana yang paling mungkin. Akhirnya, kami membandingkan hasil algoritma kami dengan jawaban yang benar.

Di bagian keluaran, kami mendapat bagian jawaban yang benar - 0,986. Secara umum, ini sangat bagus, mengingat kami melakukan pelatihan tanpa menggunakan jawaban yang benar.

Mari kita lihat tabel. Apa yang direkomendasikan untuk diperhatikan? Lihat di mana produk yang tidak diidentifikasi dengan benar oleh algoritma berada.

Pertama, kami melihat kesalahan algoritma hanya di area di mana produk yang diproduksi pada mesin yang berbeda memiliki diameter yang sama. Terlihat cukup logis.

Kedua, kita melihat bahwa algoritma hanya keliru pada objek-objek yang paling jauh dari nilai rata-rata sebenarnya dari objek dan pada saat yang sama cukup dekat dengan pusat salah.

Tetapi hal yang paling menarik adalah bahwa masalah dimulai terutama setelah melewati nilai yang kira-kira sama

μ + - 2 σ

$\mu +- 2\sigma$ , khususnya dalam kasus kami di persimpangan

μ_{1} - 2 σ_{1}

$\mu_1 - 2\sigma_1$ dan

μ_{2} + 2 σ_{2}

$\mu_2 + 2\sigma_2$

Mereka yang ingin dapat "memutar" parameter dan melihat bagaimana kualitas algoritma berubah. Ini akan berguna untuk asimilasi materi yang lebih baik.

Baiklah, kita beralih ke contoh berikut

Contoh ketiga. Kasus dua dimensi

Kali ini kami akan mempertimbangkan kasus ketika kami memiliki data tidak hanya pada diameter produk, tetapi juga, misalnya, pada berat masing-masing produk. Biarkan mesin No. 1 menghasilkan bagian dengan berat 12 g dan standar deviasi 1 g, mesin No. 2 menghasilkan produk dengan berat 10 g dan standar deviasi 0,8 mm.

Untuk sisanya, kita akan menggunakan data dari contoh sebelumnya, kita akan menuliskannya.

N 1 = 6000, μ_{11} = 64, σ_{11} = 4, μ_{12} = 14, σ_{12} = 1

$N1 = 6000, \mu_{1 1} = 64, \sigma_{1 1} = 4, \mu_{1 2} = 14, \sigma_{1 2} = 1$

N 2 = 4000, μ_{21} = 52, σ_{21} = 2, μ_{22} = 10, σ_{22} = 0.9

$N2 = 4000, \mu_{2 1} = 52, \sigma_{2 1} = 2, \mu_{2 2} = 10, \sigma_{2 2} = 0.9$

Situasinya sama - semua detail dicampur dalam satu tumpukan besar dan kita perlu memilahnya dan menentukan untuk setiap bagian kemungkinan produksinya pada mesin No. 1 dan mesin No. 2.

Solusi

Secara umum, pendekatan solusi untuk contoh sebelumnya tidak berubah, dengan pengecualian rumus untuk menentukan kepadatan probabilitas variabel acak. Untuk kasus multidimensi, rumus berikut digunakan:

p_{B_{j}} (x_{i}) = \frac{1}{(2 π)^{k / 2} \cdot ∣ Σ_{j} ∣^{1 / 2}} \cdot \exp (- \frac{1}{2} \cdot (x_{i} - μ_{j})^{T} Σ_{j}^{- 1} (x_{i} - μ_{j}))

$p_{B_j}(x_i) = \frac{1}{(2\pi)^{k/2} \centerdot \mid \Sigma_j\mid ^{1/2}} \centerdot \exp(- \frac{1}{2} \centerdot (x_i-\mu_j)^T \Sigma_j^{-1} (x_i-\mu_j))$

dimana

Σ_{j}

$\Sigma_j$ Merupakan matriks kovarian dari variabel acak

j

$j$ hipotesis (dalam kasus kami, matriks kovarian dari parameter produk yang dihasilkan

j

$j$ alat mesin)

∣ Σ_{j} ∣

$\mid \Sigma_j\mid$ Merupakan penentu matriks kovarians variabel acak

j

$j$ hipotesa

μ_{j}

$\mu_j$ - Ekspektasi matematis variabel acak

j

$j$ hipotesis (dalam kasus kami, nilai rata-rata produk yang diproduksi

j

$j$ alat mesin)

x_{i}

$x_i$ - variabel acak (dalam kasus kami, parameternya

i

$i$ -th product)

Kali ini, untuk memahami kode, pembaca yang terhormat akan membutuhkan pengetahuan tentang operasi dengan matriks

Fungsi yang Digunakan dalam Contoh Ketiga

#      ,    
#      : .,  . 
def gaus_func_02(k, m, x, mu, sigma):
    pj_xi = []
    for j in range(k):
        det_sigma_j = np.linalg.det(sigma[j])
        factor_1 = 1 / (((2 * math.pi)**(k/2)) * ((det_sigma_j)**0.5))
        factor_2 = []
        for i in x:
#    ,            
#          
#   .    ,      
#     
            factor_2.append(math.e**float(
            -0.5 * np.matrix(i - mu[j]) * np.matrix(np.linalg.inv(sigma[j])) * np.matrix(i - mu[j]).T))
        pj_xi.append(factor_1 * np.array(factor_2))
    return np.array(pj_xi)

#     ,       №1  №2
def proba_func_02(pjxi, w, k):
    #          
    P_X = []
    for j in range(k):
        P_X.append(w[j] * pjxi[j])
    P_X = np.sum(np.array(P_X), axis = 0)
    #    ,       №1  №2
    P_J_X = []
    for j in range(k):
        P_J_X.append(w[j] * pjxi[j] / P_X)
    return np.array(P_J_X)

#        
def pred_x_02(proba_X, limit_proba):
    pred_X = []
    for x in proba_X[0]:
        if x >= limit_proba:
            pred_X.append(1)
        else:
            pred_X.append(2)
    return np.array(pred_X)

#             
def graph_02_algorithm(pred_X, mu):
    #   
    pred_X = np.array(pred_X)

    #   ,         
    answers_1 = []
    answers_2 = []

    for i in range(pred_X.shape[0]):
        if pred_X[i] == 1:
            answers_1.append(X[i])
        else:
            answers_2.append(X[i])
    
    print ' "     "'
    plt.figure(figsize=(16, 6))  
    plt.plot(
        np.array(answers_1)[:,0], np.array(answers_1)[:,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        np.array(answers_2)[:,0], np.array(answers_2)[:,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'slateblue', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()
    
#          
def graph_02_true(X, mu):
    print ' "  "'
    plt.figure(figsize=(16, 6))  
    plt.plot(
        X[0:N1,0], X[0:N1,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        X[N1:N,0], X[N1:N,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'slateblue', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()

Tentukan parameternya

#  
k = 2

#      №1
N1 = 6000
#      №2
N2 = 4000
N = N1+N2

#   (  )
m = 2

#    №1
mu_1_1 = 64.
#    №1
mu_1_2 = 14.
#  
sigma_1_1 = 3.5
sigma_1_2 = 1.

#    №2
mu_2_1 = 52.
#    №2
mu_2_2 = 9.5
#  
sigma_2_1 = 2.
sigma_2_2 = 0.7

Kami akan menghitung

X = np.zeros((N, m))
np.random.seed(seed=42)
#    ,   №1
X[:N1, 0] = np.random.normal(loc=mu_1_1, scale=sigma_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_1_2, scale=sigma_1_2, size=N1)
#  ,   №2
X[N1:N, 0] = np.random.normal(loc=mu_2_1, scale=sigma_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_2_2, scale=sigma_2_2, size=N2)

#    (   ,    )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#           (  )
mu  = np.array(([mu_1_1, mu_1_2], [mu_2_1, mu_2_2]))

#        (  )
sigma = np.array(([sigma_1_1, 0.],[0., sigma_1_2], [sigma_2_1, 0.],[0., sigma_2_2]))
sigma = sigma.reshape(k, m, m)

#     ,    №1  №2
w = np.array([float(1./k), float(1./k)])

#   
pj_xi = gaus_func_02(k, m, X, mu, sigma)
proba_X = proba_func_02(pj_xi, w, m)

#   ,   ,        
limit_proba = 0.5

pred_X = pred_x_02(proba_X, limit_proba)
        
#    
print '   :', round(accuracy_score(y, pred_X),3)
print

graph_02_algorithm(pred_X, mu)

graph_02_true(X, mu)

Grafik No. 2.1 "Distribusi produk sesuai dengan algoritma"

Grafik No. 2.2 "Distribusi produk yang sebenarnya"

Dengan analogi dengan contoh sebelumnya, kami menghasilkan 10.000 nilai sesuai dengan parameter di atas

μ

$\mu$ dan

σ

$\sigma$ , menulis beberapa fungsi untuk pengoperasian algoritme kami dan meluncurkannya. Perbedaan mendasar antara kode ini dan kode dari contoh sebelumnya adalah bahwa kali ini kami menggunakan ekspresi matriks untuk perhitungan.

Hasilnya, algoritma kami menunjukkan proporsi jawaban yang benar sama dengan 0,998, yang sebenarnya cukup baik.

Masalahnya, seperti yang kita lihat, semuanya sama - kesalahan dalam detail yang dihasilkan pada mesin yang berbeda dan pada saat yang sama memiliki dimensi dan bobot yang serupa.

Kode ini ditulis sehingga pembaca dapat mengganti nilai parameter apa pun

μ

$\mu$ dan

σ

$\sigma$ dan lihat bagaimana algoritme akan bekerja: dalam hal ini kualitasnya akan menurun, di mana ia akan meningkat. Hal utama adalah jangan lupa mempelajari jadwal. Yah, kita sedang bergerak ke titik akhir kita dalam mem-parsing algoritma EM, sebenarnya algoritma EM itu sendiri.

Memenuhi algoritma EM

Kami melanjutkan contoh kami dengan bagian-bagian yang tiba di gudang. Tetapi kali ini kita hanya akan tahu bahwa produk-produk tersebut diproduksi pada dua mesin yang berbeda, ada 10.000 di antaranya, masing-masing bagian memiliki diameter dan ukuran, dan kita tidak tahu apa-apa lagi. Tetapi tugasnya tidak berubah - kita, seperti sebelumnya, dari tumpukan besar, barang-barang campuran secara acak, kita perlu menentukan mesin mana yang termasuk bagian ini atau itu.

Pada pandangan pertama, ini terdengar hampir tidak realistis, tetapi pada kenyataannya kita memiliki di tangan kita toolkit yang kuat: rumus Bayes dan rumus kepadatan probabilitas dari variabel acak. Mari kita ambil semua kebaikan ini.

Solusi

Apa yang akan kita lakukan? Seperti seharusnya dalam algoritma-EM, kita pertama-tama menginisialisasi parameter:

Probabilitas hipotesis adalah untuk mengekstraksi bagian yang dihasilkan pada mesin No. 1 -

w_{1}

$w_1$ kami menentukan kemungkinan hipotesis untuk mengekstraksi bagian yang diproduksi pada mesin No. 2 -

w_{2}

$w_2$ . Hanya ada dua hipotesis, yang berarti masing-masing pada langkah pertama akan sama dengan 0,5.

Ekspektasi matematis dari variabel acak

μ

$\mu$ tentukan sebagai berikut. Kami mencampur semua produk menggunakan fungsi acak, membagi set sama menjadi dua bagian, untuk setiap bagian untuk setiap parameter (diameter, berat) kami menentukan nilai rata-rata.

Simpangan baku

σ

$\sigma$ mengambil apa yang disebut dari langit-langit - mengaturnya sama dengan persatuan dalam semua hal. Kami menulis dalam format matriks kovarians.

Kami siap membuat E-langkah pertama dari algoritma. Menggunakan parameter yang diinisialisasi dari variabel acak, kami menentukan probabilitas setiap bagian yang ditugaskan ke mesin No. 1 atau mesin No. 2.

Sebenarnya, dengan cara ini, kami mengambil E-langkah pertama.

Sekarang terserah langkah-M. Semuanya sederhana di sini. Setelah kami menentukan probabilitas setiap bagian yang diproduksi pada mesin tertentu, kami dapat menghitung ulang probabilitas setiap hipotesis -

w_{1}

$w_1$ ,

w_{2}

$w_2$ , dan

μ

$\mu$ dan

σ

$\sigma$ .

Kami akan melakukan 15 iterasi seperti itu, dalam dua langkah masing-masing.

Kami melihat kode.

Fungsi yang Digunakan dalam Algoritma EM

#   E-
def e_step(x, k, m, n, w, mu, sigma):
    #      i-     j-  
    pj_xi = []
    for j in range(k):
        det_sigma_j = np.linalg.det(sigma[j])
        factor_1 = 1 / (((2 * math.pi)**(k/2)) * ((det_sigma_j)**0.5))
        factor_2 = []
        for i in x:
            factor_2.append(math.e**float(
                -0.5 * np.matrix(i - mu[j]) * np.matrix(np.linalg.inv(sigma[j])) * np.matrix(i - mu[j]).T))
        pj_xi.append(factor_1 * np.array(factor_2))
    pj_xi = np.array(pj_xi)
    
    #     ,  i-    j- 
    pj_xi_w = []
    for j in range(k):
        pj_xi_w.append(pj_xi[j] * w[j])
    pj_xi_w = np.array(pj_xi_w)
    
    #     i-      
    sum_pj_xi_w = np.sum(pj_xi_w, axis = 0)
    
    #         
    proba_xi = []
    for j in range(k):
        proba_xi.append(pj_xi_w[j] / sum_pj_xi_w)
    
    return np.array(proba_xi)

#  ,    ,            ,
#       
def x_new(proba_xi):
    X1_new_ind = []
    X2_new_ind = []
    X_answers = []

    count = 0
    for x in proba_xi[0]:
        if x >= 0.5:
            X1_new_ind.append(count)
            X_answers.append(1)
        else:
            X2_new_ind.append(count)
            X_answers.append(2)
        count += 1
    #      ,    №1, №2   
    return X1_new_ind, X2_new_ind, X_answers


#   M-
def m_step(x, proba_xi,N):
    w_new = np.sum(proba_xi, axis = 1) / N
    
    #   
    mu_new = (np.array((np.matrix(proba_xi) * np.matrix(X))).T / np.sum(proba_xi, axis = 1)).T
    
    #  
    cov_new = []
    for mu in range(mu_new.shape[0]):
        X_cd = []
        X_cd_proba = []
        count = 0
        for x_i in x:
            cd = np.array(x_i - mu_new[mu])
            X_cd.append(cd)
            X_cd_proba.append(cd * proba_xi[mu][count])
            count += 1
        X_cd = np.array(X_cd)
        X_cd = X_cd.reshape(N, m)
        X_cd_proba = np.array(X_cd_proba)
        X_cd_proba = X_cd_proba.reshape(N, m)

        cov_new.append(np.matrix(X_cd.T) * np.matrix(X_cd_proba))
    cov_new = np.array((np.array(cov_new) / (np.sum(proba_xi, axis = 1)-1)))
    #             ,    
    #       ,        
    if cov_new[0][0][1] < 0:
        cov_new[0][0][1] = 0
    if cov_new[0][1][0] < 0:
        cov_new[0][1][0] = 0
    
    if cov_new[1][0][1] < 0:
        cov_new[1][0][1] = 0
    if cov_new[1][1][0] < 0:
        cov_new[1][1][0] = 0
    
    #   
    sigma_new = cov_new**0.5
    return w_new, mu_new, sigma_new

Tentukan parameternya

#   
#      №1 ( 1)
N1 = 6000
#      №2 ( 2)
N2 = 4000
#       
N = N1 + N2

#  
k = 2

#    №1
mu_samples_1_1 = 64.
#    №1
mu_samples_1_2 = 14.

#    №2
mu_samples_2_1 = 52.
#    №2
mu_samples_2_2 = 9.5

#      №1
sigma_samples_1_1 = 3.5
#      №1
sigma_samples_1_2 = 1.

#      №2
sigma_samples_2_1 = 2.
#      №2
sigma_samples_2_2 = 0.7

Kami akan menghitung

X = np.zeros((N, 2))

np.random.seed(seed=42)
#    ,    №1
X[:N1, 0] = np.random.normal(loc=mu_samples_1_1, scale=sigma_samples_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_samples_1_2, scale=sigma_samples_1_2, size=N1)
#    ,    №2
X[N1:N, 0] = np.random.normal(loc=mu_samples_2_1, scale=sigma_samples_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_samples_2_2, scale=sigma_samples_2_2, size=N2)

#   
m = X.shape[1]

#   
n = X.shape[0]

#        (   )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#     ,    №1  №2
w = np.array([float(1./k), float(1./k)])

np.random.seed(seed = None)
#        (   )
mu  = np.array(
    (np.mean(X[np.random.choice(n, n/k)], axis = 0), np.mean(X[np.random.choice(n, n/k)], axis = 0)))
# mu = np.array(([mu_samples_1_1, mu_samples_1_2],[mu_samples_2_1, mu_samples_2_2]))

#         (    )
sigma = np.array(([1., 0.],[0., 1.], [1., 0.],[0., 1.]))
# sigma = np.array(([sigma_samples_1_1, 0.],[0., sigma_samples_1_2], [sigma_samples_2_1, 0.],[0., sigma_samples_2_2]))
sigma = sigma.reshape(k, m, m)

#    EM-
steps = 15
#   EM-
for i in range(steps):
    proba_xi = e_step(X, k, m, n, w, mu, sigma)
    w, mu, sigma = m_step(X, proba_xi,N)
    X1_new_ind, X2_new_ind, X_answers = x_new(proba_xi)
    print ' №', i+1
    print
    print '   '
    print mu
    print
    print '   '
    print sigma
    print
    print '   '
    print round(accuracy_score(y, X_answers),3)
    
    plt.figure(figsize=(16, 6))  
    plt.plot(
        X[X1_new_ind,0], X[X1_new_ind,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        X[X2_new_ind,0], X[X2_new_ind,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'purple', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()

Bagan 3.1 "Distribusi bagian pada setiap iterasi dari algoritma-EM"

Semua iterasi

Apa yang terjadi pada akhirnya? Dan ternyata yang berikut, tidak memiliki informasi awal tentang kemungkinan hipotesis untuk ekstraksi bagian yang diproduksi pada mesin No. 1 -

w_{1}

$w_1$ atau nomor mesin 2

w_{2}

$w_2$ , ekspektasi matematis dan standar deviasi variabel acak -

μ

$\mu$ dan

σ

$\sigma$ , kami mencapai kualitas yang sebanding dengan contoh sebelumnya, di mana semua parameter di atas diketahui.

Bandingkan parameter yang kami atur saat menghasilkan set produk, yang disebut parameter benar dan parameter yang ditentukan algoritma EM untuk kami:

Parameter sejati

μ_{11} = 64, μ_{12} = 14, σ_{11} = 3.5, σ_{12} = 1.

$\mu_{11} = 64, \mu_{12} = 14, \sigma_{11} = 3.5, \sigma_{12} = 1.$

μ_{21} = 52, μ_{12} = 9.5, σ_{21} = 2.0, σ_{22} = 0.7

$\mu_{21} = 52, \mu_{12} = 9.5, \sigma_{21} = 2.0, \sigma_{22} = 0.7$

Parameter khusus

μ_{11} \approx 64.03, μ_{12} \approx 13.99, σ_{11} \approx 3.44, σ_{12} \approx 1.23

$\mu_{11} \approx 64.03, \mu_{12} \approx 13.99, \sigma_{11} \approx 3.44, \sigma_{12} \approx 1.23$

μ_{21} \approx 52.06, μ_{12} \approx 9.54, σ_{21} \approx 1.65, σ_{22} \approx 0.77

$\mu_{21} \approx 52.06, \mu_{12} \approx 9.54, \sigma_{21} \approx 1.65, \sigma_{22} \approx 0.77$

Artikel ini berakhir. Sebagai kesimpulan, saya akan menambahkan bahwa untuk keperluan analisis data, tidak masuk akal untuk menemukan kembali roda dan menulis algoritma EM sendiri, Anda cukup menggunakan fungsi perpustakaan sklearn yang disediakan.

Mari kita lihat bagaimana pustaka GaussianMixture sklearn bekerja contoh kita.

Modul impor

from sklearn.mixture import GaussianMixture as GMM

Tentukan parameternya

#   
#      №1 ( 1)
N1 = 6000
#      №2 ( 2)
N2 = 4000
#       
N = N1 + N2

#    №1
mu_samples_1_1 = 64.
#    №1
mu_samples_1_2 = 14.

#    №2
mu_samples_2_1 = 52.
#    №2
mu_samples_2_2 = 9.5

#      №1
sigma_samples_1_1 = 3.5
#      №1
sigma_samples_1_2 = 1.

#      №2
sigma_samples_2_1 = 2.
#      №2
sigma_samples_2_2 = 0.7


#  
k = 2

X = np.zeros((N, 2))

np.random.seed(seed=42)
#    ,    №1
X[:N1, 0] = np.random.normal(loc=mu_samples_1_1, scale=sigma_samples_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_samples_1_2, scale=sigma_samples_1_2, size=N1)
#    ,    №2
X[N1:N, 0] = np.random.normal(loc=mu_samples_2_1, scale=sigma_samples_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_samples_2_2, scale=sigma_samples_2_2, size=N2)

#        (   )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

Jalankan algoritma EM

np.random.seed(1)
model = GMM(n_components=k, covariance_type='full')
model.fit(X)


temp_predict_X = model.predict(X)
X_answers = []
for i in range(X.shape[0]):
    if temp_predict_X[i] == 0:
        X_answers.append(1)
    else:
        X_answers.append(2)
        

print '   '
print round(accuracy_score(y, X_answers),3)

Bekerja dengan keras! Tidak ada kesalahan! Seperti yang Anda harapkan, pustaka GaussianMixture skussian lebih cepat dan lebih baik daripada algoritma kami. Namun, tujuan artikel ini, seperti yang kita ingat, adalah untuk memahami cara kerja algoritma EM. Saya harap tujuannya telah tercapai :)

Sumber informasi

Sastra

1) Statistik terapan: Klasifikasi dan pengurangan dimensi, S.A. Ayvazyan, V.M. Buchstaber, I.S. Enyukov, L.D. Meshalkin, Moskow, Keuangan dan Statistik,

Sumber daya Internet 1989

1) Artikel 'Tutorial Lembut Algoritma EM dan Penerapannya pada Estimasi Parameter untuk Campuran Gaussian dan Model Markov Tersembunyi', International Computer Science Institute, Jeff A. Bilmes, 1998

2) Presentasi 'Model Campuran & Algoritma EM', NY University, David Sontag

3) Kursus 'Cari struktur data'

4) Formula untuk probabilitas penuh dan formula Bayes

5) Fungsi 'sklearn.mixture.GaussianMixture'

6) Distribusi normal, Wikipedia

7) Multidimensional distribusi normal, Wikipedia

8)github.com/diefimov/MTH594_MachineLearning/blob/master/ipython/Lecture10.ipynb

9) Artikel tentang algoritma EM

10) Artikel lain tentang algoritma EM

11) Github github.com/AlexanderPetrenko83/Articles

Kami membongkar algoritma-EM menjadi batu bata kecil

Algoritma EM secara umum

Ingat formula Bayes

Teorema Contoh Satu Oleh Bayes

Contoh kedua untuk teorema Bayes menggunakan parameter distribusi normal dari variabel acak: ekspektasi matematis μ\mudan standar deviasi σ\sigma

Contoh ketiga. Kasus dua dimensi

Memenuhi algoritma EM

Sumber informasi

More articles:

Contoh kedua untuk teorema Bayes menggunakan parameter distribusi normal dari variabel acak: ekspektasi matematis $\mu$ dan standar deviasi $\sigma$