LoginSignup
1
3

More than 1 year has passed since last update.

@nijigen_plotin株式会社リサーチ・アンド・イノベーション

成功確率の観測データにはベイズ更新が便利

Last updated at Posted at 2021-06-20

はじめに

成功 or 失敗のベルヌーイ試行において、観測データの成功回数 / 試行回数の算術平均の値を成功確率として求めることが多いと思いますが
ベイズ更新を用いるともっと便利ですよ という話です。

例

算術平均で成功確率を求めると、以下のようになります。

ユーザーAさん

試行回数 試行結果 成功確率(成功回数 / 試行回数)
0 -- --
1 成功 100%
2 失敗 50%
3 成功 66.66%
4 成功 75%
5 失敗 60%

ユーザーBさん

試行回数 試行結果 成功確率(成功回数 / 試行回数)
0 -- --
1 成功 100%
2 成功 100%
3 成功 100%
... (4 ~ 99はすべて成功) 100%
100 成功 100%

困る点

  • 試行回数が0の時は成功確率が入らない
    • データ持ちの関係でNULLにできない場合等
  • 試行回数が加味されていない
    • 試行回数が1回 かつ 成功回数が1回 の人と試行回数が100回 かつ成功回数が100回の人はどちらも成功確率100%
    • 沢山試行していてずっと成功している人のほうが良い という風にしたい

ベイズ更新

事前共役分布をベータ分布、ベルヌーイ試行の観測データを用いて事後分布（ベータ分布）を作成する、ベイズ更新を用いると、困る点の問題が解決します。
詳しい内容は後述するとして、先にベーシックなベイズ更新での成功確率の求め方を記載します。

成功回数 + 1 = $\alpha$
失敗回数 + 1 = $\beta$
とすると

  • 成功確率(期待値)
\frac{\alpha}{\alpha + \beta}
  • 成功確率(分散) 
\frac{\alpha \beta}{(\alpha + \beta)^{2} \cdot (\alpha + \beta + 1)}

ん？分散？となりますが、今回のベイズ更新で返ってくるのはベータ分布という連続確率分布ですので、ばらつきを示す分散の値も算出されます。
先ほどのユーザーAさん、ユーザーBさんのテーブルにベイズ更新の情報も付け足してみます。

ユーザーAさん_2

試行回数 試行結果 成功確率(成功回数 / 試行回数) 成功確率(期待値) 成功確率(分散)
0 -- -- 50% 0.0833
1 成功 100% 66.66% 0.0555
2 失敗 50% 50% 0.05
3 成功 66.66% 60% 0.04
4 成功 75% 66.66% 0.0317
5 失敗 60% 57.14% 0.0306

ユーザーBさん_2

試行回数 試行結果 成功確率(成功回数 / 試行回数) 成功確率(期待値) 成功確率(分散)
0 -- -- 50% 0.0833
1 成功 100% 66.66% 0.0555
2 成功 100% 75% 0.0375
3 成功 100% 80% 0.0266
... (4 ~ 99はすべて成功) 100% 83.33% ~ 99% 0.0198 ~ 0.00096
100 成功 100% 99.01% 0.000094

試行回数が0の時は成功確率が入らない について

試行回数が0の時も期待値が50%、分散が0.0833と値を入れることができます。
事前分布というのは任意のパラメータを設定することができるので、成功確率0 ~ 100%まで一様になる設定にしています。
「最初から期待値を70%にしておく」 といった設定も可能です。
事前分布の設定が違っていても、試行回数を重ねるごとに成功確率(期待値)は1つの値へと収束していきます。（真の確率pがある場合）

試行回数が加味されていない について

成功回数 / 試行回数 の成功確率ではユーザーAさんの1回目とユーザーBさんの100回目の値が同じです。
が、成功確率(期待値)では、ユーザーAさんが66.66%、ユーザーBさんが99.01%でBさんのほうが値が高くなっています。

更新の様子

数字だけだと分かりづらいと思うので、試行結果を受けて分布がどう変わっていくかをグラフにしてみました。

ユーザーAさん
A_sample_animation.gif

ユーザーBさん
B_sample_animation.gif

このように成功or失敗の観測データを受けて分布が形を変えて更新されていきます。
更新されるたびに、分布の山は鋭くなり、成功確率(期待値)まわりのばらつきが減っていきます。
成功回数 / 試行回数 = 成功確率の算術平均と違い、回数が加味できるのが良い所であり、計算も簡単で理解しやすいです。

データ生成→グラフ描画コード
import numpy as np
from scipy.stats import beta, bernoulli
import matplotlib.pyplot as plt
import matplotlib.animation as ArtistAnimation
plt.style.use('seaborn-darkgrid')
plt.rcParams['font.family'] = 'Yu Gothic'
plt.rcParams['font.size'] = 17

# ユーザーA,Bの観測データ
A_sample = bernoulli.rvs(p=0.6, size=100, random_state=3)
B_sample = np.ones(100)

