はじめに
ベイズ更新は、新しいデータが利用可能になるときに仮説の確率を更新することができる統計的アプローチです。この実験では、Matplotlib を使ってベイズ更新がどのように機能するかを示すアニメーションを作成します。具体的には、コイントス実験をシミュレートし、ベイズ更新を使ってコインが表を向いて止まる確率を推定します。
VM のヒント
VM の起動が完了したら、左上隅をクリックしてノートブックタブに切り替え、Jupyter Notebook を使って練習しましょう。
時々、Jupyter Notebook が読み込み終わるまで数秒待つ必要があります。Jupyter Notebook の制限により、操作の検証は自動化できません。
学習中に問題に遭遇した場合は、Labby にお問い合わせください。セッション後にフィードバックを提供してください。すぐに問題を解決いたします。
必要なライブラリをインポートする
まず、この実験で使用するライブラリをインポートします。具体的には、可視化に matplotlib.pyplot を、数値計算に numpy を、数学関数に math を使用します。
import math
import matplotlib.pyplot as plt
import numpy as np
ベータ分布の確率密度関数を定義する
ベータ分布は、確率の分布を表すために頻繁に使用される連続確率分布です。ベイズ更新では、ベータ分布を事前分布として使用して、データを観察する前の仮説の確率に関する信念を表現します。その後、新しいデータを観察するときにベータ分布を更新します。
ベイズ更新をシミュレートするには、ベータ分布の確率密度関数 (PDF) を計算する関数を定義する必要があります。ベータ分布の PDF で使用されるガンマ関数を計算するために、math.gamma 関数を使用できます。
def beta_pdf(x, a, b):
return (x**(a-1) * (1-x)**(b-1) * math.gamma(a + b)
/ (math.gamma(a) * math.gamma(b)))
UpdateDist クラスを定義する
次に、新しいデータが観察されるたびにベータ分布を更新するために使用される UpdateDist という名前のクラスを定義します。UpdateDist クラスは 2 つの引数を取ります：Matplotlib の軸オブジェクトと成功の初期確率。
class UpdateDist:
def __init__(self, ax, prob=0.5):
self.success = 0
self.prob = prob
self.line, = ax.plot([], [], 'k-')
self.x = np.linspace(0, 1, 200)
self.ax = ax
## Set up plot parameters
self.ax.set_xlim(0, 1)
self.ax.set_ylim(0, 10)
self.ax.grid(True)
## This vertical line represents the theoretical value, to
## which the plotted distribution should converge.
self.ax.axvline(prob, linestyle='--', color='black')
__init__ メソッドは、成功の初期数をゼロに設定 (self.success = 0)、成功の初期確率を引数として渡された値に設定 (self.prob = prob) することで、クラスインスタンスを初期化します。また、ベータ分布を表す線オブジェクトを作成し、グラフのパラメータを設定します。
__call__ メソッドは、アニメーションが更新されるたびに呼び出されます。これはコイントス実験をシミュレートし、それに応じてベータ分布を更新します。
def __call__(self, i):
## This way the plot can continuously run and we just keep
## watching new realizations of the process
if i == 0:
self.success = 0
self.line.set_data([], [])
return self.line,
## Choose success based on exceed a threshold with a uniform pick
if np.random.rand() < self.prob:
self.success += 1
y = beta_pdf(self.x, self.success + 1, (i - self.success) + 1)
self.line.set_data(self.x, y)
return self.line,
アニメーションの最初のフレームの場合 (if i == 0)、成功の数をゼロにリセットし、線オブジェクトをクリアします。それ以外の場合、0 から 1 の間の乱数 (np.random.rand()) を生成し、それを成功の確率 (self.prob) と比較することでコイントス実験をシミュレートします。乱数が成功の確率より小さい場合、それを成功として数え、beta_pdf 関数を使ってベータ分布を更新します。最後に、新しいデータで線オブジェクトを更新して返します。
アニメーションを作成する
これで UpdateDist クラスを定義したので、Matplotlib の FuncAnimation クラスを使ってアニメーションを作成できます。グラフオブジェクトと軸オブジェクトを作成し、軸オブジェクトを UpdateDist クラスに渡して、そのクラスの新しいインスタンスを作成します。
fig, ax = plt.subplots()
ud = UpdateDist(ax, prob=0.7)
anim = FuncAnimation(fig, ud, frames=100, interval=100, blit=True)
plt.show()
FuncAnimation クラスはいくつかの引数を取ります。
fig：グラフオブジェクトud：UpdateDistインスタンスframes：アニメーションするフレーム数interval：フレーム間の時間（ミリ秒）blit：変更された部分のみを更新するかどうか
結果を解釈する
アニメーションは、新しいデータが観察されるにつれてベータ分布がどのように更新されるかを示しています。黒い破線は成功の真の確率（すなわち、コインが表を向く確率）を表しています。アニメーションが進行するにつれて、ベータ分布は成功の事前確率（0.7）でピークを持ち始め、より多くのデータが観察されるにつれて徐々に成功の真の確率に向かってシフトすることがわかります。
まとめ
この実験では、Matplotlib を使ってベイズ更新を示すアニメーションを作成しました。ベータ分布の確率密度関数 (PDF) を計算する関数と、新しいデータが観察されるたびにベータ分布を更新するクラスを定義しました。その後、Matplotlib の FuncAnimation クラスを使ってアニメーションを作成し、結果を解釈しました。