2017-11-29

幾何分布と負の二項分布

python 確率

幾何分布 (geometric distribution)

あっかりーんという存在感が薄い人がいたとします。あっかりーんが他人とすれ違うときに認識される確率が毎回30%だったとして、 3人目(3回目)に初めて認識される確率はいくらでしょうか。

この場合、2人に認識されない確率 $(1-0.3)^{2}$ と3人目に認識される確率 $0.3$ をかけて以下のように

$0.3\cdot(1-0.3)^{2}$

と計算すると思います。

これを抽象的に言うと、ベルヌーイ試行（※二項分布の記事参照）において初めてある事象（当たりなど）が起こるまでに要した回数をXとして、 $x$ 回目が当たり(確率 $p$ )で $x-1$ 回目までははずれ(確率 $1-p$ )であるので確率関数は

$f(x) = p(1-p)^{x-1} \ (x=1,2,3,···)$

となります。以下、python(+ jupyter notebook)であっかりーんが認識される確率 $p$ を0.1~0.7まで振って10回目までのレンジで分布をみてみました。

# geometric distribution
def geom_(x,p):
    return (1-p)**(x-1)*p

range_x = range(1, 11)

dict_geom= {}
for p in [n*0.1 for n in range(1,8)]:
    dict_geom[p]= [geom_(x, p) for x in range_x]

import pandas as pd

df_geom=pd.DataFrame(dict_geom)
df_geom

	0.1	0.2	0.3	0.4	0.5	0.6	0.7
0	0.100000	0.200000	0.300000	0.400000	0.500000	0.600000	0.700000
1	0.090000	0.160000	0.210000	0.240000	0.250000	0.240000	0.210000
2	0.081000	0.128000	0.147000	0.144000	0.125000	0.096000	0.063000
3	0.072900	0.102400	0.102900	0.086400	0.062500	0.038400	0.018900
4	0.065610	0.081920	0.072030	0.051840	0.031250	0.015360	0.005670
5	0.059049	0.065536	0.050421	0.031104	0.015625	0.006144	0.001701
6	0.053144	0.052429	0.035295	0.018662	0.007812	0.002458	0.000510
7	0.047830	0.041943	0.024706	0.011197	0.003906	0.000983	0.000153
8	0.043047	0.033554	0.017294	0.006718	0.001953	0.000393	0.000046
9	0.038742	0.026844	0.012106	0.004031	0.000977	0.000157	0.000014

matplotlibで可視化してみましょう。

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(range_x, df_geom, 'o')
plt.legend(df_geom.columns, title='p =')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

scipyパッケージにもあるので同じく確認してみましょう。

from scipy.stats import geom

list_geom= [geom(p).pmf(range(1,11)) for p in [n*0.1 for n in range(1,8)]]
df_geom_=pd.DataFrame(list_geom)
plt.plot(range_x, df_geom_.T, 'o')

負の二項分布 (negative binomial distribution)

幾何分布を拡張して、 $r$ 回当たりがでるまで試行を行なったとします。当たりが $r$ 回でるまでに要した回数を $X$ とすると確率関数は

$f(x) = {}_{x-1} C _{r-1} p^r(1-p)^{x-r}$

となり、これは負の二項分布(negative binomial distribution)またはPascal分布と呼ばれています。 $r=1$ のときは幾何分布になります。

また、 $r$ 回当たりがでるまでのはずれの回数を $X$ とした確率関数は

$f(x) = {}_{r+x-1} C _x p^r(1-p)^x$

となります。式の通り負の二項分布のパラメータは、当たりの確率 $p$ と当たり回数 $x$ の2つです。
(例)
あっかりーんが認識される確率 $p$ が0.3のとき、5人に認識される( $=r$ )まで出歩いたときの確率分布はどのような形になるでしょうか。以下でみてみましょう。

# r回当たるまで要した回数をX

from scipy.special import comb
def nebinom_s(x, r, p):
    return comb(x-1,r-1)*p**r*(1-p)**(x-r)

range_x2 = range(1,40)
r = 5     # 当てる回数
p = 0.3 # 当たる確率

nbi_s= [nebinom_s(x, r, p) for x in range_x2]
plt.plot(range_x2, nbi_s, 'o')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

# r回当たるまでのはずれの回数をX

def nebinom_f(x, r, p):
    return comb(r+x-1,x)*p**r*(1-p)**x

nbi_f= [nebinom_f(x, r, p) for x in range_x2]
plt.plot(range_x2, nbi_f, 'o')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

scipyパッケージにもあるのでそれでも試してみましょう。scipyのはデフォではずれの回数を $X$ とした確率関数の式を用いています。

from scipy.stats import nbinom

plt.plot(range_x2, nbinom.pmf(range_x2, r,p), 'o')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

確率 $p$ が0.3のとき、当たり回数(認識した人の人数) $r$ を1~10まで振ってみましょう。

for temp_r in range(1,11):
    plt.plot(range_x2, nbinom.pmf(range_x2, temp_r, p), 'o', label=temp_r)
plt.legend(title='r =')
plt.ylabel('probability', fontsize=15)
plt.xlabel('x', fontsize=15)