# alpha, betaを作成
def make_ab(sample_data):
    a = 1
    b = 1
    sample_a = []
    sample_b = []
    sample_a.append(a)
    sample_b.append(b)
    for x in sample_data:
        if x == 1:
            a += 1
        elif x == 0:
            b += 1
        sample_a.append(a)
        sample_b.append(b)
    return np.array(sample_a), np.array(sample_b)

A_beta_process = make_ab(A_sample)
B_beta_process = make_ab(B_sample)
x = np.linspace(0.01, 1, 100)

# gifアニメーション描画
%matplotlib nbagg
fig, ax = plt.subplots(figsize=(5, 5))
ims = []
for i in range(len(A_beta_process[0])):
    lines = ax.plot(x, beta.pdf(x, a=A_beta_process[0][i], b=A_beta_process[1][i]), color='royalblue')
    title = ax.text(0.5, 1.01, 
                    f'試行 {i} 回目 平均={round(beta.mean(A_beta_process[0][i], A_beta_process[1][i]), 2)} 分散={round(beta.var(A_beta_process[0][i], A_beta_process[1][i]), 3)}'
                    , ha='center', va='bottom', transform=ax.transAxes)
    # fig.fill_between(x, np.zeros(len(x)), beta.pdf(x, a=A_beta_process[0][i], b=A_beta_process[1][i]), alpha=0.5)
    # plt.ylim(0, 4)
    ims.append(lines + [title])
    if i == 20:
        break
ani = animation.ArtistAnimation(fig, ims, interval=500)
ani.save('A_sample_animation.gif', writer='imagemagick')
plt.show()

ちなみに、ユーザーAさんの真の成功確率pを60%としたとき、沢山の試行回数をこなせば算術平均でもベイズ更新を用いた成功確率(期待値)でも最終的に示す値はとても近くなります。

算術平均と成功確率(期待値)の比較

ユーザーAさん
A_sample_comparison.gif

ユーザーBさん

B_sample_comparison.gif

ユーザーあたりの試行回数が多い観測データであればどちらを使ってもそこまで差は生じないと思いますが、試行回数が少数のデータが多い場合は、おすすめです。

グラフ描画コード
# 算術平均とalpha, betaを作る
A_sample_1 = bernoulli.rvs(p=0.6, size=1000, random_state=3)
A_beta_process_1 = make_ab(A_sample_1)
A_percent = np.insert(np.cumsum(A_sample_1) / np.array(range(1,len(A_sample_1) + 1)),0,np.nan)
B_sample_1 = np.ones(1000)
B_beta_process_1 = make_ab(B_sample_1)
B_percent = np.insert(np.cumsum(B_sample_1) / np.array(range(1,len(B_sample_1) + 1)),0,np.nan)


# 描画
%matplotlib nbagg
fig, ax = plt.subplots(figsize=(5, 5))
ims = []
for i in range(len(A_beta_process_1[0])):
    lines = ax.plot(x, beta.pdf(x, a=A_beta_process_1[0][i], b=A_beta_process_1[1][i]), color='royalblue')
    lines2 = ax.plot([A_percent[i], A_percent[i]], [0, 15], color='red')
    # bars = [ax.bar(A_percent[i], 4, width=0.02, color='red', alpha=0.5)]
    plt.xlim(0, 1)
    title = ax.text(0.5, 1.01, 
                    f'試行 {i} 回目 算術平均={round(A_percent[i],2)} 期待値={round(beta.mean(A_beta_process_1[0][i], A_beta_process_1[1][i]), 2)} '
                    , ha='center', va='bottom', transform=ax.transAxes)
    # fig.fill_between(x, np.zeros(len(x)), beta.pdf(x, a=A_beta_process[0][i], b=A_beta_process[1][i]), alpha=0.5)
    # plt.ylim(0, 4)
    ims.append(lines + lines2 + [title])
    # ims.append([title])
    if i == 70:
        break
ani = animation.ArtistAnimation(fig, ims, interval=250)
ani.save('A_sample_comparison.gif', writer='imagemagick')
plt.show()


fig, ax = plt.subplots(figsize=(5, 5))
ims = []
for i in range(len(B_beta_process_1[0])):
    lines = ax.plot(x, beta.pdf(x, a=B_beta_process_1[0][i], b=B_beta_process_1[1][i]), color='royalblue')
    lines2 = ax.plot([B_percent[i], B_percent[i]], [0, 70], color='red', lw=3, alpha=0.8)
    # bars = [ax.bar(A_percent[i], 4, width=0.02, color='red', alpha=0.5)]
    plt.xlim(0, 1)
    title = ax.text(0.5, 1.01, 
                    f'試行 {i} 回目 算術平均={round(B_percent[i],2)} 期待値={round(beta.mean(B_beta_process_1[0][i], B_beta_process_1[1][i]), 2)} '
                    , ha='center', va='bottom', transform=ax.transAxes)
    # fig.fill_between(x, np.zeros(len(x)), beta.pdf(x, a=A_beta_process[0][i], b=A_beta_process[1][i]), alpha=0.5)
    # plt.ylim(0, 4)
    ims.append(lines + lines2 + [title])
    # ims.append([title])
    if i == 70:
        break
