"Introduction à l'apprentissage automatique par inférence bayésienne" Inférence approximative du modèle mixte de Poisson implémenté uniquement avec Python numpy

J'ai lu le mois dernier ["Introduction to Machine Learning by Bayesian Inference"](https://www.amazon.co.jp/ Machine Learning Startup Series-Introduction to Machine Learning by Bayesian Inference-KS Information Science Special Book-Suyama-Atsushi / dp / J'ai implémenté 4061538322) en Python.

Le contenu est le chapitre 4, «4.3 Raisonnement dans un modèle mixte de Poisson». Un livre vert avec un motif de fleurs. Ceci est une introduction ... L'apprentissage automatique est trop profond. Lol

Créer des données

Cette fois, nous allons créer une distribution de Poisson bimodale sous forme de données multimodales unidimensionnelles.

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

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

#Créer des exemples de données
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)

#Combiner des données
X = np.concatenate([x1, x2])

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

Le résultat est le suivant.

• Je pense que la forme du graphique changera légèrement en raison du caractère aléatoire de la distribution de Poisson. Pour expliquer brièvement Moyenne: 30, Nombre de données: 300 Moyenne: 50, Nombre de données: 200 La distribution est un mélange des deux distributions de Poisson.

C'est une coïncidence complète, mais ce graphique ne semble pas intuitivement bimodal, c'est donc une preuve parfaite de l'efficacité de l'inférence bayésienne!

Créer une fonction logumexp

Tout d'abord, comme préparation, créez une fonction qui calcule logsumexp.

… Et ici, «logsumexp» n'est pas apparu dans le livre, non? Je t'y pensais! C'est vrai. Ça ne sort pas. Lol

Cependant, ce "logsumexp" joue un rôle important dans cet algorithme!

Premièrement, cette fonction est utilisée dans l'équation (4.38) du livre. Puisqu'il y a une expression conditionnelle de η dans la ligne inférieure, logsumexp est requis.

En fait, si vous calculez normalement η selon la formule de la ligne supérieure, η ne satisfait pas la formule conditionnelle de la ligne inférieure, vous devez donc la "normaliser". Aussi, dans cette normalisation, il s'agit d'une opération pour "aligner les valeurs totales sur 1", il est donc nécessaire de "diviser chaque valeur par la valeur totale". Par conséquent, la formule décrite ci-dessus est transformée comme suit. La deuxième formule est exp⁡(logx) = x exp⁡(-logx) = -x Il peut être transformé à partir de la formule de.

En outre, la troisième étape est facile si vous utilisez la loi exponentielle.

• Exp (n) signifie "e à la nième puissance".

En conséquence, un terme de normalisation a été ajouté à la fin de la troisième ligne. Si vous regardez attentivement ce terme de normalisation, il contient "log", "sum [= Σ]" et "exp [= η]"!

Et lorsqu'il s'agit de cette "exp de logsum", il semble que cela puisse provoquer un débordement ou un sous-dépassement ...

Par conséquent, afin d'éviter les débordements et sous-débordements à l'avance, nous avons besoin d'une fonction logsumexp à implémenter à l'avenir!

Cela fait longtemps, mais la fonction elle-même est facile, alors implémentons-la immédiatement.

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

Je ne peux pas l'expliquer davantage, donc veuillez consulter l'article "Distribution gaussienne mixte et logsumexp" pour plus de détails.

Inférence approximative du modèle mixte de Poisson par échantillonnage de Gibbs

Nous allons enfin implémenter l'algorithme d'ici!