あっかりーんを認識する人が増えるとそれまでの認識しなかった人が増えて行くことが分かりました(直感的にそうですね)。それと同時に分布の裾野が広がって曖昧さも増加しているようです。

2017-10-29

二項分布(Binomial distribution)

python 確率

f:id:udonista:20171029040332p:plain

今ここにSレアであるあっかりーんが当たるガチャがあるとします。ネットの情報からあっかりーんが当たる確率は75%(!?)だとの情報を運良くゲットしたとしましょう。あなたはこのSレアあっかりーんがとても欲しくて、1回あたり200円のそのガチャに5000円をつっこもうと考えています。
〜〜
ガチャを回し、当たり・はずれのような結果が2パターンになるものをベルヌーイ試行と言います。
※各回の試行は独立です

当たる確率を $p$ , はずれる確率を $(1-p)$ であるガチャを $n$ 回引くとすると、各あっかりーん数 $x$ の得られる確率は

${\begin{align} f(x)&= {}_n C _xp^x(1-p)^{n-x} \\ \end{align}}$

で与えられます。この分布を二項分布(Binomial distribution)と呼びます。ベルヌーイ試行をn回行なったときにある事象(当たりなど)が何回起こるかの確率分布です。この一般化された式を見ての通り、二項分布は試行回数nと当たる確率pにより決まることが分かります。つまり、二項分布のパラメータは $n$ と $p$ なので $Bi(n,p)$ と書いたりします。

それでは最初に書いたあっかりーんガチャを例に、二項分布をみてみましょう。あっかりーんが当たる確率は75%、引く回数は5000円/200円で25回です。
つまり、 $Bi(25, 0.75)$ なので、このときの確率分布は、

$f(x) = {}_{25} C _x0.75^x(1-0.75)^{25-x}$

となります。以下のようにpythonのscipyパッケージを用いると簡単に求めることができます。

from scipy.stats import binom

p = 0.75             # あっかりーんが当たる確率
n = 5000/200         # 200円のガチャに5000円をつっこむので25回まわせる
x = range(int(n)+1)  # 0~25回まで当たるときに着目

binomial = binom.pmf(x,n,p)

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(x, binomial, 'o')
plt.xlabel('number of success', fontsize=15)
plt.ylabel('probability', fontsize=15)

f:id:udonista:20171029040549p:plain

分布がでました。
では、25回まわした時に得られるあっかりーん数 $x$ のうち、最も確率が高いのはなんでしょうか？

binomial.argmax()

考察

このガチャを25回まわすと19個ものあっかりーんが得られる可能性が最も高いことがわかりました。よかったですね。
ちなみに二項分布は期待値 $\mu=np$ , 分散 $\sigma^{2} = np(1-p)$ となります。

n*p

18.75

n*p*(1-p)

4.6875

※平均分散の導き方
・平均

$\begin{align} \mu&= \sum_{x=0}^n x\cdot f(x) \\ &= \sum_{x=0}^n x\cdot \frac{n!}{x!(n-x)!} p^x(1-p)^{n-x} \\ &= \sum_{x=0}^n x\cdot \frac{n!}{x(x-1)!(n-x)!} p^x(1-p)^{n-x} \\ &= \sum_{x=0}^n \frac{n(n-1)!}{(x-1)!(n-x)!} p\cdot p^{x-1}(1-p)^{n-x} \\ &= np\sum_{x=0}^n \frac{(n-1)!}{(x-1)!(n-x)!} p^{x-1}(1-p)^{n-x} \\ &= np \\ \end{align}$

※ $x! = x(x-1)!$
5行目の和の部分は $Bi(n-1, p)$ の確率分布の和であるので1になります。

・分散

$\begin{align} \sigma^{2}&= \sum_{x=0}^n (x-\mu)^{2}\cdot f(x)\\ &= \sum_{x=0}^n (x^{2}-2\mu x +\mu^{2})\cdot f(x)\\ &= \sum_{x=0}^n x^{2}\cdot f(x)-2\mu\sum_{x=0}^n x \cdot f(x) + \mu^{2}\sum_{x=0}^n f(x)\\ &= \sum_{x=0}^n x^{2}\cdot f(x) -2\mu \cdot \mu +\mu^{2}\cdot 1 \\ &= \sum_{x=0}^n x^{2}\cdot f(x) -(\sum_{x=0}^n x\cdot f(x))^{2}\\ &= n(n-1)p^{2}+np-(np)^{2}\\ &= (np)^{2} -np^{2} + np -(np)^{2}\\ &= np(1-p)\\ \end{align}$

※

$\begin{align} \sum_{x=0}^n x^{2}\cdot f(x)&= \sum_{x=0}^n (x(x-1)+x) \frac{n!}{x!(n-x)!}p^{x}(1-p)^{n-x}\\ &= \sum_{x=0}^n x(x-1) \frac{n(n-1)(n-2)!}{x(x-1)(x-2)!(n-x)!}p^2 \cdot p^{x-2}(1-p)^{n-x}\\\ &+ \sum_{x=0}^n x\cdot \frac{n(n-1)!}{x(x-1)!(n-x)!}p \cdot p^{x-1}(1-p)^{n-x}\\ &= n(n-1)p^{2}+np \end{align}$

