[PYTHON] Bayesianische Schätzung am Beispiel des Münzwurfs Siehe die hintere Verteilung für jeden Versuch

Einführung

Ich habe gehört, dass die Bayes'schen Statistiken in letzter Zeit heiß sind, deshalb dachte ich, ich sollte ["Einführung in die Bayes'schen Statistiken"] studieren (https://www.diamond.co.jp/book/9784478013328.html). Es ist auf halbem Weg leicht zu verstehen, aber ich fand es schwierig zu verstehen, sobald die kontinuierliche Verteilung herauskam, weil es ein Stil ist, der mit Zahlen erklärt, ohne so viele mathematische Formeln wie möglich zu verwenden. Zum Lernen werde ich das Münzwurf als einfaches Beispiel verwenden und bei jedem Werfen die hintere Verteilung betrachten. Ich verstehe die MCMC-Methode noch nicht, daher werde ich die Wahrscheinlichkeitsdichtefunktion in kleinen Stücken berechnen.

Referenz

Hiroyuki Kojima "Einführung in die Bayes'sche Statistik", Diamond, November 2015

Die Erklärung basiert auf der gleichen Abbildung wie im obigen Buch am Beispiel des Münzwurfs.

Dies ist das gleiche Thema, aber es handelt sich um einen erweiterten Inhalt, der mit der MCMC-Methode unter Verwendung von pyMC berechnet wurde. Ich möchte auch pyMC verwenden können. Sehen Sie sich die hintere Verteilung der Vorder- und Rückseite der Münze für jeden Versuch an (PyMC)

Dies ist das erste Kapitel der "Einführung in die Bayes'sche Statistik". Die Figur wird oft verwendet und ist leicht zu verstehen. Bayes Schätzung, um in 5 Minuten klar zu verstehen

Bayesianische Schätzung

Die Stärke der Bayes'schen Statistik ist

  1. Die Eigenschaft, dass es erraten werden kann, selbst wenn nur wenige Daten vorhanden sind, und die genauer wird, wenn mehr Daten vorhanden sind.
  2. Lernfunktion, die Vermutungen als Reaktion auf eingehende Informationen sofort automatisch aktualisiert Es scheint in [^ 1] zu sein.

Die Bayes'sche Schätzung basiert auf dem Bayes'schen Theorem. $p(\theta|x)=p(\theta) \frac{p(x|\theta)}{\int p(x|\theta)p(\theta)d\theta}$ Hier\thetaIst ein Parameter,p(\theta)Ist vorherige Verteilung,p(x|\theta)Ist\thetaWannxBedingte Wahrscheinlichkeit von(Haftung)、p(\theta|x)Ist事後分布です。 Im Beispiel des Münzwurfs ist $ \ theta $ die Vorder- oder Rückseite, $ p (\ theta) $ die Wahrscheinlichkeitsverteilung von $ \ theta $, $ x $ die Daten darüber, ob die Münze vorne oder hinten ist und der Nenner ist Stellt die Wahrscheinlichkeit dar, $ x $ zu erhalten. Die ursprüngliche Wahrscheinlichkeit des Auftretens von $ \ theta $ war die Vorwahrscheinlichkeit $ p (\ theta) $, aber mit der Information, dass $ x $ auftreten wird, ist die Nachwahrscheinlichkeit $ p (\ theta | x) $. .. Dies wird als Bayes'sches Update bezeichnet.

Die hintere Verteilung $ p (\ theta | x) $, die mit Informationen aktualisiert wurde, kann erneut als vorherige Verteilung $ p (\ theta) $ aktualisiert werden. Die Antwort ändert sich auch dann nicht, wenn die Reihenfolge der Multiplikation geändert wird. Im Beispiel des Münzwurfs ist also $ x = $ {vorne, vorne, hinten, hinten} oder $ x = $ {vorne, hinten, vorne, hinten} das Finale Die hintere Verteilung ist die gleiche. Dies wird als sequentielle Rationalität der Bayes'schen Schätzung bezeichnet. Da sich die bisherigen Ergebnisse der Daten in der vorherigen Verteilung widerspiegeln, kann von einer Art Lernfunktion gesprochen werden [^ 2].

Programm

