[PYTHON] Estimation bayésienne en utilisant le lancer de pièces comme exemple Voir la distribution postérieure pour chaque essai

introduction

Récemment, j'ai entendu dire que les statistiques bayésiennes sont chaudes, alors j'ai pensé que je devrais étudier "Introduction to Completely Self-Study Bayesian Statistics". C'est facile à comprendre à mi-parcours, mais j'ai senti que c'était difficile à comprendre dès la sortie de la distribution continue car c'est un style qui explique avec des chiffres sans utiliser autant que possible des formules mathématiques. Par conséquent, pour étudier, je vais utiliser le lancer de pièces de monnaie comme exemple simple et regarder la distribution postérieure à chaque fois que je le lance. Je ne comprends pas encore la méthode MCMC, je vais donc calculer la fonction de densité de probabilité en petits morceaux.

référence

Hiroyuki Kojima "Introduction to Bayesian Statistics", Diamond, novembre 2015

L'explication est basée sur le même chiffre que le livre ci-dessus, en utilisant le lancer de pièces comme exemple. Kitashukyo 108th Academic Education Practice Study Group samedi 26 janvier 2019 "La pensée de l'estimation bayésienne"

C'est le même sujet, mais c'est un contenu avancé calculé par la méthode MCMC à l'aide de pyMC. Je souhaite également pouvoir utiliser pyMC. Regardez la distribution postérieure du recto et du verso de la pièce pour chaque essai (PyMC)

Ceci est le premier chapitre de "Introduction aux statistiques bayésiennes". Le chiffre est souvent utilisé et il est facile à comprendre. Estimation Bayes pour comprendre clairement en 5 minutes

Estimation bayésienne

La force des statistiques bayésiennes est

  1. La propriété qu'il peut être deviné même s'il y a peu de données, et il devient plus précis car il y a plus de données.
  2. Fonction d'apprentissage qui met à jour automatiquement les estimations en réponse aux informations entrantes instantanément Il semble que ce soit dans [^ 1].

L'estimation bayésienne est basée sur le théorème bayésien. $p(\theta|x)=p(\theta) \frac{p(x|\theta)}{\int p(x|\theta)p(\theta)d\theta}$ ici\thetaEst un paramètre,p(\theta)Est une distribution antérieure,p(x|\theta)Est\thetaQuandxProbabilité conditionnelle de(Responsabilité)、p(\theta|x)Est事後分布です。 Dans l'exemple du lancer de pièces, $ \ theta $ est l'avant ou l'arrière, $ p (\ theta) $ est la distribution de probabilité de $ \ theta $, $ x $ est les données indiquant si la pièce est avant ou arrière, et le dénominateur est Représente la probabilité d'obtenir $ x $. La probabilité originale que $ \ theta $ se produise était la pré-probabilité $ p (\ theta) $, mais avec l'information que $ x $ se produira, la post-probabilité est $ p (\ theta | x) $. .. C'est ce qu'on appelle une mise à jour bayésienne.

La distribution postérieure $ p (\ theta | x) $, qui a été mise à jour avec des informations, peut être mise à jour l'une après l'autre comme la distribution précédente $ p (\ theta) $. La réponse ne change pas même si l'ordre de multiplication est changé, donc dans l'exemple du lancer de pièces, $ x = $ {avant, avant, arrière, arrière} ou $ x = $ {avant, arrière, avant, arrière} est le dernier La distribution postérieure est la même. C'est ce qu'on appelle la rationalité séquentielle de l'estimation bayésienne. Puisque les résultats des données jusqu'à présent sont reflétés dans la distribution antérieure, on peut dire que c'est une sorte de fonction d'apprentissage [^ 2].

programme