※ $x^{2}=x(x-1)+x$

おまけ
当たる確率0.1%のSSレア京子ちゃんを狙った場合：

p_kyoko = 0.1/100     #京子ちゃんが当たる確率 

binomial_k = binom.pmf(x,n,p_kyoko)

plt.plot(x, binomial_k, 'o')
plt.xlabel('number of success', fontsize=15)
plt.ylabel('probability', fontsize=15)

f:id:udonista:20171029040623p:plain

なんと25回まわしたぐらいではあたらなそうです。では10倍の5万円つぎこんだらどうでしょうか

n = 50000/200        # 200円のガチャに50000円をつっこむ
x = range(int(n)+1)  # 0~250回まで当たるときに着目

binomial_k = binom.pmf(x,n,p_kyoko)

plt.plot(x, binomial_k, 'o')
plt.xlabel('number of success', fontsize=15)
plt.ylabel('probability', fontsize=15)

f:id:udonista:20171029040639p:plain

5回までのところを拡大してみましょう。

plt.plot(range(6), binomial_k[:6], 'o')
plt.xlabel('number of success', fontsize=15)
plt.ylabel('probability', fontsize=15)

f:id:udonista:20171029040653p:plain

京子ちゃんを当てに行く場合、250回まわしても1回も当たらない確率は80%弱もあることが分かりました。 250回まわして当たる回数の期待値は、

n*p_kyoko

 0.25

分布をみると1回も当たらない確率が80%近いのに、期待値を計算すると当たる確率は25%でした。
このような偏った分布の場合では、平均値の計算は適さないです。

今回は以上になります。

2017-10-28

摂氏、華氏の換算

python

f:id:udonista:20171029035957p:plain

華氏はアメリカで使われていますね。ですが実はアメリカ以外の国では今やほとんど使われていないようです(華氏, wikipedia, 10.27.2017)
華氏の考案者はドイツの物理学者Fahrenheitさんなのですがドイツでは華氏ではなく摂氏が使用されているみたいです。

華氏は真水の凝固点を32度(32°F)、沸点を212度(212°F)としていて、摂氏0-100度を180等分しています。
つまり摂氏と華氏の1度の幅は異なっていて、その関係(比)は

100/180 = 5/9

だから、華氏が1度上がると摂氏は5/9度増えて、摂氏が1度上がると華氏は9/5度増えるわけですね。
では摂氏(C)→華氏(F)、華氏(F)→摂氏(C)の換算式を考えてみます。

${ \displaystyle \begin{align} F=\frac{9}{5}C\\ C=\frac{5}{9}F\\ \end{align} }$

これでいいでしょうか？　華氏は摂氏0度のとき32度から始まっているので、上の式にその要素を加えてあげないといけませんでした。

${ \displaystyle \begin{align} F&=\frac{9}{5}C + 32\\ C&=\frac{5}{9}(F-32)\\ \end{align} }$

換算式が完成しました。 ※32度は華氏の値なので二つ目の式では括弧内で先に引く必要があります。

それではpythonで摂氏⇄華氏換算するコードを以下に書いて見ました。
※コメントアウトした部分はお好みで

def convertCtoF(degrees_c):
    f= 9/5*degrees_c + 32
    #f= str(f) + ' °F'
    return f

def convertFtoC(degrees_f):
    c= 5/9*(degrees_f -32)
    #c= round(c, 1)
    #c= str(c) + ' °C'
    return c

convertCtoF(14.5)

58.1

convertFtoC(80)

26.666666666666668

次に、摂氏と華氏換算グラフを作成してみます。

import numpy as np
list_c = np.linspace(0, 100, num=101)
list_f = convertCtoF(list_c)
list_tmp = np.vstack((list_c, list_f))

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

df_tmp = pd.DataFrame(list_tmp.T)

plt.figure(figsize=(4,6))
plt.plot(list_c, list_f)

list_idx_c = [20, 40, 60, 80]
plt.vlines(list_idx_c, ymin=min(list_f), ymax=[df_tmp.iloc[20, 1], df_tmp.iloc[40, 1], df_tmp.iloc[60, 1], df_tmp.iloc[80, 1]], linestyles='dashed', alpha=0.7)
plt.hlines([df_tmp.iloc[20, 1], df_tmp.iloc[40, 1], df_tmp.iloc[60, 1], df_tmp.iloc[80, 1]], xmin=0, xmax=list_idx_c, linestyles='dashed', alpha=0.7)

plt.xlim(min(list_c), max(list_c))
plt.ylim(min(list_f), max(list_f))
plt.xlabel('°C', fontsize=15)
plt.ylabel('°F', fontsize=15)

f:id:udonista:20171028041457j:plain

早見表的な感じでグラフを作成しましたが、分かりやすいのだろうか。。
以上になります。

勉強したことメモ

いろいろな勉強とつたわる日本語の練習をかねたブログです

幾何分布と負の二項分布

幾何分布 (geometric distribution)

負の二項分布 (negative binomial distribution)

二項分布(Binomial distribution)

考察

摂氏、華氏の換算