多項分布の最尤推定量を求める（ラグランジュの未定乗数法を使わない）

前置き

「多項分布　最尤推定量」などで検索するとラグランジュの未定乗数法を使っている解説が多い印象だった.

それらに問題があるわけではまったくないが色々な回答例があるほうが, ないよりはいいと思ったので書いた.

"N=1" のとき

$x = (x_1, \ldots, x_m)'$ をサイズパラメータ $n$, 確率パラメータ $p = (p_1, \ldots, p_m)'$ の多項分布に従う確率変数とする. このことを次のようにも書く.

$$x \sim \mathrm{Mult}(x | n, p).$$

$p_j$ は確率なので0以上で $\sum_j p_j = 1$ である.

多項分布の確率関数は

$$\mathrm{Mult}(x | n, p)= C_{n,x} \prod_{i=1}^m p_i^{x_i}$$

である. ここで $C_{n,x}$ は多項係数を表すが, このあと対数とって微分すると消える（0になる）部分なので省略した書き方にしている.

$n = \sum_j x_j$ は固定で既知. 未知パラメータの $p$ を推定したい.

$p$ についての対数尤度は

$$l(p)= \log C_{n,x} + \sum_{j=1}^m x_j \log p_j$$

これを偏微分して0となる点を求めたい.

$$p_m = 1-\sum_{j=1}^{m-1} p_j$$

に注意して,

$$\frac{\partial}{\partial p_j} l(p)= x_i / p_j - x_m/p_m.$$

右辺が0になる $p_j$ は

$$x_jp_m - p_jx_m = 0$$

を満たすので,

$$p_j = x_jp_m/x_m \tag{1}$$

を満たす.

よって,

$$\sum_{j=1}^m p_j = (p_m/x_m)\sum_{j=1}^m x_j.$$

左辺は1でなければいけないので $p_m$ の最尤推定量 $\hat{p}_m$は

$$\hat{p}_m = \frac{x_m}{\sum_{j=1}^m x_j}.$$

式(1)で $p_m = \hat{p}_m$ とすると, $p_j$ の最尤推定量 $\hat{p}_j$ は

$$\hat p_j = \frac{x_j}{\sum_{j=1}^m x_j}.$$

"N>1" のとき

独立に同じ確率パラメータ $p$ を持つ多項分布に従う確率変数 $\boldsymbol{x}_i$ ($i=1, \ldots, N$) があるとき, その同時分布の確率関数は

$$\prod_{i=1}^N \mathrm{Mult}(\boldsymbol{x}_i | n_i, p)=\prod_{i=1}^N C_{n_i,x_i} \prod_{j=1}^m p_j^{x_{ij}}$$

なので $p$ についての対数尤度は

$$l(p)= \left(\sum_{i=1}^N \log C_{n_i,\boldsymbol{x}_i} \right) + \left(\sum_{j=1}^m \left\{\sum_{i=1}^N x_{ij} \right\} \log p_j \right).$$

$\sum_{i=1}^N x_{ij}$ を新たな $x_j$ と見れば後は「$N=1$のとき」と一緒である.

Rによる実装例

p <- c(0.1,0.2,0.3,0.4) #set p
set.seed(132)
X <- rmultinom(1, 100, p) #sample (N=1 n=100)
phat <- sweep(X, 1, colSums(X), FUN = "/") #equivalent to X/sum(X)
cat("true p: ", p, "\n")
cat("estimates: ", phat, "\n")

> cat("true p: ", p, "\n")
true p:  0.1 0.2 0.3 0.4 
> cat("estimates: ", phat, "\n")
estimates:  0.11 0.26 0.25 0.38