Pythonで体験するベイズ推論 PyMCによるMCMC入門 第6章（6.4まで）写経

目次

始めに

　この記事では、Pythonで体験するベイズ推論 PyMCによるMCMC入門

の第6章の6.4までを写経した記録をまとめます。前回はこちらPythonで体験するベイズ推論 PyMCによるMCMC入門 第5章 写経。 実行はJupyter Labにて実行しています。

　写経は以下のようにまとめました。

• 写経ではありますが、関数を作成し、コードの再実行をやりやすいようにした部分があります。

• より良いと思われるコードを考えた場合は書き換えてコメントを添えるようにし、変更点をなるべく明示するようにしてあります。

• 個人的に気になった点のメモを残すようにしてあります。同じような疑問を持った方の助けになれば幸いです。

• 以前書いたコードと同じようなコード（例えばグラフの描写等）は効率化のために飛ばしているところもあります。

• 記事内で使用するモジュールなどは一番最初に宣言するようにしてあります。

　Pythonで体験するベイズ推論 PyMCによるMCMC入門を読んでいている際のちょっとした正誤表代わりになればと思います。

この記事で使用する主なモジュール、設定

　この記事では、主に以下のモジュールや設定を使用しています。

from pymc import rbeta
import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
from urllib.request import urlretrieve

6 事前分布をハッキリさせよう

6.4 例題 : ベイズ多腕バンディット

class Bandits(object):
    def __init__(self, p_array):
        self.p = p_array
        self.optimal = np.argmax(p_array)
        
    def pull(self, i):
        return np.random.rand() < self.p[i]
    
    def __len__(self):
        return len(self.p)
    
class BayesianStrategy(object):
    def __init__(self, bandits):
        self.bandits = bandits
        n_bandits = len(self.bandits)
        self.wins = np.zeros(n_bandits)
        self.trials = np.zeros(n_bandits)
        self.N = 0
        self.choices = []
        self.bb_score = []
        
    def sample_bandits(self, n=1):
        bb_score = np.zeros(n)
        choices = np.zeros(n)
        
        for k in range(n):
            choice = np.argmax(rbeta(1 + self.wins, 1 + self.trials - self.wins))
            
            result = self.bandits.pull(choice)
            
            self.wins[choice] += result
            self.trials[choice] += 1
            bb_score[k] = result
            self.N += 1
            choices[k] = choice
            
        self.bb_score = np.r_[self.bb_score, bb_score]
        self.choices = np.r_[self.choices, choices]
        return

メモ（__len__メソッド）

　Banditsクラス定義の中の__len__は特殊メソッドと呼ばれるもののひとつです。特殊メソッドはPythonであらかじめ定義されています。

参考：3. データモデル、object.__len__(self)

特殊メソッド__len__はlen(self.bandits)のようにメソッドを呼び出すことができ、ドット（.）による通常のメソッドの呼び出しとは異なります。むしろ、これが目的と言っていいのではないでしょうか？（間違っていたら指摘してください。）

メモ（np.argmax(rbeta(1 + self.wins, 1 + self.trials - self.wins))）

　np.argmax(rbeta(1 + self.wins, 1 + self.trials - self.wins))の補足をしておきます。

　まず、rbeta関数はベータ分布に従う確率変数です：

rbeta(1, 1)

Output:

    0.52678815388438394

　さらに、np.argmax(rbeta(1 + self.wins, 1 + self.trials - self.wins))のように、同じサイズのarrayを引数に渡すと、各要素をパラメータに持つベータ分布に従う確率変数となります：

beta_dist = rbeta(np.array([1, 1, 1]), np.array([1, 2, 3]))
beta_dist

Output:

    array([ 0.93062631,  0.73661115,  0.5891043 ])

　np.argmax(rbeta(1 + self.wins, 1 + self.trials - self.wins))は上のようなarrayの要素の中で最も値の大きい（もっとも当たりやすい）要素のインデックスを返します：

np.argmax(beta_dist)

Output:

    0

figsize(22, 20)
plt.rcParams['font.size'] = 17

beta = stats.beta
x = np.linspace(0.001, .999, 200)

def plot_priors(bayesian_strategy, prob, lw=3, alpha=0.2, plt_vlines=True):
    wins = bayesian_strategy.wins
    trials = bayesian_strategy.trials
    for i in range(prob.shape[0]):
        y = beta(1 + wins[i], 1 + trials[i] - wins[i])
        p = plt.plot(x, y.pdf(x), lw=lw)
        c = p[0].get_markeredgecolor()
        plt.fill_between(x, y.pdf(x), 0, color=c, alpha=alpha,
                         label=f"underlying probability: {prob[i]}")
        if plt_vlines:
            plt.vlines(prob[i], 0, y.pdf(prob[i]),
                       colors=c, linestyles="--", lw=2)
        plt.autoscale(tight="True")
        plt.title(f"Posteriors After {bayesian_strategy.N} pull"\
                  + "s" * (bayesian_strategy.N > 1))
        plt.autoscale(tight=True)
    return

