Introduction à la modélisation statistique bayésienne avec python ~ Essai de régression linéaire avec MCMC ~

introduction

L'approche complète basée sur les données est-elle correcte? Utilisez plus de savoir-faire humain pour la modélisation!

«Je pense que l'approche du big data et du machine learning de nos jours est d'apprendre toutes les règles à partir des données. «Cependant, en fait, il est rare que ** une grande quantité de données complètes sur tous les événements soit accumulée **.

• Par conséquent, je voudrais prêter attention aux ** savoir-faire (connaissances) qui n'ont pas été convertis en données mais qui ont été accumulés à partir de nombreuses années d'expérience et qui existent dans l'esprit humain ** ―― En utilisant ces savoir-faire, il est possible de compléter des données insuffisantes par des statistiques bayésiennes (modélisation de Bayes). `` En raison des limites d'une approche basée sur les données pour tous les événements, les statistiques bayésiennes attirent à nouveau l'attention, et récemment, des termes tels que ** bayésien deep learning ** et ** bayésien machine learning ** sont apparus. Masu --Ainsi, dans cet article, je vais vous expliquer les bases des statistiques bayésiennes en utilisant la bibliothèque python (pyMC). «Le concept des statistiques bayésiennes est difficile à« comprendre »pour les débutants, mais il est facile d'obtenir une image lors de la saisie à partir de« utiliser »en utilisant des données réelles.

Statistiques bayésiennes (de base)

Ici, je n'introduirai que l'idée minimale (image) et le théorème de base. ――Il existe de nombreux livres faciles à comprendre pour des explications précises et détaillées des statistiques bayésiennes, veuillez donc vous y référer.

Pensée bayésienne

Les statistiques bayésiennes déduisent la «cause» du «résultat»

«Habituellement, les personnes qui commencent à étudier les statistiques commencent souvent par des statistiques de la théorie des fréquences, mais je pense que c'est l'une des raisons pour lesquelles les statistiques bayésiennes sont difficiles à comprendre.

• Dans les statistiques de la théorie des fréquences, la vraie valeur de la population est la «cause», et les données observées sont considérées comme le «résultat». ――D'autre part, les statistiques bayésiennes estiment la population (cause) à partir des données observées (résultats).
• Par conséquent, les statistiques bayésiennes traitent la population comme des variables (et non des valeurs ou des constantes qui changent en fonction des données observées).

Dans les statistiques bayésiennes, les variables sont subjectivement «déterminées provisoirement» et corrigées (mises à jour) avec les données observées.

• Lors de l'estimation des variables (paramètres) de «cause», les statistiques bayésiennes utilisent les concepts de «pré-distribution» et de «distribution postérieure».
• La distribution a priori peut être établie sur la base des connaissances que les humains connaissent déjà (c'est la raison pour laquelle le savoir-faire de l'esprit humain peut être utilisé en utilisant les statistiques bayésiennes). ―― La distribution du résultat de la correction (mise à jour) de cette distribution antérieure avec les données observées est appelée «distribution postérieure». --Dans les statistiques bayésiennes, le processus de modification de la distribution basée sur ces données est appelé "mise à jour des probabilités".

Théorème de Bayes