ani = animation.ArtistAnimation(fig, ims, interval=250)
ani.save('B_sample_comparison.gif', writer='imagemagick')
plt.show()

ベイズ推定

ベイズの定理を用いた確率分布の推定をベイズ推定と言います。
事象$B_1,B_2,...$を互いに排反な原因とします。
事象$A$を考えこれを結果とすると、各原因からAが生じる条件付き確率は
$P(A|B_1), P(A|B_2)....$となります。
事象$A$を生じさせた各原因の可能性は、条件付き確率として

P(B_i|A) = \frac{P(A|B_i)P(B_i)}{\sum_{k=1}^{\infty} P(A | B_k)P(B_k)}

で計算できます。（ベイズの定理）

分母部分は各原因かつAが生じる確率
$P(A \cap B_1), P(A \cap B_2)....$の総和であり、
条件付き確率が

P(A|B) = \frac{P(A \cap B)}{P(B)}

で表せることから分母を払って

P(A|B)P(B) = P(A \cap B)

となっています。

それぞれの確率は
$P(B_i)$ = 事前確率
$P(B_i | A)$ = 事後確率
$P(A | B_i)$ = 尤度
と呼びます。

$\sum_{k=1}^{\infty} P(A | B_k)P(B_k)$部分は定数であるため省略すると、両辺の関係としては

P(B_i | A) \propto P(A|B_i)P(B_i)

と書くことができます。

事前確率分布・事後確率分布

今回のユーザーAさん・ユーザーBさんの例の場合、
ユーザーAさんは5回までの試行で3回成功、ユーザーBさんは5回までの試行で5回成功しています。
このように成功する確率を$p$、失敗する確率を$1-p$とし、$n$回の試行中$x$回の成功を得る確率は

f(x|p) = {_n C _x}p^x(1 - p)^{n - x}\quad x = 0,1,...,n

の二項分布$Binom(n,p)$となります。(前提として、各試行については独立である必要があります)
そして、この$f(x | p)$ですが、$P(A|B_i)$の尤度の部分となります。
この二項分布における$p$の事前確率分布を設定すれば、事後確率分布をベイズの定理から計算することができます。(両辺の関係式参照)

事前確率分布

事前確率分布として、今回の例ではベータ分布を用います。

f(x | \alpha, \beta) = \frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)}

分母はベータ関数であり、

B(\alpha, \beta) = \int_{0}^{1}x^{\alpha - 1}(1 - x)^{\beta-1}du = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)}

で表します。
ベータ分布を事前確率の$P(B_i)$とし、互いに掛ける事で事後確率$P(B_i|A)$が求まります。

事後確率分布

解釈しやすくするため、事前確率の記号を$\theta$に書き直します。

  • 二項分布($p \rightarrow \theta $)
f(x|\theta) = {_n C _x}\theta^x(1 - \theta)^{n - x}\quad x = 0,1,...,n
  • ベータ分布($x \rightarrow \theta$)
f(\theta | \alpha, \beta) = \frac{\theta^{\alpha - 1}(1 - \theta)^{\beta - 1}}{B(\alpha, \beta)}

互いに掛けると

\frac{{_nC_x \theta^{\alpha + x -1}(1 - \theta)^{\beta + (n - x) - 1}}}{B(\alpha, \beta)}

両辺の関係式を用いると

P(B_i | A) \propto \theta^{\alpha + x -1}(1 - \theta)^{\beta + (n - x) - 1}

となり、事後分布は

f(\theta | \alpha + x, \beta + (n - x)) = \frac{\theta^{\alpha + x - 1}(1 - \theta)^{\beta + (n - x) - 1}}{B(\alpha + x, \beta + (n - x))}

となります。
(互いに掛けた$\frac{_nC_x}{B(\alpha, \beta)}$部分は計算すると$\frac{1}{B(\alpha + x, \beta +(n - x))}$となり、定数部分として省略しています)

今回の例では$Beta(1, 1)$を事前分布と設定し、$n$回の試行の内$x$回成功の観測データを掛ける事で事後分布$Beta(1 + x, 1 + (n - x))$を求めています。
成功確率(期待値)と成功確率(分散)の計算式は$Beta(\alpha, \beta)$の期待値と分散の計算式で、それを使用していました。

さいごに

これの使い道ですが私は回帰等の説明変数として使ったり、2群の差を比較するオッズ比検定として使うことが多いです。
誰かの参考になれば嬉しいです！
間違い等ありましたらご指摘いただけると助かります。

参考文献

この記事は以下の情報を参考にしております。

1
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
3