多次元確率分布：多項分布

$${0}$$か$${1}$$かの二項分布の拡張であり、$${d}$$面を持つサイコロの各面が出る確率が$${{\bf P}=(p_1, \cdots, p_d)^T}$$とし、それぞれが出た回数$${{\bf x}=(x_1, \cdots, x_d)^T}$$が従う確率分布である。
$${\displaystyle{\sum_{i=1}^d p_i=1, p_i \geq 0, i \forall d}}$$
$${\displaystyle{\sum_{i=1}^d x_i=n, x_i \geq 0, i \forall d}}$$
$${j}$$面が$${x_j}$$回出る確率は、$${p_j^{x_j}}$$で、$${n}$$回で各面が$${(x_1, \cdots, x_d)}$$回出る組み合わせは、
$${\displaystyle{\frac{n!}{x_1!x_2!\cdots x_d!}}}$$から、$${{\bf x}}$$の確率質量関数は、
$${f({\bf x})=\displaystyle{\frac{n!}{x_1!x_2!\cdots x_d!}p_1^{x_1}\cdots p_d^{x_d}}}$$
で与えられる。
$${d=2}$$の時は、
$${f({\bf x})=\displaystyle{\frac{n!}{x_1!x_2!}p_1^{x_1}p_2^{x_2} = \frac{n!}{x_1!(1-x_1)!}p_1^{x_1}(1-p_1)^{n-x_1} }}$$
で二項分布の確率質量関数となる。

多項定理
$${(p_1+\cdots+p_d)^n=\displaystyle{\sum_{x\in{\bf x}}\frac{n!}{x_1!x_2!\cdots x_d!}p_1^{x_1}\cdots p_d^{x_d}}}$$
から、モーメント母関数は、
$${M_{{\bf x}}({\bf t})=E[e^{{\bf t}\cdot{\bf x}}]}$$
$${=\displaystyle{\sum_{x\in{\bf x}}e^{t_1x_1}e^{t_2x_2}\cdots e^{t_dx_d}\frac{n!}{x_1!x_2!\cdots x_d!}p_1^{x_1}\cdots p_d^{x_d}=(p_1e^{t_1}+\cdots+p_de^{t_d})^n}}$$
となる。
これを用いて、$${x_j}$$の期待値
$${E[x_j]=\displaystyle{\frac{\partial M_{{\bf x}}({\bf t})}{\partial t_j}\Big|_{{\bf t}={\bf 0}}= p_je^{t_j}n(p_1e^{t_1}+\cdots+p_de^{t_d})^{n-1}\Big|_{{\bf t={\bf 0}}}=np_j}}$$
が得られる。
共分散は、２階微分を以下のように計算して得られる。
$${\displaystyle{\frac{\partial^2 M_{{\bf x}}({\bf t})}{\partial t_i \partial t_j}=\frac{\partial }{\partial t_i }p_je^{t_j}n(p_1e^{t_1}+\cdots+p_de^{t_d})^{n-1} }}$$
$${\displaystyle{=p_i p_j e^{t_i+t_j}n(n-1)(p_1e^{t_1}+\cdots+p_de^{t_d})^{n-2} }}$$
$${\displaystyle{ + \delta_{ij}p_je^{t_j}n(p_1e^{t_1}+\cdots+p_de^{t_d})^{n-1}}}$$
ここで、$${\delta_{ij}}$$は
$${\delta_{ij} = \left\{\begin{array}{ll}1 & (i=j)\\0 & (i \neq j)\end{array}\right.}$$
のデルタ関数である。
よって、
$${\displaystyle{\frac{\partial^2 M_{{\bf x}}({\bf t})}{\partial t_i \partial t_j} \Big|_{{\bf t}={\bf 0}}=\delta_{ij}p_jn + p_ip_jn(n-1) }}$$
これから、共分散
$${Cov[x^j, x^i]=\delta_{ij}p_jn + p_ip_jn(n-1) - n^2p_ip_j = np_j(\delta_{ij}-p_i)}$$
が得られる。
多項分布に従う乱数のpythonでの実装は、numpy.randomのmultinominalで行うことができる。
ここでは、$${{\bf P}=(\frac{1}{2},\frac{1}{3},\frac{1}{6})^T}$$として、$${{\bf x}}$$を20個生成した。

import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)

pvec= (1/2,1/3,1/6)
x=rand.multinomial(n=20, pvals=pvec, size=20).transpose()
plt.barh(range(20), x[0])
plt.barh(range(20), x[1], left=x[0], color='g')
plt.barh(range(20), x[2], left=x[0]+x[1], color='r')
plt.title("Multinominal sum of each x vector");