0. Module

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Times New Roman"      #Définir la police entière
plt.rcParams["xtick.direction"] = "in"               #Vers l'intérieur de la ligne d'échelle de l'axe x
plt.rcParams["ytick.direction"] = "in"               #Ligne d'échelle de l'axe Y vers l'intérieur
plt.rcParams["xtick.minor.visible"] = True           #Ajout de l'échelle auxiliaire sur l'axe x
plt.rcParams["ytick.minor.visible"] = True           #Ajout de l'échelle auxiliaire de l'axe y
plt.rcParams["xtick.major.width"] = 1.5              #Largeur de ligne de la ligne d'échelle principale de l'axe x
plt.rcParams["ytick.major.width"] = 1.5              #Largeur de ligne de la ligne d'échelle principale de l'axe y
plt.rcParams["xtick.minor.width"] = 1.0              #Largeur de ligne de la ligne d'échelle auxiliaire de l'axe x
plt.rcParams["ytick.minor.width"] = 1.0              #Largeur de ligne de la ligne d'échelle auxiliaire de l'axe y
plt.rcParams["xtick.major.size"] = 10                #longueur de la ligne d'échelle principale sur l'axe des x
plt.rcParams["ytick.major.size"] = 10                #Longueur de la ligne d'échelle principale de l'axe y
plt.rcParams["xtick.minor.size"] = 5                 #longueur de la ligne d'échelle auxiliaire sur l'axe des x
plt.rcParams["ytick.minor.size"] = 5                 #longueur de la ligne d'échelle auxiliaire sur l'axe y
plt.rcParams["font.size"] = 14                       #Taille de police
plt.rcParams["axes.linewidth"] = 1.5                 #Épaisseur du boîtier

1. Distribution antérieure uniforme

Si vous ne savez rien à l'avance, supposez une distribution préalable uniforme pour le moment. On dit que cela s'appelle le principe de la raison insuffisante. Tout d'abord, suivons ce principe.

N = 1000 #Nombre de divisions de calcul

p_x_front = np.linspace(0.0,1.0,N) #p(x|table)
p_x_back = 1.0 - p_x_front #p(x|retour)

prior_dist = np.ones(N) #p(θ)Distribution antérieure

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,prior_dist,label="Prior distribution")
axes.scatter(np.sum(prior_dist*p_x_front/N),0.5,s=100,label="expected value")
axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

L'axe vertical est la densité de probabilité. Une fois intégré, il devient 1. Le graphique rond est la valeur attendue. La valeur sur l'axe vertical de la valeur attendue est simplement faite pour être facile à voir. output_8_1.png

def update(p_x_theta,p_theta):
    return (p_x_theta * p_theta) / np.sum(p_x_theta * p_theta / N)

fig,axes = plt.subplots(figsize=(8,6))
dist_front1_back0 = update(p_x_front,prior_dist)
axes.plot(p_x_front,dist_front1_back0,label="front=1, back=0")
axes.scatter(np.sum(dist_front1_back0*p_x_front/N),0.5,s=100)

dist_front1_back1 = update(p_x_back,dist_front1_back0)
axes.plot(p_x_front,dist_front1_back1,label="front=1, back=1")
axes.scatter(np.sum(dist_front1_back1*p_x_front/N),0.5,s=100)

dist_front2_back1 = update(p_x_front,dist_front1_back1)
axes.plot(p_x_front,dist_front2_back1,label="front=2, back=1")
axes.scatter(np.sum(dist_front2_back1*p_x_front/N),0.5,s=100)

dist_front2_back2 = update(p_x_back,dist_front2_back1)
axes.plot(p_x_front,dist_front2_back2,label="front=2, back=2")
axes.scatter(np.sum(dist_front2_back2*p_x_front/N),0.2,s=100)

axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

La fonction de mise à jour est le théorème bayésien lui-même. La distribution postérieure change chaque fois que vous obtenez des informations indiquant si une pièce est recto ou verso. Vous pouvez également voir en comparant les graphiques orange et rouge que plus vous avez de données, plus elles sont précises. output_9_1.png

Regardons la distribution postérieure lorsque le front est 100 fois et le dos 100 fois.

def updates(front,back):
    p_theta = prior_dist[:]
    
    for i in range(front):
        p_theta = update(p_x_front,p_theta)
        
    for i in range(back):
        p_theta = update(p_x_back,p_theta)
        
    return p_theta

