"Einführung in das maschinelle Lernen durch Bayes'sche Inferenz" Ungefähre Inferenz eines Poisson-Mischmodells, das nur mit Python-Numpy implementiert wurde

Ich habe letzten Monat ["Einführung in das maschinelle Lernen durch Bayesianische Inferenz"] gelesen (https://www.amazon.co.jp/ Startup-Serie zum maschinellen Lernen - Einführung in das maschinelle Lernen durch Bayesianische Inferenz-KS Informationswissenschaft Spezialbuch-Suyama-Atsushi / dp / Ich habe 4061538322) in Python implementiert.

Der Inhalt ist Kapitel 4, "4.3 Argumentation in einem Poisson-Mischmodell". Ein grünes Buch mit einem Blumenmuster. Dies ist eine Einführung ... Maschinelles Lernen ist zu tief. Lol

Daten erstellen

Dieses Mal werden wir eine bimodale Poisson-Verteilung als multimodale eindimensionale Daten erstellen.

import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

plt.figure(figsize=(12, 5))
plt.title('Poisson distribution', fontsize=20)

#Beispieldaten erstellen
Lambda_x1 = 30
Lambda_x2 = 50
samples_x1 = 300
samples_x2 = 200
x1 = np.random.poisson(Lambda_x1, samples_x1)
x2 = np.random.poisson(Lambda_x2, samples_x2)

#Daten kombinieren
X = np.concatenate([x1, x2])

plt.hist(X, 30, density=True, label='all data')
plt.legend()
plt.show()

Das Ergebnis ist wie folgt.

• Ich denke, dass sich die Form des Diagramms aufgrund der Zufälligkeit der Poisson-Verteilung geringfügig ändern wird. Um es kurz zu erklären Durchschnitt: 30, Anzahl der Daten: 300 Durchschnitt: 50, Anzahl der Daten: 200 Die Verteilung ist eine Mischung der beiden Poisson-Verteilungen.

Dies ist ein völliger Zufall, aber dieses Diagramm scheint nicht intuitiv bimodal zu sein, daher ist es ein perfekter Beweis für die Wirksamkeit der Bayes'schen Inferenz!

Erstellen Sie eine logumexp-Funktion

Erstellen Sie zunächst als Vorbereitung eine Funktion, die logsumexp berechnet.

… Und hier ist „logsumexp“ nicht im Buch erschienen, oder? Ich dachte du da! Korrekt. Es kommt nicht heraus. Lol

Dieses "logsumexp" spielt jedoch eine wichtige Rolle in diesem Algorithmus!

Erstens wird diese Funktion in Gleichung (4.38) im Buch verwendet. Da es in der unteren Zeile einen bedingten Ausdruck von η gibt, ist logsumexp erforderlich.

Wenn Sie η normalerweise gemäß der Formel in der oberen Reihe berechnen, erfüllt η die bedingte Formel in der unteren Reihe nicht, sodass Sie sie "normalisieren" müssen. Bei dieser Normalisierung ist es auch eine Operation, "die Gesamtwerte auf 1 auszurichten", so dass es notwendig ist, "jeden Wert durch den Gesamtwert zu teilen". Daher wird die oben beschriebene Formel wie folgt transformiert. Die zweite Formel lautet exp⁡(logx) = x exp⁡(-logx) = -x Es kann aus der Formel von transformiert werden.

Auch die dritte Stufe ist einfach, wenn Sie das Exponentialgesetz verwenden.

• Exp (n) steht für "e zur n-ten Potenz".

Infolgedessen wurde am Ende der dritten Zeile ein Begriff für die Normalisierung hinzugefügt. Wenn Sie sich diesen Normalisierungsterm genau ansehen, enthält er "log", "sum [= Σ]" und "exp [= η]"!

Und wenn es um diese "logsum exp" geht, scheint es, dass ein Überlauf oder Unterlauf auftreten kann ...

Um Über- und Unterlauf im Voraus zu vermeiden, benötigen wir daher eine logsumexp-Funktion, die in Zukunft implementiert werden soll!

Es ist lange her, aber die Funktion selbst ist einfach, also implementieren wir sie sofort.

def log_sum_exp(X):
    max_x = np.max(X, axis=1).reshape(-1, 1)
    return np.log(np.sum(np.exp(X - max_x), axis=1).reshape(-1, 1)) + max_x

Ich kann es nicht weiter erklären. Weitere Informationen finden Sie im Artikel "Gemischte Gaußsche Verteilung und logsumexp".

Ungefähre Inferenz des Poisson-Mischmodells durch Gibbs-Probenahme

Wir werden den Algorithmus endlich von hier aus implementieren!

Diesmal ["Einführung in das maschinelle Lernen durch Bayesianische Inferenz"](https://www.amazon.co.jp/ Startup-Serie zum maschinellen Lernen - Einführung in das maschinelle Lernen durch Bayesianische Inferenz-KS Informationswissenschaft Spezialbuch-Suyama-Atsushi / dp / Implementieren Sie "Algorithmus 4.2 Gibbs Sampling für Poisson Mixed Model" (beschrieben in 4061538322) basierend auf Python numpy.

#Bereiten Sie eine Liste für die Probenahme vor
sample_s = []
sample_lambda = []
sample_pi = []

#Konstanten einstellen(Anzahl der Wiederholungen,Anzahl von Beispielen,Anzahl der Cluster)
MAXITER = 100
N = len(X)
K = 2

#Parameter Anfangswert
init_Gam_param_a = 1
init_Gam_param_b = 1
init_Dir_alpha = np.ones(K)

# λ,Anfangswerteinstellung von π
Lambda = np.ones(K)
Pi = np.ones(K)

#Normalisiert gemäß der Bedingung von π(Teilen Sie jeden Wert durch den Gesamtwert)
if np.sum(Pi) != 1:
    Pi = Pi / np.sum(Pi)
    
#Wiederholte Abtastung für die Anzahl der MAXITER
for i in range(MAXITER):
    
    #Bereiten Sie die Datenbank von s vor
    s = np.zeros((N, K))
    
    #Berechnen Sie η zu Probe s
    log_eta = np.dot(X.reshape(N, 1), np.log(Lambda.reshape(1, K))) - Lambda.reshape(1, K) + np.log(Pi.reshape(1, K))
    
    #Normalisieren Sie η(Verwenden Sie logsumexp, um Über- und Unterlauf zu verhindern)
    logsumexp_eta = -1 * log_sum_exp(log_eta)
    eta = np.exp(log_eta + logsumexp_eta)

    #Stichprobe s aus der Kategorieverteilung mit η als Parameter
    for n in range(N):
        s[n] = np.random.multinomial(1, eta[n])   
    #Zur Beispielliste hinzufügen
    sample_s.append(np.copy(s))
    
    
    #a, um λ abzutasten,Berechnen Sie b
    Gam_param_a = (np.dot(s.T, X.reshape(N, 1)) + init_Gam_param_a).T[0]
    Gam_param_b = np.sum(s, axis=0).T + init_Gam_param_b
    
    # a, 1/Probe λ aus der Gammaverteilung mit b als Parameter
    Lambda = np.random.gamma(Gam_param_a, 1/Gam_param_b)
    #Zur Beispielliste hinzufügen
    sample_lambda.append(np.copy(Lambda))
    
        
    #Berechnen Sie α, um π abzutasten
    Dir_alpha = np.sum(s, axis=0) + init_Dir_alpha
    
    #Probe π aus der Richtungsverteilung mit α als Parameter
    Pi = np.random.dirichlet(Dir_alpha)
    #Zur Beispielliste hinzufügen
    sample_pi.append(np.copy(Pi))
    
#Cluster in keiner bestimmten Reihenfolge(Weil die Reihenfolge nicht definiert ist)
sample_s_ndarray = np.array(sample_s)
sample_lambda_ndarray = np.array(sample_lambda)
sample_pi_ndarray = np.array(sample_pi)

Es wird grundsätzlich gemäß dem Buch implementiert, beachten Sie jedoch die folgenden Punkte.

• Aktualisieren Sie jeden Parameter η, a, b, α, bevor Sie jeweils s, λ, π abtasten ・ Π muss normalisiert werden, wenn es sich um den Anfangswert handelt (dividieren Sie einfach jeden Wert durch den Gesamtwert). ・ Η wird durch Kombinieren mit logsumexp in np.exp () normalisiert. ・ In dem Buch wird die for-Anweisung von N und K verwendet. Wenn jedoch die Matrix berechnet wird, ist die for-Anweisung nur erforderlich, wenn s abgetastet wird.
• Zu jeder entsprechenden Liste werden jedes Mal abgetastete Werte hinzugefügt

Der Anfangswert jedes Parameters kann auch 1 oder so sein. (Achten Sie auf die Normalisierung nur für π!)

Bestätigung des Probenahmeergebnisses

Lassen Sie uns die Ergebnisse der Reihe nach überprüfen!

… Und vorher beachten Sie, dass die Cluster nicht in Ordnung sind. Dieses Mal wurden die Ergebnisse zufällig in der Reihenfolge der Clustereinstellungen erhalten, aber die Reihenfolge der Einstellungen und die erhaltenen Clusterergebnisse stimmen möglicherweise nicht überein. Bitte seien Sie jedoch versichert, dass auch in diesem Fall kein besonderes Problem vorliegt.

Beginnen wir mit λ!

#Durchschnittswert für jeden Cluster
ave_Lambda = list(np.average(sample_lambda_ndarray, axis=0))

print(f'prediction: {ave_Lambda}')

Das Ergebnis ist wie folgt. prediction: [29.921538459827033, 49.185569726045905]

λ entspricht dem Durchschnittswert jedes Clusters. Beim Erstellen der Daten Lambda_x1 = 30 Lambda_x2 = 50 Ich habe es eingestellt, also ist es ziemlich genau!

Als nächstes schauen wir uns π an.

#Prozentsatz der Cluster-Stichproben in allen Daten
ave_Pi = list(np.average(sample_pi_ndarray, axis=0))

all_samples = samples_x1 + samples_x2

print(f'prediction: {ave_Pi}')

Die Ergebnisse sind wie folgt. prediction: [0.5801320180878531, 0.4198679819121469]

In Bezug auf die Anzahl der Daten samples_x1 = 300 samples_x2 = 200 Ich habe es eingestellt. Teilen Sie jeden Wert durch die Gesamtzahl der Daten, 500, ist [0,6, 0,4], also sind sie fast gleich!

Zum Schluss s überprüfen und fertig! Danke für deine harte Arbeit.

#Anzahl der Proben in jedem Cluster
sum_s = np.sum(np.sum(sample_s_ndarray, axis=0), axis=0) / MAXITER

print(f'prediction: {sum_s}')

Ergebnis ist, prediction: [291.18 208.82] ist geworden.

Da die Anzahl der erhaltenen Proben durch die Anzahl der MAXITER verdoppelt wird, kann sie erhalten werden, indem die Gesamtprobe jedes Clusters durch MAXITER dividiert wird.

Referenzen / Referenzlinks

["Einführung in das maschinelle Lernen durch Bayes'sche Inferenz"](https://www.amazon.co.jp/ Startup-Serie zum maschinellen Lernen - Einführung in das maschinelle Lernen durch Bayes'sche Inferenz-KS Informationswissenschaft Spezialbuch-Suyama-Atsushi / dp / 4061538322) "Gemischte Gaußsche Verteilung und logsum exp"