0. Modul

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Times New Roman"      #Stellen Sie die gesamte Schriftart ein
plt.rcParams["xtick.direction"] = "in"               #Nach innen die x-Achsen-Skalierungslinie
plt.rcParams["ytick.direction"] = "in"               #Y-Achsen-Skalierungslinie nach innen
plt.rcParams["xtick.minor.visible"] = True           #Hinzufügen einer x-Achsen-Hilfsskala
plt.rcParams["ytick.minor.visible"] = True           #Hinzufügen einer Hilfsskala der y-Achse
plt.rcParams["xtick.major.width"] = 1.5              #Linienbreite der x-Achsen-Hauptskalenlinie
plt.rcParams["ytick.major.width"] = 1.5              #Linienbreite der Hauptskalenlinie der y-Achse
plt.rcParams["xtick.minor.width"] = 1.0              #Linienbreite der Hilfsskalenlinie der x-Achse
plt.rcParams["ytick.minor.width"] = 1.0              #Linienbreite der Hilfsskalenlinie der y-Achse
plt.rcParams["xtick.major.size"] = 10                #Länge der x-Achsen-Hauptskalenlinie
plt.rcParams["ytick.major.size"] = 10                #Länge der Hauptmaßstablinie der y-Achse
plt.rcParams["xtick.minor.size"] = 5                 #Länge der Hilfsskala der x-Achse
plt.rcParams["ytick.minor.size"] = 5                 #Länge der Hilfsskala der y-Achse
plt.rcParams["font.size"] = 14                       #Schriftgröße
plt.rcParams["axes.linewidth"] = 1.5                 #Gehäusedicke

1. Gleichmäßige vorherige Verteilung

Wenn Sie im Voraus nichts wissen, gehen Sie vorerst von einer einheitlichen Vorverteilung aus. Es wird gesagt, dass dies [das Prinzip der unzureichenden Vernunft] genannt wird (https://clover.fcg.world/2016/07/30/6249/). Folgen wir zunächst diesem Prinzip.

N = 1000 #Anzahl der Berechnungsabteilungen

p_x_front = np.linspace(0.0,1.0,N) #p(x|Tabelle)
p_x_back = 1.0 - p_x_front #p(x|zurück)

prior_dist = np.ones(N) #p(θ)Vorherige Verteilung

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))

Die vertikale Achse ist die Wahrscheinlichkeitsdichte. Wenn es integriert ist, wird es 1. Das runde Diagramm ist der erwartete Wert. Der Wert auf der vertikalen Achse des erwarteten Werts ist nur leicht erkennbar. 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))

Die Aktualisierungsfunktion ist das Bayes'sche Theorem selbst. Die hintere Verteilung ändert sich jedes Mal, wenn Sie Informationen darüber erhalten, ob eine Münze vorne oder hinten ist. Sie können auch durch Vergleichen der orangefarbenen und roten Diagramme feststellen, dass die Daten umso genauer sind, je mehr Daten Sie haben. output_9_1.png

Schauen wir uns die hintere Verteilung an, wenn die Vorderseite 100-mal und die Rückseite 100-mal ist.

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))

Die Dispersion ist viel kleiner geworden. Der erwartete Wert beträgt 0,5. output_12_1.png

2. Voreingenommene vorherige Verteilung

Die Bayes'sche Schätzung kann die Subjektivität widerspiegeln. Dieses Mal denke ich, dass die Tabelle leicht herauszukommen sein wird, also werde ich die Beta-Verteilung von $ \ alpha = 3, \ beta = 2 $ für die vorherige Verteilung verwenden.

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))

Der erwartete Wert liegt bei ca. 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))

Im Vergleich zu einer einheitlichen vorherigen Verteilung ist der erwartete Wert aufgrund der Subjektivität, mit der die Tabelle wahrscheinlich erscheint, höher. output_16_1.png

Schauen wir uns die hintere Verteilung an, wenn Vorder- und Rückseite wieder 100-mal sind.

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))

Der erwartete Wert betrug 0,502, was etwas größer als 0,5 war, aber mit zunehmender Anzahl von Daten nahm der Einfluss der anfänglichen Subjektivität allmählich ab und es konnte eine genaue Verteilung erhalten werden.

output_17_1.png

Zusammenfassung

Ich habe eine ungefähre Vorstellung von der Bayes'schen Schätzung. Ich verstehe die MCMC-Methode nicht, daher folgt als Nächstes [Einführung in die Datenanalyse mit Bayesian Statistical Modeling, beginnend mit R und Stan](https://logics-of-blue.com/r-stan-bayesian-model-intro-book Ich hoffe, ich kann diesmal -support /) lesen und über fortgeschrittenere Inhalte schreiben.

[^ 1]: Hiroyuki Kojima "Einführung in die vollständige Bayes-Statistik zum Selbststudium" Diamond, November 2015, S. 6 [^ 2]: Hiroyuki Kojima "Einführung in die vollständige Bayes-Statistik zum Selbststudium" Diamond, November 2015, S.145

Recommended Posts

Bayesianische Schätzung am Beispiel des Münzwurfs Siehe die hintere Verteilung für jeden Versuch
Verstehen Sie die Funktion der Faltung am Beispiel der Bildverarbeitung