Scipyで標本分布を計算する(離散分布)

目的

母集団の確率分布が与えられたとき，標本分布を計算します。
今回は離散分布を扱います。

方法

方針

• pythonの科学技術計算パッケージ集であるscipyの，統計パッケージscipy.statsと，数値計算ライブラリのnumpyを使います。
• 次の3ステップで求めます。
  1. 母集団の離散確率分布をscipy.statsの上で離散分布として定義して
  2. サンプリングを行い
  3. 得られたサンプルから標本分布を計算します

離散分布の定義

scipy.stats.rv_discrete
以下の内容はscipyのドキュメンテーションからの抜粋です。
次のようにして読めるヘルプのExampleのところに載っています。

from scipy import stats
help(stats.rv_discrete)

離散分布の定義

$xk$に離散確率変数を，$pk$に対応する確率を入れて，stats.rv_discreteに渡すだけです。
7個の離散値をとる確率変数の場合の例になっています。

from scipy import stats
xk = np.arange(7)
pk = (0.1, 0.2, 0.3, 0.1, 0.1, 0.0, 0.2)
custm = stats.rv_discrete(name='custm', values=(xk, pk))

離散分布の可視化

離散分布の確率質量関数(pmf)を可視化するには次のようにすれば良いです。

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)
ax.plot(xk, custm.pmf(xk), 'ro', ms=12, mec='r')
ax.vlines(xk, 0, custm.pmf(xk), colors='r', lw=4)
plt.show()

サンプリング

作成した離散分布からサンプリングを行います。
たとえば1000回ランダムサンプリングする場合は次ようにすればよいです。
(rvs: random variables)
再現性を確保するために乱数のシードを固定したい場合は，一行目のコメントを外してください。

import numpy as np
# np.random.seed(0)
sampled_rvs = custm.rvs(size=1000)

標本分布の計算

numpy のヒストグラム関数を利用して計算します。
binのところを+1することを忘れないように注意。

f = np.histogram(sampled_rvs, bins = np.arange(7 + 1), density=True)[0]

まとめて関数にする

以上の内容をまとめます。
先ほどまでは，一次元の確率分布の場合を扱っていましたが，多次元の同時確率分布にも対応できるように，
最初にp.ravel()で一次元に潰してから標本分布を計算し，reshapeしてから返します。
引数のpはnupmyの配列で渡してください。

import numpy as np
from scipy import stats

def sampling_dist(p, size, seed = 0):
    xk = np.arange(p.size)
    pk = p.ravel()
    true_dist = stats.rv_discrete(name='true_dist', values=(xk, pk))

    np.random.seed(seed)
    sampled_rvs = true_dist.rvs(size=size)

    f = np.histogram(sampled_rvs, bins = np.arange(p.size + 1), density=True)[0]
    f.reshape(p.shape)

    return(f.reshape(p.shape))

p = np.array([[0.10,0.10],[0.10,0.15],[0.15,0.10],[0.10,0.20]])
p_est = sampling_dist(p, 1000)
print(p)
print(p_est)