【数学溢れ話】【Token】二項分布(Binomial Distribution)の平均と分散と最尤値

まず「ある事象が毎回確率pで起こるとしてn回試行した時にその事象がk回起こる」確率密度関数の式形$P_n(k)$から出発しましょう。その導出には二項定理(Binomial Theorem)を用います。
[【Token】階乗と順列と組み合わせ](https://qiita.com/ochimusha01/items/90d9ef0a224bdd604721)

```math
(a+b)^n={}_nC_0a^nb^0+{}_nC_1a^{n-1}b^1+{}_nC_2a^{n-2}b^2+…

```

```math
+{}_nC_ra^{n-r}b^r+…+{}_nC_na^0b^n=\sum_{n=0}^n{}_nC_ra^{n-r}b^r
```
高校数学で暗記させられた以下の公式の大元ですね。

```math
(a+b)^2=a^2+2ab+b^2
```

```math
(a+b)^3=a^3+3a^2b+3ab^2+b^3
```

```math
(a+b)^4=a^4+4a^3b+6a^2b^2+4ab^3+b^4
```

各項の係数(組み合わせ上の出現数)は、いわゆる「パスカルの三角形」に従います。
[【Token】パスカルの三角形と二項定理または二項展開](https://qiita.com/ochimusha01/items/d668b19d5389b52e0dd2)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/388bfa4a-3fd5-6fff-004e-6f7bb944d3af.png)
これを以下の様に捉えるところから二項分布の概念は出発する訳です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/4a7ecb1c-7fae-a001-4bed-f33b76768a41.png)

```math
P_n(k)={}_nC_kP^kq^{n-k}(ただしq=1-p)
```
![スクリーンショット 2024-05-24 15.11.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/ba39195d-9d5e-e23a-f09c-de7bd8640c7b.png)

すなわち

```math
\sum_{k=0}^nP_n(k)=\sum_{k=0}^n{}_nC_kP^kq^{n-k}=(p+q)^n=1
```

* 自明の場合としてn=0の時は「1通りの結果しかない=いわゆる情報量0の状態」なのでn=1の場合から考えると…


```math
\sum_{k=0}^nP_1(k)=(p+q)^1={}_1C_0p^1q^0+{}_1C_1p^0q^1
```

```math
=1*p*1+1*1*q=p+q=1
```

```math
コイン(p=q=\frac{1}{2})で考えると\frac{1}{2}+\frac{1}{2}=1
```

* n=2の場合…


```math
\sum_{k=0}^nP_1(k)=(p+q)^2={}_2C_0p^2q^0+{}_2C_1p^1q^1+{}_2C_2p^0q^2
```

```math
=p^2+2pq+q^2=1
```

```math
コイン(p=q=\frac{1}{2})で考えると\frac{1}{4}+\frac{1}{2}+\frac{1}{4}=1
```

* n=3の場合…


```math
\sum_{k=0}^nP_1(k)=(p+q)^3={}_3C_0p^3q^0+{}_3C_1p^2q^1+{}_3C_2p^1q^2+{}_3C_3p^0q^3
```

```math
=p^3+3p^2q+3pq^2+q^3=1
```

```math
コイン(p=q=\frac{1}{2})で考えると\frac{1}{8}+\frac{3}{8}+\frac{3}{8}+\frac{1}{8}=1
```

* n=4の場合…


```math
\sum_{k=0}^nP_1(k)=(p+q)^4={}_4C_0p^4q^0+{}_4C_1p^3q^1+{}_4C_2p^2q^2+{}_4C_3p^1q^3+{}_4C_4p^0q^4
```

```math
=p^4+4p^3q+6p^2q^2+4pq^3+q^4=1
```

```math
コイン(p=q=\frac{1}{2})で考えると\frac{1}{16}+\frac{1}{4}+\frac{3}{8}+\frac{1}{4}+\frac{1}{16}=1
```

この様に、高校数学においては公式の暗記に終始した「二項式の展開」と「組み合わせ問題」が「パスカルの三角形」概念によって統合され、しかも確率論における中心極限定理に迫る勢いを見せる様子を目の当たりにすると不思議な感慨にとらわれます。
![20191002222245.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/7f909cc2-3d12-f148-6548-0215f815d0f4.gif)

同じ数理が幾何学や集合論の背景にも透けて見えるのが興味深いところですね。
[【Token】デカルト座標系概念の呪い？](https://qiita.com/ochimusha01/items/46b503c5679ba65270b5)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/9e64b4a9-145d-34a5-aad7-74df77c8be2a.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/0cb41bc2-fcae-7dcc-e2f5-f843f57be8ad.png)
もっと早く裏では全てがこうやって一つにつながってると気づいていたら、もっと早く数学に興味を持てていたのに…
[二項分布](https://jp.mathworks.com/help/stats/binomial-distribution.html)

* 確率質量関数(PMF=Probability Mass Function)

```math
f(x;n,p)= {}_nC_x p^x(1-p)^{n-x}
```

* 累積分布関数

```math
f(x;n,p)=\int_{x=0}^n {}_nC_x p^x(1-p)^{n-x}
```
とりあえず、基礎の基礎？

# 二項分布の平均

以下のサイトによれば「場合の数」は2009年度カリキュラムより小学校6年から教えている様です。
[親子中学数学：場合の数](https://oyako-cyugaku.com/category6/entry32.html)

「異なるn個の中から異なるr個を取り出して１列に並べた順列(Sequence Without Repetition)」${}_nP_r=\frac{n!}{(n-r)!}$や「異なるn個の中から選んだ異なるr個との組み合わせ(Combination)」${}_nC_r=\frac{{}_nP_r}{n!}$あるいは$\frac{n!}{n!(n-r)!}$自体は、それだけ初学者向けという事ですね。

```math
{}_3C_0=\frac{3!}{0!*(3-0)!}=\frac{3!}{0!3!}=1
```

```math
{}_3C_1=\frac{3!}{1!(3-1)!}=\frac{3!}{1!2!}=\frac{3*2*1}{1*2}=\frac{6}{2}=3
```

```math
{}_3C_2=\frac{3!}{2!(3-2)!}=\frac{3!}{2!1!}=\frac{3*2*1}{2*1}=\frac{6}{2}=3
```

```math
{}_3C_3=\frac{3!}{3!(3-3)!}=\frac{3!}{3!0!}=1
```

よって

```math
{}_3C_0p^3q^0+{}_3C_1p^2q^1+{}_3C_2p^1q^2+{}_3C_3p^0q^3
```

```math
=p^3+3p^2q+3pq^2+q^3
```

上掲の計算を見てもわかる通り、二項係数には左右対称性があります(そもそもpとqが可換な対称式)。
[二項係数の有名公式一覧と2つの証明方針](https://manabitimes.jp/math/588)

```math
{}_nC_r=\frac{n!}{r!(n-r)!}=\frac{n!}{(n-r)!(n-(n-r))!}={}_nC_{n-r}
```

だから直感的に期待値(平均)はnp(n=1の場合だけを抽出したベルヌーイ分布の場合はp)とわかりますが、実際の計算は以下の様になります。
[二項分布の平均と分散の二通りの証明](https://manabitimes.jp/math/913)
![スクリーンショット 2024-05-25 0.35.16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/c099a871-d0e6-fe6b-1fcc-557ed521b3e3.png)

### 期待値(平均)の線形性を使った証明
i回目に当たった時1、当たらない時0となる確率変数Xを$x_i=(x_1,x_2,…,x_n)$と定める。期待値の線形性$E[\sum_{i=1}^n x_i]=\sum_{i=1}^nE[x_i]$より

```math
E[X]=E[x_1]+E[x_2]+…+E[x_n]
```
右辺の各項はいずれも当選確率pと等しいので$E[X]=np$

### 二項分布の期待値計算による証明

期待値の定義$E[X]=\sum_{k=0}^n kP(X=k)$より

```math
\sum_{k=0}^n kP(X=k)=\sum_{k=1}^n k_nC_kP^kq^{n-k}
```
ここで二項係数の以下の公式を使う、

```math
k_nC_k=k\frac{n!}{k!(n-k)!}=\frac{n(n-1)!}{(k-1)!(n-k)!}={}_{n-1}C_{k-1}
```

代入して

```math
\sum_{k=1}^n {}_{n-1}C_{k-1}P^kq^{n-k}=np\sum_{k=1}^n {}_{n-1}C_{k-1}P^{k-1}q^{n-k-1}
```

```math
=np\sum_{k=1}^{n-1} {}_{n-1}C_{k}P^{k}q^{n-k}=np(p+q)^{n-1}
```
ここでp+q=1なのでnp

# 二項分布の分散

二項分布の分散np(1-p)=npq(n=1の場合だけを抽出したベルヌーイ分布の場合はp(1-p)=pq)の直感的イメージは以下。
![スクリーンショット 2024-05-25 0.35.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/3c08ed18-a58f-dedd-66f7-c5d6fbf910c3.png)
相加相乗平均の不等式$\frac{a+b}{2}≧\sqrt{ab}$の等号条件a=bの様に$p=q=\frac{1}{2}$の時の$\frac{1}{4}$が最大でpないしはqが０(もう一方が1)の時に最小0となります、
[相加相乗平均の不等式：意味・例題・おもしろい証明](https://manabitimes.jp/math/566)

* 相加相乗平均の不等式$\frac{a+b}{2}≧\sqrt{ab}$の等号条件a=bの様に…直接対応するのは、以下の確率質量関数で表現されるベルヌーイ分布の標準偏差$\sqrt{pq}=\sqrt{p(1-p)}$。

```math
f(k;p)=p^k(1-p)^{1-k}
```

実際の計算は以下の様になります。

### 分散の線形性を使った証明

期待値(平均)同様、その線形性を利用して

```math
V[X]=V[x_1]=V[x_2]+…+V[x_n]=nV[x_1]
```
そして$V[x_1]$について、公式$V[aX]=a^2V[X]$を用いて

```math
p(1-p)^2+(1-p)p^2=p(1-p)=pq
```

# 二項分布の分散の計算による証明

分散の定義$E[X^2]-E[X]^2$と、上で求めた 二項分布の期待値E[X]=npを用いて

```math
E[X^2]=\sum_{k=0}^n k^2P(X=k)
```

```math
=n(n-1)\sum_{k=2}^n {}_nC_k P^k(1-P)^{n-k}+\sum_{k=1}^n k {}_nC_k P^k(1-P)^{n-k}
```

```math
=n(n-1)\sum_{k=1}^n {}_{n-2}C_{k-2} P^k(1-P)^{n-k}+E[X]
```

```math
=n(n-1)p^2\sum_{k=1}^n {}_{n-2}C_{k-2} P^{k-2}(1-P)^{n-k}+np
```

```math
=n(n-1)p^2\sum_{k=1}^n {}_{n-2}C_{k} P^{k-2}(1-P)^{n-k-2}+np
```

```math
=n(n-1)p^2(p+q)^{n-2}+np=n(n-1)p^2+np
```
最初の分散の定義に立ち返って代入

```math
E[X^2]-E[X]^2=n(n-1)p^2+np-(np)^2
```

```math
=np((n-1)p+1)-np(np)=np(np-p+1-np)
```

```math
=np(1-p)=npq
```

こうして求めた平均と分散は、サンプル数Nが十分大きく、正規分布による矜持が可能な条件下において、その二項分布を近似する正規分布にパラメーターとして引き渡されるのです。
[二項分布の正規近似条件について](https://data-everyday.com/stat_ml/post-248/)

# 二項分布の最尤推定

小学校から習う「場合の数」や、高校から習う「反復試行」問題においては、単純化の為、以下の二項分布の確率質量関数の確率pが最初から分かっているのがセオリーとなっている様です。
[反復試行と期待値（二項分布）](https://examist.jp/mathematics/probability/hanpukusikou-kitaiti/)

```math
f(x;p)={}_nC_x p^x(1-p)^{n-x}
```

そりゃ確かに「コインを1枚投げるとする(表と裏が出る確率がそれぞれ$\frac{1}{2}$)」「サイコロを1個振るとする(1から6までの出目が出る確率がそれぞれ$\frac{1}{6}$)」という前提から出発しても誰も疑問なんて持ちませんから、説明の手間が大いに省けます。これが「子供用プール」の発想で、「大人用海水浴場」においては直接観察された結果以外は何も信じられないのです。この目的の為に研ぎ澄まされた方法論を最尤推定(Maximum Likelihood Estimation)あるいは最尤法(Method of maximum likelihood)といいます。

```math
L(θ;x)={}_nC_x θ^x(1-θ)^{n-x}
```

尤度関数導入にあたって記号を「(コインの表裏やサイコロの出目の様に)あらかじめ自明の場合として分かっている確率p」から「未知の分布の未知のパラメーターθ」に置き換えました。
[最尤法によるパラメータ推定の基礎を理解する(二項分布と正規分布のパラメータの最尤推定量の導出)](https://qiita.com/g-k/items/4b2d557c24558074f07f)

> あるパラメータθに従う確率密度関数をf(x;θ)とした時、尤度関数l(θ;x)=f(x;θ)が最大となるような推定量$θ=\hat{θ}$を最尤推定量と呼びます。

* とりあえず未知の分布を二項分布X~B(n,p)と推定する。
* すると自明の場合として確率pが推定対象パラメーターとして選択される。

この時点で既にこれだけの判断が下されてる訳ですね。これが連続一様分布関数$\int_a^b \frac{1}{b-a}$の場合だと、パラメーターはルベール測度b-aとなり、その最尤推定は「(単純化の為に$a=0,b=x_{max}$と置いた上で)$x_{max}$を求める」内容となる訳です。
[【Token】連続一様分布の平均と分散と最尤値](https://qiita.com/ochimusha01/items/ca01c8671522ede6bdda)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/e01f31d3-4cc3-a803-2d5a-c56f7630b602.png)

* 、「(直接観察された結果以外は何も信じられない)大人用海水浴場の世界」…例えば連続一様分布の最も単純な形態たる$f(x)=\left[1\right]_0^1$の正体はパラメーターがα=β=1の特殊なベータ分布X~B(x;α,β)だったりする。
[【徹底解説】ベータ分布と一様分布](https://academ-aid.com/statistics/beta-uniform)
* ベータ分布X~B(x;α,β)を用いれば二項分布の累積分布関数も連続的に表現可能となる。
[Wikipedia「不完全ベータ関数」](https://ja.wikipedia.org/wiki/%E4%B8%8D%E5%AE%8C%E5%85%A8%E3%83%99%E3%83%BC%E3%82%BF%E9%96%A2%E6%95%B0)

* 「半径r⇄直径2r⇄円周2πr⇄球の表面積$4πr^2$⇄球の体積$\frac{4}{3}πr^3$⇄円の面積$πr^2$」の推移を完全記述する為には

そう、数理の深淵は常にあなたの足元にごく当たり前の様に隠れ潜んでいて、あなたに発見されて結縁に至る機会を虎視眈々と狙っているという訳です…


### 実際の計算過程

最尤推定の計算は非常に煩雑な内容となるので簡略化の為、対数をとって偏微分するのがセオリーです。理由の詳細についてあちこちに様々な記載がありますが、とりあえずChatGPTにまとめさせた「優等生的回答」が以下。

>尤度関数を対数変換する理由は、対数変換することで計算が簡単になり、数値的な安定性が向上するためです。具体的には以下のような利点があります：
>
>計算の簡略化:
>尤度関数は通常、多くの確率密度関数や確率質量関数の積として表現されます。特に観測データが大量にある場合、積の形は非常に複雑になります。しかし、対数を取ることで積は和に変換され、計算が大幅に簡略化されます。
>
>数値的な安定性:
>尤度関数の値は非常に小さくなることが多いため、直接計算するとアンダーフローの問題が発生する可能性があります。対数変換を行うことで、非常に小さな値の積を計算する際の数値的な問題を回避できます。対数を取ることで、非常に小さな値でも適切に扱うことができます。
>
>微分計算の容易さ:
最尤推定量を求めるためには、尤度関数をパラメータで微分し、その導関数を0に設定する必要があります。対数変換を行うことで、積の微分ではなく和の微分を計算すればよくなり、計算が簡単になります。
>
>解析的な扱いやすさ:
尤度関数自体は非線形な形状をとることが多いため、その解析が困難な場合があります。対数を取ることで、関数の形状が単純化され、解析的に扱いやすくなります。
>
>これらの理由により、尤度関数を対数変換して対数尤度関数を用いることが一般的です。対数尤度関数を最大化することは、元の尤度関数を最大化することと等価であり、計算の効率性と安定性が向上します。

それでは実際に取り組んでみましょう。尤度関数L(θ;x)におけるθの最大値を求める為に微分する訳ですが、上掲の理由により尤度関数L(θ;x)そのものでなく尤度関数L(θ;x)の対数を取った対数尤度関数l(θ;x)=LogL(θ;x)を使います。

```math
l(θ;x)=LogL(θ;x)=\log(\sum_{x=0}^n f(x;n,p))
```

```math
=\log({}_nC_xp^x(1-p)^{n-x})
```

```math
=\log{}_nC_x+x\log p+(n-x)\log(1-p)
```
ただし二項係数項$\log{}_nC_x$はパラメーターを含まない=最尤推定に関係してこない定数項なので除去(つまり最尤推定は二項分布と、そのn=1の場合たるベルヌーイ分布を区別しない？)

```math
=x\log p+(n-x)\log(1-p)
```
これをθで微分します。

```math
\frac{l(θ)}{∂θ}=\frac{l(θ)}{∂θ}x\log p+\frac{l(θ)}{∂θ}(n-x)\log(1-p)
```

```math
=\frac{x}{θ}+\frac{n-x}{1-θ}
```
ここで$\frac{l(θ)}{∂θ}=0$と置いて解くと最尤推定値$θ=\frac{x}{n}$に到達。

```math
\frac{x}{θ}+\frac{n-x}{1-θ}=0
```

```math
\hat{θ}=\frac{x}{n}
```

辿り着いてみればなんのこともない…「成功した試行回数を試行回数全体で割る($frac{x}{n}$)」のが二項分布における最尤推定で、「コインを連続3回投げて(３枚同時に投げて)偶然3枚とも表が出た場合、このコインには表しかない(p=1,q=0)と思い込むのを防ぐ特別な秘策」など別に存在せず、ただその状態が数学的により厳密に表現可能となっただけという次第。
[最尤推定量とは？初めての人にも分かりやすく解説]()

そんな感じで以下続報…