超幾何分布の期待値、分散、モーメント母関数
非復元抽出の回数を重ねるごとに確率分布が変化する確率変数の期待値と分散、モーメント母関数を求める。
Aが$${M}$$個、Bが$${N-M}$$個入っている箱から、$${n}$$個を取り出した中にAが$${x}$$個入っている確率の分布$${HG(N,M,n)}$$は以下のように求められる。
$${N}$$から$${n}$$個取り出す組み合わせは$${{}_NC_n}$$、$${M}$$個から$${x}$$個取り出される組み合わせは$${{}_MC_x}$$、同様に、$${N-M}$$個から$${n-x}$$個取り出される組み合わせは$${{}_{N-M}C_{n-x}}$$、よって、$${f(x)=\displaystyle{\frac{{}_MC_x{}_{N-M}C_{n-x}}{{}_NC_n}}}$$、
ここで、$${x\in\{\max\{0,n-(N-M)\},\min\{n,M\}\}}$$であることに注意する。
HG(N,M,n)のモーメント母関数は、
$${M_x(t)=E[e^{tx}]=\displaystyle{\frac{1}{{}_NC_n}\sum_x e^{tx}{}_MC_x{}_{N-M}C_{n-x} }}$$と与えられ、これは超幾何級数であることから、これを超幾何分布と呼ぶ。
期待値と分散は、モーメント母関数を使わないで求める。
$${E[x]=\displaystyle{\frac{1}{{}_NC_n} \sum^{n}_{x=0}x{}_MC_x{}_{N-M}C_{n-x}=\frac{1}{{}_NC_n} \sum^{n}_{x=1}x{}_MC_x{}_{N-M}C_{n-x}}}$$
ここで、$${\displaystyle{{}_MC_x=\frac{M}{x}{}_{M-1}C_{x-1}}}$$を使って、
$${\displaystyle{\frac{1}{{}_NC_n} \sum^{n}_{x=1}x{}_MC_x{}_{N-M}C_{n-x}=\frac{M}{{}_NC_n}\sum^{n}_{x=1}{}_{M-1}C_{x-1}{}_{N-M}C_{n-x}}}$$
$${\displaystyle{= \frac{M}{{}_NC_n} \sum^{n}_{x=0}{}_{M-1}C_{x}{}_{N-M}C_{n-x-1}}}$$
同様に、$${\displaystyle{{}_NC_n=\frac{N}{n}{}_{N-1}C_{n-1}}}$$から、
$${\displaystyle{\frac{M}{{}_NC_n} \sum^{n}_{x=0}{}_{M-1}C_{x}{}_{N-M}C_{n-x-1}=\frac{nM}{N}\frac{1}{{}_{N-1}C_{n-1}}\sum^{n}_{x=0}{}_{M-1}C_{x}{}_{N-M}C_{n-x-1} }}$$
$${\sum_xf(x)=1}$$から、$${{}_NC_n=\displaystyle{\sum_{x=0}^{n}{}_MC_x{}_{N-M}C_{n-x}}}$$
よって、
$${{}_{N-1}C_{n-1}=\displaystyle{\sum_{x=0}^{n-1}{}_MC_x{}_{N-1-M}C_{n-1-x}=\sum_{x=0}^{n-1}{}_{M-1}C_x{}_{N-M}C_{n-1-x}}}$$（最後に$${M+1}$$を$${M}$$に変換している）。
これらから、
$${E[x]=\displaystyle{\frac{nM}{N}\frac{1}{{}_{N-1}C_{n-1}}\sum^{n}_{x=0}{}_{M-1}C_{x}{}_{N-M}C_{n-x-1}=\frac{nM}{N}}}$$
分散は、$${V[x]=E[x^2]-E^2[x]=E[x(x-1)]+E[x]-E^2[x]}$$を用いて得る方が簡単で、
$${E[x(x-1)]=\displaystyle{\frac{1}{{}_NC_n} \sum^{n}_{x=0}x(x-1){}_MC_x{}_{N-M}C_{n-x}=\frac{1}{{}_NC_n} \sum^{n}_{x=2}x(x-1){}_MC_x{}_{N-M}C_{n-x}}}$$
$${\displaystyle{= \frac{M}{{}_NC_n} \sum^{n}_{x=2}(x-1){}_{M-1}C_{x-1}{}_{N-M}C_{n-x}= \frac{M(M-1)}{{}_NC_n} \sum^{n}_{x=2}{}_{M-2}C_{x-2}{}_{N-M}C_{n-x}}}$$
$${\displaystyle{=\frac{M(M-1)}{{}_NC_n} \sum^{n-2}_{x=0}{}_{M-2}C_{x}{}_{N-M}C_{n-x-2} }}$$
$${\displaystyle{{}_NC_n=\frac{N(N-1)}{n(n-1)}{}_{N-2}C_{n-2}}}$$から、
$${E[x(x-1)]=\displaystyle{\frac{M(M-1)n(n-1)}{N(N-1)}\frac{1}{{}_{N-2}C_{n-2}} \sum^{n-2}_{x=0}{}_{M-2}C_{x}{}_{N-M}C_{n-x-2} }}$$
$${\displaystyle{=\frac{M(M-1)n(n-1)}{N(N-1)} }}$$
と、得られる。
超幾何分布に従う乱数の発生は、pythonのnumpy.randomのhypergeometricで実装できる。
ここでは、全体個数$${N=100}$$中にに含まれているAの個数$${M=2}$$とし、これを非復元抽出で$${n=10}$$個取り出す試行を3000回行った時の、各試行における$${n}$$中の$${A}$$の個数$${x}$$が従う乱数$${s}$$である。
numpy.random.hypergeometricでは、A以外が含まれている個数を返すことから、$${n-s}$$としている。
import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)
N=100
M=2
n=10
hygeo = rand.hypergeometric(N-M, M, n, 1000)
plt.title(f"Hylergeometric N={N},M={M},n={n}")
plt.hist(n-hygeo, 'auto')
plt.ylabel('Probability')
plt.xlabel('Data');
明らかに、抽出サンプルの中に入っているAの個数は少ない。