hidden_prob = np.array([0.85, 0.60, 0.75])
bandits = Bandits(hidden_prob)
bayesian_strat = BayesianStrategy(bandits)

draw_samples = [1, 1, 3, 10, 10, 25, 50, 100, 200, 600]

for j, i in enumerate(draw_samples):
    plt.subplot(5, 2, j + 1)
    bayesian_strat.sample_bandits(i)
    plot_priors(bayesian_strat, hidden_prob)
    plt.autoscale(tight=True)
plt.tight_layout()

6.4.3 良さを測る

urlretrieve("https://git.io/vXL9A", "other_strats.py")

Output:

    ('other_strats.py', <http.client.HTTPMessage at 0x1cc24eb0ac8>)

figsize(25, 12)

from other_strats import *

hidden_prob = np.array([0.15, 0.2, 0.1, 0.05])
bandits = Bandits(hidden_prob)

def regret(probabilities, choices):
    w_opt = probabilities.max()
    return (w_opt - probabilities[choices.astype(int)]).cumsum()

strategies = [upper_credible_choice,
              bayesian_bandit_choice,
              ucb_bayes,
              max_mean,
              random_choice]
algos = []
for strat in strategies:
    algos.append(GeneralBanditStrat(bandits, strat))

for strat in algos:
    strat.sample_bandits(10000)

for i, strat in enumerate(algos):
    _regret = regret(hidden_prob, strat.choices)
    plt.plot(_regret, label=strategies[i].__name__, lw=3)

plt.title("Total Regret of Bayesian Bandits Strategy vs. Random guessing")
plt.xlabel("Number of pulls")
plt.ylabel("Regret after $n$ pulls");
plt.legend(loc="upper left");

メモ（cumsumメソッド）

　arrayのcumsumは、与えられたarrayの要素を足し合わせたarrayを返します

参考：numpy.cumsum

例えば、以下のようになります：

np.array([i for i in range(11)]).cumsum()

Output:

    array([ 0,  1,  3,  6, 10, 15, 21, 28, 36, 45, 55], dtype=int32)

regret関数の(w_opt - probabilities[choices.astype(int)]).cumsum()では、w_optと実際の行動であるchoicesの差を足し合わせてregretを計算しています。

trials = 200
expected_total_regret = np.zeros((1000, 3))

for i_strat, strat in enumerate(strategies[:-2]):
    for i in range(trials):
        general_strat = GeneralBanditStrat(bandits, strat)
        general_strat.sample_bandits(1000)
        _regret = regret(hidden_prob, general_strat.choices)
        expected_total_regret[:, i_strat] += _regret
    plt.plot(expected_total_regret[:, i_strat] / trials, lw=3, label=strat.__name__)

plt.title("Expected Total Regret of Multi-armed Bandit strategies")
plt.xlabel("Number of pulls")
plt.ylabel("Exepected Total Regret \n after $n$ pulls");
plt.legend(loc="upper left");

plt.figure()
[pl1, pl2, pl3] = plt.plot(expected_total_regret[:, [0, 1, 2]], lw=3)
plt.xscale("log")
plt.legend([pl1, pl2, pl3],
           ["Upper Credible Bound", "Bayesian Bandit", "UCB-Bayes"],
           loc="upper left")
plt.ylabel("Exepected Total Regret \n after $\log{n}$ pulls");
plt.title("log-scale of above");
plt.ylabel("Exepected Total Regret \n after $\log{n}$ pulls");

figsize(24, 16)
beta = stats.beta
hidden_prob = beta.rvs(1, 13, size=35)
print(hidden_prob)
bandits = Bandits(hidden_prob)
bayesian_strat = BayesianStrategy(bandits)

for j, i in enumerate([100, 200, 500, 1300]):
    plt.subplot(2, 2, j + 1)
    bayesian_strat.sample_bandits(i)
    plot_priors(bayesian_strat, hidden_prob, lw=2, alpha=0.0, plt_vlines=False)
    # plt.legend()
    plt.xlim(0, 0.5)

Output:

    [ 0.03470568  0.0838137   0.06460904  0.00633369  0.01256643  0.07361407
      0.00803128  0.0054401   0.0757274   0.06154705  0.05111673  0.17662841
      0.13431165  0.05624876  0.01879327  0.07184054  0.10648022  0.03216907
      0.10596969  0.10426467  0.07715617  0.00247444  0.14421856  0.10918469
      0.08518456  0.05307192  0.02775096  0.13794162  0.0101437   0.04057421
      0.0352338   0.06571338  0.12687146  0.05107877  0.01482898]
    

　今回の写経は以上です。 ここまで読んでいただきありがとうございました。