確率
確率とその数学的基礎に対するある種の理解を欠いたままで、データサイエンスを行うのは困難です。
目標のために、現実の事象の不確実さを定量化する方法の１つとして確率を考えることは重要です。表記はP(E)で、「事象Eの発生する確率」を表します。
確率論を使ってモデルを構築し、確率論を使ってモデルを評価します。確率はあらゆる場面で使用します。
従属と独立
大まかに言うと、Eの発生に関する何らかの情報がFの発生についての情報を与える場合、２つの事象EとFは従属関係にあるといいます。
例えば、コインを２回投げたとします。１回目でオモテが出たか否かを知っていたとしても、２回目の結果には影響がありません。そのためこれらの事象は独立です。
一方、１回目がオモテであったことは２回とも裏が出る事象について明らかな情報を与えます。この事象は従属です。
数学的に、事象EとFが独立である場合、その両方が発生する確率は、それぞれの事象が発生する確率の積となります。
P(E, F) =P(E)P(F)
この例では、１回目に表が出る確率は1/2であり、２回とも裏が出る確率は1/4です。しかし、１回目に表が出て、２回目とも裏がでる確率は０となります。
条件付き確率
独立なら定義は、P(E, F) =P(E)P(F)
両者が独立である必要がない（そして、Fの発生する確率が0でない）のであれば、FにおけるEの条件付き確率は次の式で定義できます。
P(E|F) = P(E, F)/P(F)
これは、事象Fが発生したことを知っている状況でEが発生する確率と考える事ができます。次のように置き換えられることもあります。
P(E, F) = P(E|F)P(F)
EとFが独立なら、次の式が成り立ちます。
P(E|F) = P(E)
これはFが発生したことを知っていても、Eの発生に何ら影響を与えないことを数学的に表したものです。
よく知られた例として、2人の子供（性格はわからない）がいる家族を考えます。次のことを仮定した場合、
- 子供がそれぞれ男の子か女の子である可能性は等しい
- 2人目の子供の性別は、1人目の子供の性別とは独立している。
女の子がいない確率は1/4、1人が女の子で1人が男の子である確率は1/2、2人とも女の子である確率は1/4となります。
1人メが女の子である（G）場合に、2人ともおんなのこである（B)確率はどうなるでしょうか。事象BかつGは事象Bに等しいので、条件付き確率の定義によると次のようになります。
P(B|G) = P(B,G)/P(G) =P(B)/P(G) =1/2
おそらくこの理解は直感的理解と一致してます。
同じく、少なくとも1人が女の子である（L）に、2人とも女の子である確率も求められます。値は異なります。
先程と同様に、事象BかつLは事象として等しいので、次のようになります。
P(B|L) = P(B,L)/P(L) =P(B)/P(L) = 1/3
これはどういうことでしょうか。つまり、少なくとも1人が女子出会った場合、2人とも女の子である場合よりも男の子と女の子が1人づつである可能性は2倍あるということになります。
家族構成を大量に生成して、この状況を確認します。
def random_kid():
return random.choice(["boy", "girl"])both_girls = 0
older_girl = 0
either_girl = 0rendom.seed(0)
for _ in range(1000):
yonger = random_kid()
older = random_kid()
if older == "girl":
older_girl += 1
if older == "girl" and younger == "girl":
both_girls += 1
if older == "girl" or younger =="girl":
eigher_girl += 1print(both_girls / older_girl) #2人とも女の子｜1人目が女の子 #0.2
print(both_girls / either_girl) #2人とも女の子｜どちらか1人が女の子 #0.3
ベイズの定理
データサイエンスで最良のパートナーの1つがベイズの定理です。
これは、条件付き確率を裏返しにする手法です。
事象Fが発生した状況で、それとは独立した事象Eが起きる確率を求めるとしましょう。しかし、事象Eが発生した状況で、事象Fの発生する確率だけが既知であるとします。次の式になります。
P(E|F) = P(E, F)/P(F) = P(F) = P(F|E)P(E)/P(F)
事象Fは、相互に排他的な2つの事象「FかつE」と「FかつnotE」に分割できます。「not E」（つまり、「Eが発生しない」）を¬Eと表記して次の式になります。
P(F) =P(F, E)+P(F, ¬E）
以上より、次の式が導き出されます。
P(E|F) = P(F|E)P(E)/[P(F|E)P(E)+P(F|¬E)P(¬E)]
これがベイズの定理です。
例えば、10,000人あたり1人が発症する疾患があるとしましょう。そしてこの疾患を99%の正確さで検出できる検査があるとします。（陽性と陰性）
検査の陽性が意味することは何でしょう。ここで、「検査が陽性である」事象をTと、「疾患を持っている」事象をDとします。ベイズの定理では検査が陽性であった場合に、疾患を持っている確率は以下になります。
P(D|T) = P(T|D)P(D) / [P(T|D)P(T) + P(T|¬D)P(¬D)]
ここで、P(T|D)つまり「疾患を持つ人が検査で陽性になる確率」は0.99%であることが分かっています。ある人が疾患を持つ確率P(D)は、1/10,000 = 0.00001。
疾患を持っていないが、テストで陽性となる確率P(T|¬D)は0.01。そしてある人が疾患を持たない確率P(¬D)は、0.9999です。これらを代入すると、
P(D|T) = 0.98%
つまり、検査で陽性が出た人が実際に疾患を持っている確率は1%以下であることになります。
確率変数
確率変数とは、確率分布に関連づいた値を持つ変数です。コインのオモテがでたら1、裏が出たら0となるような確率変数が非常に単純な例です。
より複雑なものは、コインを10回投げた際にオモテが出た回数をとるものや、range(10)から均等に確からしさで値を取り出すものなどが考えられます。
関連する分布は、確率変数が取りうる値それぞれの起こりやすさを与えます。コイン投げでは、値が0となる確率が0.5、値が1となる確率が0.5です。range(10)は0から9までそれぞれの正確さが0.1となる分布を持ちます。
確率変数の値を確率の重み付き平均で計算される期待値を話題にすることがある。例えば、コイン投げの期待値は1/2（0*1/2 + 1*1/2）であり、range(10)の期待値は4.5になります。
他の事象と同様に条件付きの事象についても定義できます。先程の子供のレを使うと、Xを女の子の数を表す確率変数とすると、Xが0の場合は1/4 、1の場合は1/2、2の場合は1/4。2人のうち数なくとも1人が女の子であった場合の女の子の数を新しい確率変数Yとして定義できます。Yが1の場合の確率は2/3、Yが2の場合は1/3です。1人目の子供が女の子出会った場合の、女の子の数を確率変数Zとすると、Zが１の場合の確率は1/2、Zが2の場合は1/2となります。
このあと、多くの場面で特別な配慮を払わずに確率変数を使うことになります。しかし、注意深く掘り下げてみれば、そこに確率変数が使われていることが分かるでしょう。
連続確率分布
コインなげは離散型分布に相当し、離散的な結果に対する生の確率に関連付けられます。ときに連続した結果の分布をモデル化する必要性が生じます。（この結果は常に実数で得られますが、すべて実生活上の事象を表しているわけではありません）例えば、一様分布派、0から1のすべての数に対して等しい重みを与えます。
0から1の間には無限の数が存在することを考えると、ここの数に対する重みは0とする必要があります。このため、連続分布は確率密度関数（probability density function:pdf）で表現し、確率変数がある範囲の値をとる確率は、その範囲で確率密度関数を積分することで得られます。
一様分布の密度関数は次のように表します。
def uniform_pdf(x):
return 1 if x >= 0 and x < 1 else 0この分布に従う確率変数が0.2から0.3の値を取る確率は1/10となります。Pythonのrandom.random()は一様分布で（疑似）乱数です。
確率変数の値がある値以下となる確率を表す累積分布関数（cumulative distribution function:cdf）の方を使うこともたまにあります。
def uniform_cdf(x):
if x < 0: return 0 #一様分布は0を下回らない
elif x < 1: return x #例えばP(X <= 0.4) = 0.4となる
else: return 1正規分布
正規分布は、あらゆる分布の中で最も重要な存在です。この釣鐘型の分布は2つのパラメータ、平均µと標準偏差σで定義されます。平均は釣鐘の中心を表し、標準偏差は釣り鐘の横幅を表します。
確率密度関数は次の式で与えられます。
実装は次のようになります。
def nomal_pdf(x, mu=0, sigma=1):
sqrt_two_pi = math.sqrt(2 * math.pi)
return (math.exp(-(x - mu) ** 2 / 2 / sigma ** 2) / (sqrt_two_pi * sigma))グラフにするとこのようになります。
µ=0, σ=1の場合を、標準正規分布といいます。
Zが標準正規分布に従う確率変数であった場合、
X = σZ + µ
確率変数Xは平均µ、標準偏差σの正規分布となります。
逆にXが平均µ標準偏差σの正規分布に従う確率変数であるなら、
Z = (X — µ) / σ
確率変数Zは標準正規分布に従います。
正規分布の累積分布関数は初歩的な方法では実装できませんが、Pythonのmath.erfを使えば次のようになります。
def normal_cdf(x, mu=0, sigma=1):
return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2
特定の確率となる値を見つけるためにnormal_cdfの逆関数が必要となる場合があります。逆関数を単純に計算する方法はありませんが、nomal_cdfは連続で単調増加であるため、二分探索が使えます。
def inverse_normal_cdf(p, mu=0, sigma=1, tokerance=0.00001):
#二分探索を用いて、逆関数の近似値を計算する#標準正規分布でない場合、標準正規分布からの差分を求める
if mu != 0 or sigma != 1:
return mu + sigma * inverse_nomal_cdf(p, tolerance=tolerance)low_z, low_p = -10.0, 0 #nomal_cdf(-10)は、0に近い値である
hi_z, hi_p = 10.0, 1 #nomal_cdf(10)は1に近い値であるwhile hi_z - low_z > tolarence:
mid_z = (low_z + hi_z) / 2 #中央の値および
mid_p = nomal_cdf(mid_z) #その地点でのcdfの値
if mid_p < p:
#中央値はまだ小さいので更に上を使う
low_z, low_p = mid_z, mid_p
elif mid_p > p:
#中央値はまだ大きいのでさらに下を使う
hi_z, hi_p = mid_z, mid_p
else:
break
return mid_z
この関数は、目的の確率に十分近づくまでZの区間の二等分を繰り返します。
中心極限定理
正規分布が有用である1つの理由が中心極限定理です。
簡潔に説明すると、常に多数の独立で同一の分布に従う確率変数の平均として定義される確率変数が、およそ正規分布となるというのが中心極限定理です。
例えば、平均がµ、標準偏差がσの確率変数 x1, … , xn があるとします。ｎは十分に大きいものとします。この時、1/n(x1 + … + xn)はおよそ平均μ、標準偏差σ/√nの正規分布となります。次の式は同様に平均0、標準偏差1の正規分布です。
((x1 + … + xn) / µn) / σ√n
これを簡単に説明するには、nとpで表される二項確率変数を見ましょう。確率pで1、確率(1-p)で0となるn個の独立した確率変数Bernoulli(p)を合計したものがBinomial(n, p)確率変数です。
def bernoulli_trial(p):
return 1 if random.randm() < p else 0def binomial(n, p):
return sum(benoulli_trial(p) for _ in range(n))
Bernoulli(p)の平均はp,標準偏差は(p(1-p)**0.5です。中心極限定理に夜と、nが大きければBinomial(n, p)はおおよそ平均 µ=np,標準偏差σ = (np(1-p)**0.5の正規分布となります。並べてプロットすれば、その類似性が把握できるできます。
def make_hist(p, n, num_points):
data = [binomial(n, p) for _ in range(num_points)]
histogram = Counter(data)
plt.bar([x - 0.4 for x in histogram.keys()],
[v / num_points for v in histgram.values()],
0.8,
color = '0.75')
mu = p * n
sigma = math.sqrt(n * p * (1 - p))
#正規分布の近似を折れ線グラフでプロットする
xs = range(min(data), max(data) + 1)
ys = [nomaL_cdf(i + 0.5, mu, sigma) - nomal_cdf(i - 0.5, mu, sigma) for i in xs]
plt.plot(xs, ys)
plt.title("Binomial istribution vs. Normal Approximation")
plt.show()
この近似から得られる教訓は、歪みがないとされているコインを１００回投げた際に、オモテは６０回以上出る確率を求めるには、Nomal(50, 5)が６０市場となる確率を求めれば良いことになります。これはBinomial(100, 0.5)の累積分布関数を計算するよりも簡単です。
補足
- scipy.statsは大抵の確率分布に対するPDF, CDFをもとめる関数を提供しています。