• Ici, je voudrais introduire le théorème le plus important de la statistique bayésienne, le "théorème de bases". ――Cependant, il n'y a pas de nouvelle information car les contenus ci-dessus sont simplement exprimés dans des formules mathématiques. ――Ce théorème lui-même peut être facilement dérivé de la transformation d'équation de la probabilité conditionnelle, donc si vous êtes intéressé, veuillez le vérifier. ――Je pense que vous devriez au moins obtenir l'image que la distribution antérieure de l'équation de droite est mise à jour avec les données et la distribution postérieure de l'équation de gauche est dérivée.
• Les statistiques de Bayes utilisent essentiellement ce théorème pour déduire des valeurs approximatives par analytiquement ou (lorsqu'il ne peut pas être résolu analytiquement) le paramètre θ en utilisant des simulations de nombres aléatoires et des méthodes d'optimisation.

A part: Interpréter le théorème de Bayes à partir de termes qui apparaissent dans l'apprentissage automatique

――Ce qui suit est un aparté, donc si vous n'êtes pas intéressé, veuillez l'ignorer. ---- P (X | θ) est appelé vraisemblance --La ** méthode d'estimation la plus probable **, qui apparaît également dans l'apprentissage automatique, recherche θ qui maximise ce P (X | θ).

• En d'autres termes, la méthode d'estimation la plus probable ne prend pas en compte la distribution a priori et peut être considérée comme une procédure qui recherche simplement θ qui maximise P (X | θ). Il s'agit d'une approche qui adapte parfaitement les paramètres aux données disponibles, donc elle provoque un surentraînement (surapprentissage).
• De plus, la méthode de recherche de paramètres appelée ** estimation MAP ** recherche θ qui maximise P (X | θ) × P (θ). ――Cela peut être considéré comme une recherche de paramètre qui considère la condition de contrainte (régularisation) de P (θ), qui est une distribution a priori. ――Comme la distribution a priori est également prise en compte, elle ne surajuste pas autant les données disponibles que l'estimation la plus probable, mais le calcul devient compliqué car la fonction de densité de probabilité de la distribution a priori doit être calculée.
• Dans le théorème de Bayes, la probabilité postérieure est calculée en considérant la vraisemblance périphérique (probabilité que les données puissent être obtenues en moyenne) dans l'estimation MAP.

Méthode de calcul

• Ici, nous allons expliquer comment calculer les paramètres. ――Il existe différentes méthodes de calcul, mais ici nous allons introduire trois méthodes typiques.

1. Distribution a priori conjuguée naturelle --Il existe des combinaisons de distributions de probabilité qui ont les mêmes distributions a priori et postérieur (par exemple, distribution de Bernoulli et distribution bêta). ――Il s'agit d'une méthode pour calculer analytiquement des paramètres en tirant parti de ses caractéristiques. ――Cependant, la distribution est très limitée pour satisfaire la condition de combinaison, elle est donc rarement utilisée dans la pratique à l'heure actuelle.

• Dans les temps modernes, lorsque les ordinateurs se sont développés, l'inférence de valeur approximative à l'aide des simulations et optimisations de nombres aléatoires suivantes est la solution principale.

2. MCMC

• Le nom officiel est une méthode appelée méthode de Monte Carlo en chaîne de Markov (en bref, c'est une méthode pour générer des "nombres aléatoires" par simulation). --Cependant, MCMC peut générer des nombres aléatoires qui suivent la distribution postérieure, ce qui évite des calculs intégraux compliqués. --Il existe divers algorithmes dans MCMC (par exemple, échantillonnage de Gibbs, algorithme M-H).
• En python, il peut être exécuté en utilisant ** pyMC **

3. Raisonnement variable ――Il s'agit d'une méthode pour rechercher une distribution proche de la distribution postérieure par la méthode du gradient. --Définir la fonction de perte (distance) de la distribution approximative et de la distribution postérieure avec divergence KL, et rechercher la valeur approximative en minimisant la fonction de perte (c'est une approche souvent utilisée dans les réseaux de neurones).

• En python, il peut être exécuté en utilisant ** pyro **

――Dans cet article, nous présenterons l'estimation des paramètres par MCMC.

Modèle de régression démo-linéaire-

HardC'est difficile à comprendre dans les phrases, alors je vais l'essayer avec python immédiatement.

• Ici, je vais résoudre un modèle de régression linéaire très simple avec des statistiques bayésiennes. «Honnêtement, ce problème n'a pas besoin d'être des statistiques bayésiennes, mais si vous maintenez la procédure, vous pouvez résoudre des modèles complexes avec la même approche. ―― Tout d'abord, créez des exemples de données à utiliser dans cette démo

Créer des exemples de données

#Créer des exemples de données
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)

#Exemples de paramètres de génération de données
N = 100
alpha_real = 2.5
beta_real = 0.9
eps_real = np.random.normal(0, 0.5, size=N)

#modèle: y = α + βx + Noise
x = np.random.normal(10, 1, N)
y_real = alpha_real + beta_real *x
y = y_real + eps_real

#Visualisation de données
plt.plot(x, y, '.')
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x, y_real, 'orange')
plt.show()

• Le paramètre problématique est que vous souhaitez estimer le modèle linéaire (ligne orange) lorsque les données (point bleu) sont fournies.
• Le modèle linéaire est défini par y = α + β * x, et α = 2,5 et β = 0,9 sont les valeurs que vous voulez déduire à l'aide des statistiques bayésiennes.

Modèle de régression linéaire avec MCMC

• Tout d'abord, installez pyMC

conda install pymc3==3.6

• Pour le moment (à partir de novembre 2020), la dernière version est la 3.9, mais il y avait un problème avec certaines fonctions de mon environnement, donc je rétrograde et je l'utilise. ――Veuillez installer tout en modifiant cette zone selon les besoins

Procédure de base de pyMC

1. Définition du modèle et exécution MCMC
2. Confirmation / évaluation de la distribution postérieure de chaque paramètre
3. Inférence utilisant la distribution postérieure (visualisation de la prédiction / intervalle de confiance)

Comment écrire un modèle avec pyMC

--Dans pyMC, vous devez décrire le modèle comme suit ――J'ai écrit diverses choses, mais le point est ** 1. Définition du modèle, 2. Calcul des probabilités, 3. Exécution de MCMC **. -Si vous cochez Official Document, diverses autres options (réglage de la valeur initiale et autres algorithmes MCMC) sont introduites. ―― Tout d'abord, il est préférable de comprendre la méthode de description minimale du modèle ci-dessous, puis d'essayer d'autres options si nécessaire.

#Dans PyMC, le modèle est défini par with ward
with pm.Model() as model:
    # 1.Définition de la distribution antérieure (définissez la variable que vous souhaitez utiliser comme paramètre)
    '''
Diverses distributions de probabilité sont préparées par pyMC, et la distribution normale est utilisée ici.
Définissez le nom de la distribution sous forme de chaîne de caractères, puis définissez les paramètres de distribution précédents.
    '''
    α = pm.Normal('alpha', mu=0, sd=10)#Distribution normale (moyenne=0,écart-type=10)
    β = pm.Normal('beta', mu=0, sd=1)#Distribution normale (moyenne=0,écart-type=1)
    noise = pm.Normal('noise', mu=0, sd=1)
    # 2.Calcul de la probabilité
    '''
La mise à jour des probabilités est effectuée en donnant la valeur observée à
    '''
     y_pred = pm.Normal('y_pred', mu=α*x+β, sd=noise, observed=y)
    # 3.Exécutez MCMC
    '''
Il existe différents paramètres tels que les valeurs initiales et les algorithmes, mais c'est basique et polyvalent."NUTS"Utilisez l'algorithme
    draws=Nombre de nombres aléatoires à générer, chains=Nombre de processus parallèles
    '''
    trace = pm.sample(draws=5000, chains=1) #5000 échantillons*1 parallèle=Générer 5000 nombres aléatoires

Effectuer une modélisation bayésienne

• Effectuons une modélisation bayésienne selon la procédure ci-dessus.

import pymc3 as pm    
    
#modèle
with pm.Model() as model:
    # 1.Distribution antérieure
    alpha = pm.Normal('alpha', mu=0, sd=10)#0 en moyenne pour α,En supposant une distribution normale avec un écart type de 10
    beta = pm.Normal('beta', mu=0, sd=1)#0 en moyenne pour β,En supposant une distribution normale avec un écart type de 1
    epsilon = pm.HalfCauchy('epcilon', 5)#En supposant une distribution semi-cauchy avec une fréquence maximale de 5 pour le bruit
    # 2.Calcul de la probabilité
    y_pred = pm.Normal('y_pred', mu=alpha+beta*x, sd=epsilon, observed=y)
    # 3.Exécutez MCMC
    trace = pm.sample(11000, chains=1)

• Ceci termine la simulation de la distribution postérieure des paramètres. ――Je vérifierai immédiatement le résultat de la simulation.

#Confirmation / évaluation de la distribution postérieure de chaque paramètre (graphique)
trace_n = trace[1000:]#Jeter les 1000 premiers cas (voir ci-dessous:"Brûler"）
pm.traceplot(trace_n)

#Confirmation / évaluation de la distribution postérieure de chaque paramètre (statistiques)
pm.summary(trace_n)

• Comment lire le résultat de sortie
• Le graphique (colonne de gauche) montre la distribution de chaque paramètre (fonction de densité de probabilité).
• La distribution est centrée sur α = environ 2,64 et β = environ 0,88, et vous pouvez voir que des paramètres proches de la bonne réponse (α = 2,5, β = 0,9) peuvent être estimés.
• Le graphique (colonne de droite) montre l'état de convergence de la simulation de nombres aléatoires, et cette fois, il converge vers une distribution stable, mais si cela ne converge pas bien, un prétraitement tel que la pré-distribution, ses paramètres, la standardisation des quantités de caractéristiques, etc. Tu as besoin de faire

--Brûler

• Dans MCMC, il faut un certain temps pour obtenir un échantillon de la distribution cible, il est donc préférable de ne pas utiliser les données au début de la génération de l'échantillon. ――Par conséquent, nous évaluons la distribution postérieure des paramètres à l'exception des 1000 premiers échantillons générés.

Enfin, utilisez les résultats de la simulation pour tracer une droite de régression et évaluer les résultats d'inférence du modèle.

# α,Acquisition de la valeur β (la valeur moyenne des résultats de la simulation est utilisée comme valeur représentative de chaque paramètre)
alpha_m = trace_n['alpha'].mean()
beta_m = trace_n['beta'].mean()

#Visualisation
plt.plot(x, y, '.')
plt.plot(x, y_real, 'orange', label='True')
plt.plot(x, alpha_m + beta_m*x, 'red', label='Estimated')
plt.legend()

• Diagramme de modèle linéaire --Une ligne de régression est tracée en utilisant la ** valeur moyenne ** du résultat de la simulation de distribution postérieure comme valeur représentative du paramètre. ――Vous pouvez voir qu'une ligne droite similaire à la droite de régression correcte tracée dans "Création d'échantillons de données" peut être dessinée.

#Génération d'échantillons à partir de la distribution postérieure
ppc = pm.sample_posterior_predictive(trace_n, samples=1000, model=model)

#Visualisation du diagramme de dispersion et de la ligne de régression
plt.plot(x, y, '.')
plt.plot(x, alpha_m + beta_m*x)

#Intervalle de confiance
idx = np.argsort(x)
x_ord = x[idx]
# 50%Section HPD
sig0 = pm.stats.hpd(ppc['y_pred'], alpha=0.5)[idx]
# 95%Section HPD
sig1 = pm.stats.hpd(ppc['y_pred'], alpha=0.05)[idx]
#Représenter l'intervalle de confiance avec un remplissage
plt.fill_between(x_ord, sig0[:,0], sig0[:,1], color='gray', alpha=1)
plt.fill_between(x_ord, sig1[:,0], sig1[:,1], color='gray', alpha=0.5)

• Le modèle de Bayes (MCMC) peut avoir le résultat d'une simulation de nombres aléatoires de distribution postérieure ―― Cela vous permet d'exprimer l'intervalle de confiance à 50% et 95%.
• Dans le modèle Bayes, l'intervalle HPD (valeur que le paramètre peut prendre) est souvent utilisé.

Qu'est-ce que la section HPD?

• Également appelée section de crédit post-densité la plus élevée
• Ensemble de largeur d'intervalle fiable afin que la valeur la plus fréquente du paramètre soit toujours incluse quelle que soit la forme de la distribution C'est compliqué quand on regarde la formule, mais c'est une définition très naturelle quand on la vérifie avec une image (graphique). --Ce site est facile à comprendre.

enfin

Possibilités de modélisation bayésienne

Dans cet article, j'ai essayé de résoudre les bases de la statistique bayésienne et du modèle de régression linéaire en utilisant MCMC.

• ** La caractéristique de la modélisation bayésienne est que le résultat de l'estimation peut être exprimé avec une distribution de probabilité au lieu d'une valeur **. ――Lors de la prise de décision, il faut tenir compte non seulement du résultat de la prédiction mais aussi de la variation (fiabilité) de l'ensemble **, donc je pense que c'est très pratique.
• Cette fois, pour des raisons de simplicité, toutes les distributions antérieures sont supposées être des distributions normales, mais pyMC prend en charge diverses distributions de probabilités, et vous pouvez définir des distributions en fonction de vos connaissances antérieures. ――En plus du modèle de régression linéaire, vous pouvez aborder divers problèmes avec la même approche, comme résoudre le modèle de régression logistique avec le modèle Basian ou visualiser la fiabilité du résultat de l'estimation en modifiant simplement la fonction de lien. Masu «En ce sens, ** la modélisation Bayes peut être considérée comme une méthode de modélisation extrêmement flexible et polyvalente **. ―― D'un autre côté, MCMC prend beaucoup de temps à simuler lorsque le modèle devient compliqué (comme vous pouvez le voir en exécutant ce qui précède). Dans ce cas, vous pouvez également le résoudre par la méthode d'optimisation par inférence variable (j'ai essayé d'expliquer la solution par inférence variable par la bibliothèque python pyro, mais comme la quantité d'articles a augmenté, j'écrirai un autre article Je vais le créer)

Pour promouvoir l'utilisation des données

Dans cet article, j'ai essayé d'expliquer les bases des statistiques bayésiennes en utilisant python, qui est l'un des langages les plus utilisés dans la pratique de nos jours. «À l'ère actuelle où le DX et l'utilisation des données sont nécessaires, j'entends souvent les inquiétudes selon lesquelles" je veux promouvoir les efforts mais il n'y a pas de données appropriées !! " «Il est bien sûr important de créer une base de données à partir de zéro et de commencer à collecter des données, mais cela ne suffira pas à rattraper les entreprises mondiales qui ouvrent la voie dans les efforts d'utilisation des données. «Par conséquent, je pensais que la modélisation bayésienne ne pouvait être utilisée que par ** des entreprises japonaises qui n'ont pas de données mais un savoir-faire commercial abondant cultivé depuis de nombreuses années **, j'ai donc écrit cet article. J'espère que cela vous donnera un petit indice pour les entreprises qui souhaitent promouvoir l'utilisation des données mais qui n'ont pas réussi.