Cette fois, ["Introduction à l'apprentissage automatique par l'inférence bayésienne"](https://www.amazon.co.jp/ Série de démarrage d'apprentissage automatique - Introduction à l'apprentissage automatique par Bayesian Inference-KS Livre spécial sur les sciences de l'information-Suyama-Atsushi / dp / Implémenter "Algorithm 4.2 Gibbs Sampling for Poisson Mixed Model" décrit dans 4061538322) basé sur Python numpy.

#Préparez une liste pour l'échantillonnage
sample_s = []
sample_lambda = []
sample_pi = []

#Définition des constantes(Nombre de répétitions,Le nombre d'échantillons,Nombre de clusters)
MAXITER = 100
N = len(X)
K = 2

#Valeur initiale du paramètre
init_Gam_param_a = 1
init_Gam_param_b = 1
init_Dir_alpha = np.ones(K)

# λ,Réglage de la valeur initiale de π
Lambda = np.ones(K)
Pi = np.ones(K)

#Normalisé selon la condition de π(Divisez chaque valeur par la valeur totale)
if np.sum(Pi) != 1:
    Pi = Pi / np.sum(Pi)
    
#Échantillonnage répété pour le nombre de MAXITER
for i in range(MAXITER):
    
    #Préparer la base de données de s
    s = np.zeros((N, K))
    
    #Calculer η pour échantillonner 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))
    
    #Normaliser η(Utilisez logsumexp pour éviter les débordements et sous-débordements)
    logsumexp_eta = -1 * log_sum_exp(log_eta)
    eta = np.exp(log_eta + logsumexp_eta)

    #Échantillon s de la distribution des catégories avec η comme paramètre
    for n in range(N):
        s[n] = np.random.multinomial(1, eta[n])   
    #Ajouter à la liste d'échantillons
    sample_s.append(np.copy(s))
    
    
    #a pour échantillonner λ,Calculer 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/Échantillon λ de la distribution gamma avec b comme paramètre
    Lambda = np.random.gamma(Gam_param_a, 1/Gam_param_b)
    #Ajouter à la liste d'échantillons
    sample_lambda.append(np.copy(Lambda))
    
        
    #Calculer α pour échantillonner π
    Dir_alpha = np.sum(s, axis=0) + init_Dir_alpha
    
    #Échantillon π de la distribution Diricle avec α comme paramètre
    Pi = np.random.dirichlet(Dir_alpha)
    #Ajouter à la liste d'échantillons
    sample_pi.append(np.copy(Pi))
    
#Clusters sans ordre particulier(Parce que l'ordre n'est pas défini)
sample_s_ndarray = np.array(sample_s)
sample_lambda_ndarray = np.array(sample_lambda)
sample_pi_ndarray = np.array(sample_pi)

Il est essentiellement implémenté selon le livre, mais notez les points suivants.

-Mettre à jour chaque paramètre η, a, b, α avant d'échantillonner chacun des s, λ, π ・ Π doit être normalisé lorsqu'il s'agit de la valeur initiale (il suffit de diviser chaque valeur par la valeur totale) ・ Η est normalisé en combinant avec logsumexp dans np.exp () ・ Dans le livre, l'instruction for de N et K est utilisée, mais si vous effectuez un calcul matriciel, l'instruction for n'est requise que lors de l'échantillonnage de s. -Les valeurs échantillonnées sont ajoutées à chaque liste correspondante à chaque fois

En outre, la valeur initiale de chaque paramètre peut être 1 ou quelque chose. (Attention à la normalisation uniquement pour π!)

Confirmation du résultat de l'échantillonnage

Vérifions les résultats dans l'ordre!

… Et avant cela, notez que les clusters sont en panne. Cette fois, les résultats ont été obtenus par hasard dans l'ordre des paramètres de cluster, mais l'ordre de configuration et les résultats de cluster obtenus peuvent ne pas correspondre. Cependant, soyez assuré qu'il n'y a pas de problème particulier dans ce cas également.

Commençons par λ!

#Valeur moyenne pour chaque cluster
ave_Lambda = list(np.average(sample_lambda_ndarray, axis=0))

print(f'prediction: {ave_Lambda}')

Le résultat est le suivant. prediction: [29.921538459827033, 49.185569726045905]

λ correspond à la valeur moyenne de chaque cluster. Lors de la création des données Lambda_x1 = 30 Lambda_x2 = 50 Je l'ai mis, donc c'est une assez bonne précision!

Ensuite, regardons π.

#Pourcentage d'échantillons en grappes dans toutes les données
ave_Pi = list(np.average(sample_pi_ndarray, axis=0))

all_samples = samples_x1 + samples_x2

print(f'prediction: {ave_Pi}')

Les résultats sont les suivants. prediction: [0.5801320180878531, 0.4198679819121469]

Concernant le nombre de données samples_x1 = 300 samples_x2 = 200 Je l'ai mis. Diviser chaque valeur par le nombre total de données, 500, est [0,6, 0,4], donc elles sont presque les mêmes!

Enfin, vérifiez s et terminez! Je vous remercie pour votre travail acharné.

#Nombre d'échantillons dans chaque cluster
sum_s = np.sum(np.sum(sample_s_ndarray, axis=0), axis=0) / MAXITER

print(f'prediction: {sum_s}')

Le résultat est, prediction: [291.18 208.82] est devenu.

Le nombre d'échantillons obtenus étant doublé du nombre de MAXITER, il peut être obtenu en divisant l'échantillon total de chaque grappe par MAXITER.

Références / Liens de référence

[«Introduction à l'apprentissage automatique par l'inférence bayésienne»](https://www.amazon.co.jp/ Série de démarrage d'apprentissage automatique - Introduction à l'apprentissage automatique par l'inférence bayésienne - Livre spécial sur la science de l'information KS - Suyama-Atsushi / dp / 4061538322) "Distribution gaussienne mixte et exp. De la somme des journaux"