ベイズ分析の考え方
岡本安晴
2019.7、2025.3
帰無仮説の検定とか最尤法では、データに対してその生成確率を表す確率モデル
が設定される。硬貨を10回投げて7回表が出たとき、表が出る確率
の二項分布のモデルを設定したときは
である。
式（１）の確率モデルの場合、確率はデータに対して与えられている。これを拡張して、データとモデルの組に対して確率を考え、その確率を次式で
で与える。すなわち、式（１）の確率にモデルの確率を掛けてデータとモデルの同時確率とする。モデル（２）はモデル（１）の拡張と考えられる。確率モデルとしては、モデル（１）もモデル（２）も数学的な構造物として数学的に扱うことができる。統計的分析において、どちらがより有用かの問題である（岡本、2018、pp. 72-74）。
一般論として、理論が発展するとそこで用いられる数学的構造も進化したものが用いられる。ニュートン物理学ではユークリッド幾何が用いられたが、現代物理学では非ユークリッド幾何が用いられる。多項式の方程式は、自然数、整数、有理数、実数と解の空間が拡張され、複素数にまで拡張されて常に根をもつようになった。統計分析も、モデル（１）からモデル（２）に進化することにより、分析がより豊かになったといえる。モデル（２）に基づく統計分析は、ベイズ分析と呼ばれている。以下に、ベイズ分析の例を挙げて説明する。説明においてプログラム（スクリプト）が示されることがあるが、スクリプトはプログラムの一部であることが多い。プログラムファイルは圧縮ファイルbysn1files.zipとしてまとめた。
ベイズ分析を簡単な確率モデルで説明する。
例―１
適材適所という言葉があるが、人それぞれに適性があり、適性を活かした仕事に就けることが望ましい。しかし、適性は目に見えるものではないので、何らかの方法で「有り・無し」の判断が行われる。いま、適性の有無の判定が表１に示すようであるとする。
|
表１ 適正判定の確率分布 |
||
|
|
合格 |
不合格 |
|
普通 |
0.1 |
0.9 |
|
適性あり |
0.9 |
0.1 |
これを式（１）の形で書けば次式となる。
（３）
式（１）における「モデル」の部分が、「普通の人に対するモデル」および「適性のある人に対するモデル」となっている。
式（１）に対する代表的な分析法である最尤法を適用してみる。尤度関数は
であるが、「パラメータ」は、式（１）における「モデル」である。「普通」あるいは「適性あり」の条件での確率は、表１あるいは式（３）で与えられるが、これはそれぞれ「普通」あるいは「適性あり」における判定確率のモデルである。式（４）の尤度関数を図示すると図１のようになる。
図１の描画は、以下のスクリプトで行った（リスト１）。
import matplotlib.pyplot as plt
x_range = [0, 1]
# 0: Average, 1: Talent
p_pass = [0.1, 0.9]
p_failure = [0.9, 0.1]
plt.plot(x_range, p_failure, color =
'r', linewidth = 3, label = 'X = Failure')
plt.plot(x_range, p_failure, color =
'r', marker = 'o', markersize = 10)
plt.plot(x_range, p_pass, color = 'b',
linewidth = 3, label = 'X = Pass')
plt.plot(x_range, p_pass, color = 'b',
marker = 'o', markersize = 10)
plt.xlim([-0.25, 1.25])
plt.xticks([0, 1], fontsize = 14)
plt.xlabel('0: Average
1: Talent', fontsize = 14)
plt.ylim([-0.1, 1.1])
plt.yticks([0.0, 1.0], fontsize = 14)
plt.ylabel('Probability(Likelihood)',
fontsize = 14)
plt.legend(loc = 'upper center',
fontsize = 14)
plt.show()
横軸の０が「普通」、１が「適性あり」を表し、青の線分は判定が「合格」であったときの尤度関数、赤の線分は判定が「不合格」であったときの尤度関数である。最尤法では、尤度関数が最大になるパラメータ値が最尤推定値である。判定が「合格」の青い線分では、「適性あり」が最尤推定値であり、判定が「不合格」の赤い線分では、「普通」が最尤推定値である。
ベイズ分析では、式（２）の形でモデルが設定されるが、いまの場合
である。
が、「普通」の人あるいは「適性あり」の人に対する「合格／不合格」の確率であり、
は「普通」の人あるいは「適性あり」の人の確率（比率）
あるいは
を表す。このを事前確率と呼び、判定すなわちデータを条件とする確率
すなわち
を事後確率と呼ぶ。事後分布は、判定が下されたとき、その判定の下でのモデル（いまの場合「普通」および「適性あり」）の確率であり、次式で与えられる。
式（６）において、データが所与の下で考えるので
は定数であり、式（６）をモデルの関数として見るときはを省いて次式
の形で考えることが多い。
式（６）あるいは式（７）より、事後分布は事前分布の影響を受けることがわかる。まず、事前分布において「普通」と「適性あり」の確率が同じである場合について考える。
とおいて事後分布を求め、グラフを描くと図２のようになる。
図２を描くスクリプトは、以下のようである（リスト２）。
Pave = [0.1, 0.9]
# [合格, 不合格]／普通
Ptal = [0.9, 0.1]
# [合格, 不合格]／適性あり
Prior = [0.5, 0.5]
# [普通, 適性あり]
Pjoint = np.zeros((2, 2))
#
合格 不合格
# 普通 [0, 0] [0, 1]
# 才能あり [1, 0] [1, 1]
Pjoint[0, 0] = Prior[0] * Pave[0]
Pjoint[0, 1] = Prior[0] * Pave[1]
Pjoint[1, 0] = Prior[1] * Ptal[0]
Pjoint[1, 1] = Prior[1] * Ptal[1]
print(Pjoint)
Ppass = np.zeros(2)
Ppass[0] = Pjoint[0,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
Ppass[1] = Pjoint[1,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
print(Ppass)
Pfail = np.zeros(2)
Pfail[0] = Pjoint[0,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
Pfail[1] = Pjoint[1,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
print(Pfail)
x_range = [0, 1]
# 0: Average, 1: Talent
plt.plot(x_range, Pfail, color = 'r',
linewidth = 3, label = 'X = Failure')
plt.plot(x_range, Pfail, color = 'r',
marker = 'o', markersize = 10)
plt.plot(x_range, Ppass, color = 'b',
linewidth = 3, label = 'X = Pass')
plt.plot(x_range, Ppass, color = 'b',
marker = 'o', markersize = 10)
plt.xlim([-0.25, 1.25])
plt.xticks([0, 1], fontsize = 14)
plt.xlabel('0: Average
1: Talent', fontsize = 14)
plt.ylim([-0.1, 1.1])
plt.yticks([0.0, 1.0], fontsize = 14)
plt.ylabel('Posterior Distribution',
fontsize = 14)
plt.legend(loc = 'upper center',
fontsize = 14)
plt.show()
事前分布が一様分布、いまの場合、「普通」と「適性あり」が同じ確率のときは、事後分布は尤度関数（図１）と同じ形になる。判定が「合格」ときは「適性あり」の事後確率が高く、判定が「不合格」のときは「普通」の事後確率が高いので、最尤法と同じ結論になる。
「適性あり」の人が100人に１人の割合である場合について考える。これは、事前分布を
と設定することで表せる。このとき、事後分布は図３に示すようになる。
判定が「合格」・「不合格」に関わらず事後分布はほとんど同じで、いずれの場合も「普通」の人である事後確率の方が「適性あり」の確率よりも高い。これでは、「合格」・「不合格」の判定を行う意味がないと言える。
|
表２ 改良版適正判定の確率分布 |
||
|
|
合格 |
不合格 |
|
普通 |
0.0001 |
0.9999 |
|
適性あり |
0.9 |
0.1 |
事前分布が100人に１人の割合の「適性あり」の人を有効に見出すために、判定法（テストの質）の向上を図って、表２に示すものが開発されたとする。「普通」の人が「合格」と判定される確率は１万分の１になっている。「適性あり」の人に対しては「合格」・「不合格」の確率は表１の場合と同じである。表２の特性のテストを、「適性あり」の人が100人に１人の割合の事前分布に実施した場合の事後分布を図４に示す。「合格」判定の人が「適性あり」の人である事後確率が高くなっていて、表２の特性のテストは100人に１人の割合の「適性あり」の人を見出すのに有用であることが分かる。
診断検査では、
は陽性予測値（positive predictive value）と呼ばれ、
は陰性予測値（negative predictive value）と呼ばれているが、具体例を用いた説明が岡本（2014；pp.188-193）にある。
例―２
上の例では、適性の検査結果は「あり／なし」のカテゴリ判定であった。次に、検査結果がカテゴリではなくスコア（連続量）の場合について調べる。
適性について調べるテストが開発され、普通の人は平均100、標準偏差10の正規分布であり、適性ありの人は平均120、標準偏差10の正規分布になっているとする。すなわち、テストのスコア
の分布は次式で与えられる。
ここで、
は、平均
、標準偏差の正規分布確率密度関数である。
このとき、「普通」と「適性あり」をパラメータ
とするモデルの尤度関数は
となる。スコアに対するパラメータの最尤推定値
は
となる。上の式において尤度関数が等しい場合は、これは確率的に0であり、最尤推定値としていずれをとってもよい。
スコアの分布と最尤推定値のグラフを描くスクリプトを以下のように用意した（リスト３）。
muA = 100 # 普通
muB = 120 # 適性あり
sgm = 10
# 標準偏差
x_values = np.arange(70, 150, 0.1)
ScoreA =
scipy.stats.norm.pdf(x_values, muA, sgm)
ScoreB =
scipy.stats.norm.pdf(x_values, muB, sgm)
plt.plot(x_values, ScoreA, label =
'Class A')
plt.plot(x_values, ScoreB, label =
'Class B')
def Lfunc(theta, x):
if theta == 'A':
return
scipy.stats.norm.pdf(x, muA, sgm)
elif theta == 'B':
return scipy.stats.norm.pdf(x, muB, sgm)
else:
raise Exception('Invalid parameter value!')
LE = []
for x in x_values:
LE.append(0 if Lfunc('A', x)
> Lfunc('B', x) else 0.04)
plt.plot(x_values, LE, label = 'MLE')
plt.yticks([])
plt.legend(loc = 'upper left')
plt.title('Distributions and MLE',
fontsize = 16)
plt.show()
Aが「普通」を表し、Bが「適性あり」を表している。最尤推定値は、次のスクリプト
LE = []
for x in x_values:
LE.append(0 if Lfunc('A', x)
> Lfunc('B', x) else 0.04)
により、「普通」を０、「適性あり」を0.04で表している。上のスクリプトを実行すると図５のグラフが描かれる。
緑の折れ線が最尤推定値を表しているが、スコアが110以下のときは最尤推定値は「普通」、110を超えると「適性あり」であることがわかる。
ベイズ分析における事後確率は、以下のようになる。
事前分布として等確率、すなわち「普通」と「適性あり」の比率が同じとした確率を設定したとき、事後確率をスコアの関数として描くスクリプトを以下のように用意した。Aが「普通」、Bが「適性あり」を表している。
prior = {'A': 0.5, 'B': 0.5}
def postP(theta, x, prior):
vA = scipy.stats.norm.pdf(x,
muA, sgm) * prior['A']
vB = scipy.stats.norm.pdf(x,
muB, sgm) * prior['B']
if theta == 'A':
return vA / (vA + vB)
elif theta == 'B':
return vB / (vA + vB)
else:
raise Exception('Invalid parameter value!')
PostA = []
PostB = []
for x in x_values:
PostA.append(postP('A', x,
prior))
PostB.append(postP('B', x,
prior))
plt.plot(x_values, PostA, color = 'b',
label = 'P(A|Score)')
plt.plot(x_values, PostB, color = 'g',
label = 'P(B|Score)')
plt.vlines(110, 0, 1, color = 'y',
linewidth = 2, label = 'Score = 110')
plt.xlabel('Score')
plt.legend()
plt.title('Posterior Distributions at
Scores\n' +\
'Prior Distribution = [0.5, 0.5]', fontsize = 16)
plt.show()
このスクリプトを実行すると、図６のグラフが表示される。
スコアが110より小さいとき「普通」の事後確率の方が高く、スコアが110より高いときは「適性あり」の事後確率の方が高い。事後分布のMAP（Maximum A Posteriori）値（モード値）を推定値とするときは、事前分布が「普通」と「適性あり」の確率が等しいとき、すなわち一様分布のとき、ベイズ分析のMAP推定値と最尤推定値は一致する。
「普通」と「適性あり」が同じ比率であるとは限らない。「適性あり」が100人に１人の割合であるとする。このとき、事前確率は
となる。上のスクリプトを
prior = {'A':
0.99, 'B': 0.01}
とおいて実行すると図７のグラフを得る。
スコアが130ぐらいまでは「普通」の事後確率の方が高い。「適性あり」の事後確率の方が高くなるためには、スコアは130を超えていなければならない。実際、「適性あり」の人の事前確率が100人に１人の割合のときに、スコアが110の人をランダムに集めたとき、「普通」の人と「適性あり」の人の人数は、次のスクリプトで求めることができる（リスト４）。
ndensity = scipy.stats.norm.pdf
p0A = 100 * ndensity(110, 100, 10) *
0.99 /\
(ndensity(110, 100, 10) * 0.99
+ ndensity(110, 120, 10) * 0.01)
print(p0A)
p0B = 100 * ndensity(110, 120, 10) *
0.01 /\
(ndensity(110,
100, 10) * 0.99
+ ndensity(110, 120, 10) * 0.01)
print(p0B)
これを実行すると以下の結果を得る。
99.0
1.0
すなわち、99人が「普通」で、「適性あり」は１人である。スコアが110のときは「普通」と判定するのが適切である。
これに対して、スコアが133に対しては、次のスクリプトで人数を求めることができる。
p1A = 100 * ndensity(133, 100, 10) *
0.99 /\
(ndensity(133, 100, 10) * 0.99
+ ndensity(133, 120, 10) * 0.01)
print(p1A)
p1B = 100 * ndensity(133, 120, 10) *
0.01 /\
(ndensity(133, 100, 10) * 0.99
+
ndensity(133, 120, 10) * 0.01)
print(p1B)
これを実行すると以下の出力を得る。
49.87799649549956
50.122003504500434
100人中「普通」が49.9人、「適性あり」が50.1人の比率で両者ほぼ等しく、スコア133を判定の境目とすることが適切であると考えられる。
図５に示されているスコアの分布では、「普通」と「適性あり」の分布の中間点は110であり、事前分布が一様分布、すなわち「適性あり」の比率が50％のときはこの中間点110が判定の分岐点になっていた。しかし図７は、「適性あり」が100人に１人の比率のときは、事前分布を考慮したベイズ分析ではスコアが110であっても「適性あり」の事後確率は低いことが示されている。これは、図５に示されているように２つの分布に重なりがあるため、事前分布の影響が強く見られたものである。２つの分布の重なりをなくす（無視できる程度にする）と事前分布の影響も目立たなくなる。
いま、テストを改良してスコアの分布の標準偏差が２になったとする（リスト３の後半）。このときの分布のグラフは図８のようになる。
２つの分布の重なりは実質上は見られない。最尤推定値のグラフが図５と同じように描かれているが、これは図５および図８で同じである。しかし、「適性あり」の人の比率をも考慮する、すなわち事前分布もモデルに取り入れているベイズ分析では、図８の分布のテストに対する事後分布は図７とは異なったものとなる。テストスコアの分布が図８に示すものであるとき、「適性あり」の人の比率が100人に１人の場合の事後分布を図９に示す。
スコア110を少し超えたあたりで「適性あり」の事後確率の方が高くなっている。
ベイズ分析を適性検査などに用いると、「適性あり」の人の比率を考慮した判定が可能になる。単純に最尤法を用いた場合は、事前分布における偏りを無視した不適切な判断が行われる危険がある。ベイズ分析による事前分布を考慮した判定が望まれる。
参考文献
岡本安晴（2018）．心理学における量的研究法・質的研究法の統合．日本情報教育学会誌、１、70-77.
import matplotlib.pyplot as plt
x_range = [0, 1]
# 0: Average, 1: Talent
p_pass = [0.1, 0.9]
p_failure = [0.9, 0.1]
plt.plot(x_range, p_failure, color =
'r', linewidth = 3, label = 'X = Failure')
plt.plot(x_range, p_failure, color =
'r', marker = 'o', markersize = 10)
plt.plot(x_range, p_pass, color = 'b',
linewidth = 3, label = 'X = Pass')
plt.plot(x_range, p_pass, color = 'b',
marker = 'o', markersize = 10)
plt.xlim([-0.25, 1.25])
plt.xticks([0, 1], fontsize = 14)
plt.xlabel('0: Average
1: Talent', fontsize = 14)
plt.ylim([-0.1, 1.1])
plt.yticks([0.0, 1.0], fontsize = 14)
plt.ylabel('Probability(Likelihood)',
fontsize = 14)
plt.legend(loc = 'upper center',
fontsize = 14)
plt.show()
import matplotlib.pyplot as plt
import numpy as np
Pave = [0.1, 0.9]
# [合格, 不合格]／普通
Ptal = [0.9, 0.1]
# [合格, 不合格]／適性あり
Prior = [0.5, 0.5]
# [普通, 適性あり]
Pjoint = np.zeros((2, 2))
#
合格 不合格
# 普通 [0, 0] [0, 1]
# 才能あり [1, 0] [1, 1]
Pjoint[0, 0] = Prior[0] * Pave[0]
Pjoint[0, 1] = Prior[0] * Pave[1]
Pjoint[1, 0] = Prior[1] * Ptal[0]
Pjoint[1, 1] = Prior[1] * Ptal[1]
print(Pjoint)
Ppass = np.zeros(2)
Ppass[0] = Pjoint[0,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
Ppass[1] = Pjoint[1,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
print(Ppass)
Pfail = np.zeros(2)
Pfail[0] = Pjoint[0,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
Pfail[1] = Pjoint[1,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
print(Pfail)
x_range = [0, 1]
# 0: Average, 1: Talent
plt.plot(x_range, Pfail, color = 'r',
linewidth = 3, label = 'X = Failure')
plt.plot(x_range, Pfail, color = 'r',
marker = 'o', markersize = 10)
plt.plot(x_range, Ppass, color = 'b',
linewidth = 3, label = 'X = Pass')
plt.plot(x_range, Ppass, color = 'b',
marker = 'o', markersize = 10)
plt.xlim([-0.25, 1.25])
plt.xticks([0, 1], fontsize = 14)
plt.xlabel('0: Average
1: Talent', fontsize = 14)
plt.ylim([-0.1, 1.1])
plt.yticks([0.0, 1.0], fontsize = 14)
plt.ylabel('Posterior Distribution',
fontsize = 14)
plt.legend(loc = 'upper center',
fontsize = 14)
plt.show()
Prior = [0.99, 0.01] # [Average, Talent]
Pjoint[0, 0] = Prior[0] * Pave[0]
Pjoint[0, 1] = Prior[0] * Pave[1]
Pjoint[1, 0] = Prior[1] * Ptal[0]
Pjoint[1, 1] = Prior[1] * Ptal[1]
print(Pjoint)
Ppass = np.zeros(2)
Ppass[0] = Pjoint[0,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
Ppass[1] = Pjoint[1,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
print(Ppass)
Pfail = np.zeros(2)
Pfail[0] = Pjoint[0,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
Pfail[1] = Pjoint[1,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
print(Pfail)
x_range = [0, 1]
# 0: Average, 1: Talent
plt.plot(x_range, Pfail, color = 'r',
linewidth = 3, label = 'X = Failure')
plt.plot(x_range, Pfail, color = 'r',
marker = 'o', markersize = 10)
plt.plot(x_range, Ppass, color = 'b',
linewidth = 3, label = 'X = Pass')
plt.plot(x_range, Ppass, color = 'b',
marker = 'o', markersize = 10)
plt.xlim([-0.25, 1.25])
plt.xticks([0, 1], fontsize = 14)
plt.xlabel('0: Average
1: Talent', fontsize = 14)
plt.ylim([-0.1, 1.1])
plt.yticks([0.0, 1.0], fontsize = 14)
plt.ylabel('Posterior Distribution',
fontsize = 14)
plt.legend(loc = 'upper center',
fontsize = 14)
plt.show()
Pave = [0.0001, 0.9999]
# [合格, 不合格]／普通
Prior = [0.99, 0.01]
# [普通, 才能あり]
Pjoint[0, 0] = Prior[0] * Pave[0]
Pjoint[0, 1] = Prior[0] * Pave[1]
Pjoint[1, 0] = Prior[1] * Ptal[0]
Pjoint[1, 1] = Prior[1] * Ptal[1]
print(Pjoint)
Ppass = np.zeros(2)
Ppass[0] = Pjoint[0,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
Ppass[1] = Pjoint[1,0] / (Pjoint[0, 0]
+ Pjoint[1, 0])
print(Ppass)
Pfail = np.zeros(2)
Pfail[0] = Pjoint[0,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
Pfail[1] = Pjoint[1,1] / (Pjoint[0, 1]
+ Pjoint[1, 1])
print(Pfail)
x_range = [0, 1]
# 0: Average, 1: Talent
plt.plot(x_range, Pfail, color = 'r',
linewidth = 3, label = 'X = Failure')
plt.plot(x_range, Pfail, color = 'r',
marker = 'o', markersize = 10)
plt.plot(x_range, Ppass, color = 'b',
linewidth = 3, label = 'X = Pass')
plt.plot(x_range, Ppass, color = 'b',
marker = 'o', markersize = 10)
plt.xlim([-0.25, 1.25])
plt.xticks([0, 1], fontsize = 14)
plt.xlabel('0: Average
1: Talent', fontsize = 14)
plt.ylim([-0.1, 1.1])
plt.yticks([0.0, 1.0], fontsize = 14)
plt.ylabel('Posterior Distribution',
fontsize = 14)
plt.legend(loc = 'upper center',
fontsize = 14)
plt.show()
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
muA = 100 # 普通
muB = 120 # 適性あり
sgm = 10
# 標準偏差
x_values = np.arange(70, 150, 0.1)
ScoreA = scipy.stats.norm.pdf(x_values,
muA, sgm)
ScoreB =
scipy.stats.norm.pdf(x_values, muB, sgm)
plt.plot(x_values, ScoreA, label =
'Class A')
plt.plot(x_values, ScoreB, label =
'Class B')
def Lfunc(theta, x):
if theta == 'A':
return scipy.stats.norm.pdf(x, muA, sgm)
elif theta == 'B':
return scipy.stats.norm.pdf(x, muB, sgm)
else:
raise Exception('Invalid parameter value!')
LE = []
for x in x_values:
LE.append(0 if Lfunc('A', x)
> Lfunc('B', x) else 0.04)
plt.plot(x_values, LE, label = 'MLE')
plt.yticks([])
plt.legend(loc = 'upper left')
plt.title('Distributions and MLE',
fontsize = 16)
plt.show()
prior = {'A': 0.5, 'B': 0.5}
def postP(theta, x, prior):
vA = scipy.stats.norm.pdf(x,
muA, sgm) * prior['A']
vB = scipy.stats.norm.pdf(x,
muB, sgm) * prior['B']
if theta == 'A':
return vA / (vA + vB)
elif theta == 'B':
return vB / (vA + vB)
else:
raise Exception('Invalid parameter value!')
PostA = []
PostB = []
for x in x_values:
PostA.append(postP('A', x,
prior))
PostB.append(postP('B', x,
prior))
plt.plot(x_values, PostA, color = 'b',
label = 'P(A|Score)')
plt.plot(x_values, PostB, color = 'g',
label = 'P(B|Score)')
plt.vlines(110, 0, 1, color = 'y',
linewidth = 2, label = 'Score = 110')
plt.xlabel('Score')
plt.legend()
plt.title('Posterior Distributions at
Scores\n' +\
'Prior Distribution = [0.5, 0.5]', fontsize = 16)
plt.show()
prior = {'A': 0.99, 'B': 0.01}
print(prior['A'])
PostA = []
PostB = []
for x in x_values:
PostA.append(postP('A', x,
prior))
PostB.append(postP('B', x,
prior))
plt.plot(x_values, PostA, color = 'b',
label = 'P(A|Score)')
plt.plot(x_values, PostB, color = 'g',
label = 'P(B|Score)')
plt.vlines(110, 0, 1, color = 'y', linewidth
= 2, label = 'Score = 110')
plt.xlabel('Score')
plt.legend()
plt.title('Posterior Distributions at
Scores\n' +\
'Prior Distribution = [0.99, 0.01]', fontsize = 16)
plt.show()
sgm = 2 # 標準偏差を２に設定
ScoreA = scipy.stats.norm.pdf(x_values,
muA, sgm)
ScoreB =
scipy.stats.norm.pdf(x_values, muB, sgm)
plt.plot(x_values, ScoreA, label =
'Class A')
plt.plot(x_values, ScoreB, label =
'Class B')
def Lfunc(theta, x):
if theta == 'A':
return scipy.stats.norm.pdf(x, muA, sgm)
elif theta == 'B':
return scipy.stats.norm.pdf(x, muB, sgm)
else:
raise Exception('Invalid parameter value!')
LE = []
for x in x_values:
LE.append(0 if Lfunc('A', x)
> Lfunc('B', x) else ScoreB.max())
plt.plot(x_values, LE, label = 'MLE')
plt.yticks([])
plt.legend(loc = 'upper left')
plt.title('Distributions and MLE\n' +\
'sgm = 2, Prior = [0.99,
0.01]', fontsize = 16)
plt.show()
def postP(theta, x, prior):
vA = scipy.stats.norm.pdf(x,
muA, sgm) * prior['A']
vB = scipy.stats.norm.pdf(x,
muB, sgm) * prior['B']
if theta == 'A':
return vA / (vA + vB)
elif theta == 'B':
return vB / (vA + vB)
else:
raise Exception('Invalid parameter value!')
PostA = []
PostB = []
for x in x_values:
PostA.append(postP('A', x,
prior))
PostB.append(postP('B', x,
prior))
plt.plot(x_values, PostA, color = 'b',
label = 'P(A|Score)')
plt.plot(x_values, PostB, color = 'g',
label = 'P(B|Score)')
plt.vlines(110, 0, 1, color = 'y',
linewidth = 2, label = 'Score = 110')
plt.xlabel('Score')
plt.legend()
plt.title('Posterior Distributions at
Scores\n' +\
'sgm = 2, Prior Distribution
= [0.99, 0.01]', fontsize = 16)
plt.show()
print(prior)
import scipy.stats
import numpy as np
ndensity = scipy.stats.norm.pdf
p0A = 100 * ndensity(110, 100, 10) *
0.99 /\
(ndensity(110, 100, 10) * 0.99
+ ndensity(110, 120, 10) * 0.01)
print(p0A)
p0B = 100 * ndensity(110, 120, 10) *
0.01 /\
(ndensity(110, 100, 10)
* 0.99
+ ndensity(110, 120, 10) * 0.01)
print(p0B)
import scipy.optimize
def fpoint(p):
v =( 50 - 100 * ndensity(p,
100, 10) * 0.99 /\
(ndensity(p, 100, 10) * 0.99
+
ndensity(p, 120, 10) * 0.01)) ** 2
return v
z = np.array([130.0])
res = scipy.optimize.minimize(fpoint,
z, method = 'Nelder-Mead',
options = {'xatol': 1.0e-7, 'fatol': 1.0e-12})
print(res.x)
p1A = 100 * ndensity(133, 100, 10) *
0.99 /\
(ndensity(133, 100, 10) * 0.99
+ ndensity(133, 120, 10) * 0.01)
print(p1A)
p1B = 100 * ndensity(133, 120, 10) *
0.01 /\
(ndensity(133, 100, 10) * 0.99
+ ndensity(133, 120, 10) * 0.01)
print(p1B)