pythonで超幾何分布の計算をする

python で超幾何分布の計算を行う方法を調べてみました。
結果を記載します。

超幾何分布とは何か

• 一言で言うと
  超幾何分布 | 統計用語集 | 統計WEB から引用します。

  N個からなる集団が2つの性質をもつときに、各集団がM個、N－M個であるとする。ここからn個を取り出すときに、M個の中からm個が取り出されるとすると、この確率分布は超幾何分布に従う。

• 公式
  母集団の要素の数N，ある属性の要素数M，標本数のn で以下のように記載できます。

$$ f(m) = \frac{{}M\mathrm{C}_m\cdot{} {{}_N\mathrm{C}_n} $$}\mathrm{C}_{n - m}

• 超幾何分布の実用例
  捕獲再捕獲法に用いられます。以下、sec2 からの引用です。

  超幾何分布は，生態学で動物の生息数を推定するのに用いられる．m 匹の動物を捕獲して目印をつけたのち放す．十分時間が経過して，捕獲した動物とそうでない動物>が混ざり合ったと考えられてから，k 匹を再捕獲する．再捕獲で捕まった目印つきの個体数を X とすると，対象となっている動物の生息数を m + n と考えると，X は超幾何>分布に従う．これを捕獲再捕獲法（capture - recapture method），もしくは，標識・再捕獲法（mark - recapture method）という．

超幾何分布の計算ができるライブラリについて

他の分布と同様に、numpy と scipy で計算ができます。

• numpy.random.hypergeometric

• scipy.stats.hypergeom

例によって、scipy.stats.hypergeom のほうがいろいろできるので、どちらかしか使わないなら、scipy.stats.hypergeom の使い方を覚えるのがよいと思います。

numpy.random.hypergeometric を使う

超幾何分布に従う、配列を取得することができます。
パラメータは、ngood 属性Aの総数、nbad 属性Bの総数、nsamp 標本数、 size 試行回数 を指定すると、試行回数を要素数とする 属性Aの取り出された回数が返ります。

• ヒストグラム描画

# ヒストグラム描画 
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# ngood 白玉 90個、 nbad 赤玉 10、 標本数 50     
ngood, nbad, nsamp = 90, 10, 50
# 10000回試行する。   
# 標本数 50 個のうち、白い玉が何個あったかが返る
s = np.random.hypergeometric(ngood, nbad, nsamp, 10000)
# ヒストグラムに描画
plt.hist(s)

(array([    7.,    84.,   388.,  1116.,  2090.,  2524.,  2252.,  1097.,
          382.,    60.]),
 array([ 40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,  49.,  50.]),
 <a list of 10 Patch objects>)

• n回起こる確率を求める
  n回 属性Aが起こる確率を求めます。

import numpy as np
# ngood 白玉 90個、 nbad 赤玉 10、 標本数 50     
ngood, nbad, nsamp = 90, 10, 50

# 10000回試行する。   
# 標本数 50 個のうち、白い玉が何個あったかが返る
s = np.random.hypergeometric(ngood, nbad, nsamp, 10000)
print("白玉が取り出される回数45回。 赤が5つ。", (s == 45).sum() / 10000)
print("白玉が取り出される回数40回。 赤が0つ。", (s == 40).sum() / 10000)

白玉が取り出される回数45回。 赤が5つ。 0.2632
白玉が取り出される回数40回。 赤が0つ。 0.0007

scipy.stats.hypergeom を使う

scipy.stats.hypergeom を使って、負の2項分布のグラフを描画してみます。　　　

• ヒストグラムを描画する

%matplotlib inline
from scipy.stats import hypergeom
import matplotlib.pyplot as plt

# 100個 玉があって、そのうち 10個 が赤玉、 50　個 を取り出した時に、赤玉が出る確率 をグラフ出力する。
[M, n, N] = [100, 10, 50]
rv = hypergeom(M, n, N)
x = np.arange(0, n+1)
pmf_dogs = rv.pmf(x)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, pmf_dogs, 'bo')
ax.vlines(x, 0, pmf_dogs, lw=2)
ax.set_xlabel('# of dogs in our group of chosen animals')
ax.set_ylabel('hypergeom PMF')
plt.show()

• n回起こる確率を求める
  確率質量関数 probability mass function で 5回赤玉が出た時の確率を求めることができます。
  超幾何分布は離散型の離散型の確率変数なので、 確率質量関数で、 連続確率変数 の場合は、確率密度関数 probability density function になり、scipy では、分布ごとに、意識していて、当たり前かもですが、メソッド名がpmf、pdf で違います。

# 100個 玉があって、そのうち 10個 が赤玉、 50　個 を取り出す時
[M, n, N] = [100, 10, 50]
# 5回赤玉が出る確率を求める
hypergeom.pmf(5, M, n, N)

0.25933354622553628

• ランダムな配列を返す
  rvs を使うと、numpy.random.hypergeometric と同様の超幾何分布に従う、ランダムな配列を取得できます。

%matplotlib inline
import matplotlib.pyplot as plt
# 100個 玉があって、そのうち 10個 が赤玉、 50　個 を取り出す時
[M, n, N] = [100, 10, 50]
R = hypergeom.rvs(M, n, N, size=10000)
plt.hist(R)

(array([    4.,    69.,   395.,  1180.,  2157.,  2577.,  2041.,  1119.,
          388.,    70.]),
 array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.]),
 <a list of 10 Patch objects>)

引数の考え方が反対? なので指定方法を変更します。

%matplotlib inline
import matplotlib.pyplot as plt
# 100個 玉があって、そのうち 90個 が白玉、 50　個 を取り出す時
[M, n, N] = [100, 90, 50]
R = hypergeom.rvs(M, n, N, size=10000)
plt.hist(R)

(array([    5.,    62.,   373.,  1179.,  2091.,  2614.,  2090.,  1131.,
          360.,    95.]),
 array([ 40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,  49.,  50.]),
 <a list of 10 Patch objects>)

参考

以下、参考になりました。

• sec2

• 超幾何分布 | 統計用語集 | 統計WEB

• 確率質量関数 - Wikipedia

• numpy.random.hypergeometric — NumPy v1.14 Manual

• scipy.stats.hypergeom — SciPy v1.1.0 Reference Guide

以上です。