p_theta = updates(100,100)
fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,p_theta)
axes.scatter(np.sum(p_theta*p_x_front/N),0.5,s=100)
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,12.0)
axes.set_xticks(np.linspace(0.0,1.0,11))

La dispersion est devenue beaucoup plus petite. La valeur attendue est de 0,5. output_12_1.png

2. Distribution antérieure biaisée

L'estimation bayésienne peut refléter la subjectivité. Cette fois, je pense que la table sera facile à sortir, donc j'utiliserai la distribution beta de $ \ alpha = 3, \ beta = 2 $ pour la distribution précédente.

a = 3.0
b = 2.0
x = np.linspace(0.0,1.0,N)
prior_dist = x**(a-1.0) * (1.0-x)**(b-1.0)
prior_dist /= np.sum(prior_dist) / N

fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,prior_dist,label="Prior distribution")
axes.scatter(np.sum(prior_dist*p_x_front/N),0.5,s=100,label="expected value")
axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

La valeur attendue est d'environ 0,6. output_15_1.png

fig,axes = plt.subplots(figsize=(8,6))
dist_front1_back0 = update(p_x_front,prior_dist)
axes.plot(p_x_front,dist_front1_back0,label="front=1, back=0")
axes.scatter(np.sum(dist_front1_back0*p_x_front/N),0.5,s=100)

dist_front1_back1 = update(p_x_back,dist_front1_back0)
axes.plot(p_x_front,dist_front1_back1,label="front=1, back=1")
axes.scatter(np.sum(dist_front1_back1*p_x_front/N),0.5,s=100)

dist_front2_back1 = update(p_x_front,dist_front1_back1)
axes.plot(p_x_front,dist_front2_back1,label="front=2, back=1")
axes.scatter(np.sum(dist_front2_back1*p_x_front/N),0.5,s=100)

dist_front2_back2 = update(p_x_back,dist_front2_back1)
axes.plot(p_x_front,dist_front2_back2,label="front=2, back=2")
axes.scatter(np.sum(dist_front2_back2*p_x_front/N),0.2,s=100)

axes.legend(loc="best")
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,2.5)
axes.set_xticks(np.linspace(0.0,1.0,11))

Par rapport au cas de l'utilisation d'une distribution antérieure uniforme, la valeur attendue est plus élevée en raison de la subjectivité que le tableau est susceptible d'apparaître. output_16_1.png

Regardons la distribution postérieure lorsque l'avant et l'arrière sont à nouveau 100 fois.

p_theta = updates(100,100)
fig,axes = plt.subplots(figsize=(8,6))
axes.plot(p_x_front,p_theta)
axes.scatter(np.sum(p_theta*p_x_front/N),0.5,s=100)
axes.set_xlabel(r"$p(\theta)$")
axes.set_ylabel("Probability density")
axes.set_ylim(-0.1,12.0)
axes.set_xticks(np.linspace(0.0,1.0,11))

La valeur attendue était de 0,502, ce qui était légèrement supérieur à 0,5, mais à mesure que le nombre de données augmentait, l'influence de la subjectivité initiale diminuait progressivement et une distribution précise pouvait être obtenue.

output_17_1.png

Résumé

J'ai eu une idée approximative de l'estimation bayésienne. Je ne comprends pas la méthode MCMC, alors voici [Introduction à l'analyse des données avec la modélisation statistique bayésienne commençant par R et Stan](https://logics-of-blue.com/r-stan-bayesian-model-intro-book J'espère pouvoir lire -support /) et écrire sur un contenu plus avancé cette fois.

[^ 1]: Hiroyuki Kojima "Introduction aux statistiques complètes de Bayes d'auto-apprentissage" Diamond, novembre 2015, p.6 [^ 2]: Hiroyuki Kojima "Introduction aux statistiques complètes de Bayes d'auto-apprentissage" Diamond, novembre 2015, p.145

Recommended Posts

Estimation bayésienne en utilisant le lancer de pièces comme exemple Voir la distribution postérieure pour chaque essai
Comprendre la fonction de convolution en utilisant le traitement d'image comme exemple