第10章:ベイジアンML-ダイナミックSR 第1節: 共役事前分布の更新

はじめに

皆さんが投資をしている時は、株価が上がったり下がったりして今後の株価が上がるか下がるかどうかの分布について気になることもあるのではないでしょうか？

本記事で紹介することは、そういった刻一刻で変わる環境に対応するための技術です。

今回やること

・コイントスの例から、試行回数を増やしながら分布を更新する過程について話す。
・株価が上がるか下がるかのバイナリーに直してから、同様に分布を更新させていく。

共役事前分布とは？

こちらのページが分かりやすかったので引用します。

共役事前分布とは、ベイズ統計を扱う際に、複雑な計算を回避するために考えられた事前分布です。共役事前分布に尤度をかけて事後分布を求めると、その関数形が同じ分布になります。つまり、共役事前分布を使えば、事前分布と事後分布は同じ形になります。

インポートと設定

import warnings
warnings.filterwarnings('ignore')

%matplotlib inline
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
import scipy.stats as stats
from matplotlib.ticker import FuncFormatter

※ここで注意点なのですが、本docker環境ではlatexがインストールされていないため、以下のコードは一部コメントアウトさせていただきました。

import matplotlib as mpl
# mpl.rcParams['text.usetex'] = True
# mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']

np.random.seed(42)
sns.set_style('dark')

フォーマット補助関数

def format_plot(axes, i, p, y, trials, success, true_p, tmle, tmap=None):
   fmt = FuncFormatter(lambda x, _: f'{x:.0%}')
   if i >= 6:
       axes[i].set_xlabel("$p$, Success Probability")
       axes[i].xaxis.set_major_formatter(fmt)
   else:
       axes[i].axes.get_xaxis().set_visible(False)
   if i % 3 == 0:
       axes[i].set_ylabel("Posterior Probability")
   axes[i].set_yticks([], [])

   axes[i].plot(p, y, lw=1, c='k')
   axes[i].fill_between(p, y, color='darkblue', alpha=0.4)
   axes[i].vlines(true_p, 0, max(10, np.max(y)), color='k', linestyle='--', lw=1)
   axes[i].set_title(f'Trials: {trials:,d} - Success: {success:,d}')
   if i > 0:
       smle = r"$\theta_{{\mathrm{{MLE}}}}$ = {:.2%}".format(tmle)
       axes[i].text(x=.02, y=.85, s=smle, transform=axes[i].transAxes)
       smap = r"$\theta_{{\mathrm{{MAP}}}}$ = {:.2%}".format(tmap)
       axes[i].text(x=.02, y=.75, s=smap, transform=axes[i].transAxes)    
   return axes[i]

コイントスのシミュレーションと事後分布の更新(Updates of Posterior)

n_trials = [0, 1, 3, 5, 10, 25, 50, 100, 500]
outcomes = stats.bernoulli.rvs(p=0.5, size=n_trials[-1])
p = np.linspace(0, 1, 100)
# uniform (uninformative) prior
a = b = 1

fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(14, 7), sharex=True)
axes = axes.flatten()
fmt = FuncFormatter(lambda x, _: f'{x:.0%}')
for i, trials in enumerate(n_trials):
   successes = outcomes[:trials]
   theta_mle = np.mean(successes)
   heads = sum(successes)
   tails = trials - heads
   update = stats.beta.pdf(p, a + heads , b + tails)
   theta_map = pd.Series(update, index=p).idxmax()
   axes[i] = format_plot(axes, i, p, update, trials=trials, success=heads, 
                         true_p=.5, tmle=theta_mle, tmap=theta_map)

title = 'Bayesian Probabilities: Updating the Posterior'
fig.suptitle(title,  y=1.02, fontsize=14)
fig.tight_layout()

試行回数を増やせば増やすほど事後分布が形を形成しているのが分かりますね(*'ω'*)

株価の動きと事後分布の更新

バイナリー化されたS＆P 500の日次リターンのさまざまなサイズのサンプルを収集します。間隔[0, 1]の各成功確率に等しい確率を割り当てる情報のない事前分布から始めて、異なる証拠サンプルの事後確率を計算します。

sp500_returns = pd.read_hdf('../data/assets.h5', key='sp500/fred').loc['2010':, 'close']
sp500_binary = (sp500_returns.pct_change().dropna() > 0).astype(int)

sp500_binary.sum()/sp500_binary.count()
'''
0.5291411042944786
'''

次のコードサンプルは、更新が成功と失敗の観測された数を事前分布のパラメーターに単純に追加して事後分布を取得することからなることを示しています。

結果の事後分布を以下にプロットします。一様分布から次第にピークに到達するような分布へと進化しているのが分かります。 500個のサンプルの後、確率は2010年から2017年までの正の変動の実際の確率52.91％に集中しています。また、一様分布からはMLEとMAP(最大事後確率)の推定値のわずかにズレています。

n_days = [0, 1, 3, 5, 10, 25, 50, 100, 500]
# random sample of trading days
# outcomes = sp500_binary.sample(n_days[-1])

# initial 500 trading days
outcomes = sp500_binary.iloc[:n_days[-1]]
p = np.linspace(0, 1, 100)

# uniform (uninformative) prior
a = b = 1

fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(14, 7), sharex=True)
axes = axes.flatten()
for i, days in enumerate(n_days):
   successes = outcomes.iloc[:days]
   theta_mle = successes.mean()
   up = successes.sum()
   down = days - up
   update = stats.beta.pdf(p, a + up , b + down)
   theta_map = pd.Series(update, index=p).idxmax()
   axes[i] = format_plot(axes, i, p, update, trials=days, success=up, 
                         true_p=sp500_binary.mean(), tmle=theta_mle, tmap=theta_map)

title = 'Bayesian Probabilities: Updating the Posterior'
fig.suptitle(title,  y=1.02, fontsize=14)
fig.tight_layout()