Traffine I/O

日本語

2022-12-09

多項分布

多項分布とは

多項分布とは、1回の試行で K 通りの事象 X_1, X_2, ..., X_K がそれぞれ確率 p_1, p_2, ..., p_K で得られるとして、試行を N 回繰り返たときにそれぞれの事象が起こる回数 X が従う確率分布です。

多項分布は二項分布K 次元に拡張した確率分布となっています。二項分布では事象は2つ(K=2)ですが、サイコロの目の数のような事象の数が6個(K=6)、つまり多次元の場合は多項分布になります。

N=1 のときの多項分布はカテゴリカル分布になります。また、N=1 かつ K=2 のときの多項分布はベルヌーイ分布になります。

多項分布の確率は次の式で表されます。

P(X=x;N,p_1, p_2, ...,p_K) = \prod_{k=1}^K \frac{N!}{n_k!}p^{x_k}_k
x_k \in \{0,1, ...,N\}, \quad \sum_{k=1}^{K} x_k=N

多項分布は Multinomial(N,p) と表記されることもあります。

多項分布の期待値と分散

カテゴリカル分布の期待値、分散はそれぞれ以下になります。

E(X_k)=Np_k \quad (k=1,2,...,K)
V(X_k)=Np_k(1-p_k) \quad (k=1,2,...,K)

N の影響を確認

K=3、事象の確率をそれぞれ0.2、0.3、0.5として、Nを変化させて分布への影響を確認します。実装コードは以下になります。

import numpy as np
from scipy.stats import multinomial
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from matplotlib.animation import FuncAnimation

rc('animation', html='html5')
np.random.seed(5)

# Set up formatting for the movie files
Writer = animation.writers['ffmpeg']
writer = Writer(fps=15, metadata=dict(artist='Me'), bitrate=1800)

p_v = np.array([0.2, 0.3, 0.5])
N_max = 20
p_max = 0.1

plt.style.use('ggplot')
cm = plt.get_cmap('jet')
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(projection='3d')

def update(N):
    # initialize previous frame
    plt.cla()

    N += 1

    x_vals = np.arange(N + 1)
    X1, X2 = np.meshgrid(x_vals, x_vals)

    x1_vals = np.delete(X1.flatten(), obj=(X1 + X2).flatten() > N)
    x2_vals = np.delete(X2.flatten(), obj=(X1 + X2).flatten() > N)
    x3_vals = np.where(x1_vals+x2_vals <= N, N - (x1_vals+x2_vals), 0.0)

    x_points = np.stack([x1_vals, x2_vals, x3_vals], axis=1)

    probability = multinomial.pmf(x=x_points, n=N, p=p_v)

    ax.bar3d(x=x1_vals - 0.45,
             y=x2_vals - 0.45,
             z=np.zeros_like(x1_vals),
             dx=0.9,
             dy=0.9,
             dz=probability,
             color=cm(probability / p_max),
             alpha=0.5,
             shade=True
    )
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('Probability')
    ax.set_title('$p=(' + ', '.join([str(phi) for phi in p_v]) + ')' +
                 ', N=' + str(N) + '$', loc='left')
    ax.set_zlim(0.0, p_max)

anime_prob = FuncAnimation(fig, update, frames=N_max, interval=100)

anime_prob.save('multinomial_dist.gif', writer='pillow', fps=1)

Multinomial distribution

N が大きくなるにつれて x がとり得る値の範囲が広がるため、分布のピークが小さくなっていることが分かります。

Ryusei Kakujo

researchgatelinkedingithub

Focusing on data science for mobility

Bench Press 100kg!