Daniel Bernoulliの発想を継承した二項分布の正規近似の導出 -

Daniel Bernoulliの正規曲線の方程式の導出法
https://doi.org/10.34336/jhsj.25.158_89

に、Daniel Bernoulliが1770年頃に発表したらしい、二項分布の正規近似の導出法が書かれてて、ちょっと面白い議論だった。

二項分布の正規近似は、de Moivreの1733年の論文に起源があるとされる。de Moivreは、本質的には、現代の殆どの教科書で使われてるのと同様、Stirlingの公式を使った議論を行ったのだろう。このStirlingの公式も、発見したのは、de Moivreとされている。

Stirlingの公式に関するde Moivreの議論は、以下の論説に詳しい。

The Early History of the Factorial Function
https://doi.org/10.1007/BF00389433

de Moivreの元々の動機からして、二項係数に端を発していて、その途中で、de Moivreは、
$\log (m-1)! = (m-\dfrac{1}{2})\log m - m + C + \dfrac{B_2}{1 \cdot 2} \dfrac{1}{m} + \dfrac{B_4}{3 \cdot 4} \dfrac{1}{m^3} + \dfrac{B_6}{5 \cdot 6} \dfrac{1}{m^5} + \cdots$
$C = 1 - \dfrac{B_2}{1\cdot2} - \dfrac{B_4}{3\cdot 4} - \dfrac{B_6}{5 \cdot 6} - \cdots$
という形の式を証明したらしい。 $B_n$ はベルヌーイ数。この定数 $C$ が $\log(2 \pi)/2$ であることを指摘したのが、Stirlingだったそうだ。この公式が、Stirlingの公式と呼ばれるようになった経緯は知らない。

この論説の中にあるDaniel Bernoulliから、Goldbachへの1729年の手紙というのは

File:DanielBernoulliLetterToGoldbach-1729-10-06.jpg
https://commons.wikimedia.org/wiki/File%3ADanielBernoulliLetterToGoldbach-1729-10-06.jpg

で読める。簡潔で証明もないけど、ガンマ関数の無限積表示が与えられているっぽい。これはGoldbachが出した問題に対する解答っぽいけど、同じ問題に対するEulerの解答が、
De progressionibus transcendentibus seu quarum termini generales algebraice dari nequeunt
https://scholarlycommons.pacific.edu/euler-works/19/
にある(Euler archiveのE19)。この報告では、ガンマ関数の積分表示(9ページ目。式の形は、現在の標準的なものとは違い、ちょっとした変数変換がいる。あと、英訳版は、lnが誤植っぽい？)が見られる。何にせよ、彼らは、確率論とは異なる文脈で、ガンマ関数に到達したらしい。

正規近似に関するD.Bernoulliの論文は、

Mensura sortis ad fortuitam successionem rerum naturaliter contingentium applicata (1770)

Continuatio argumenti de mensura sortis ad fortuitam successionem rerum naturaliter contingentium applicata (1771)

と2本あったようだが、電子化された原文を見つけられなかった。

D.Bernoulliによる正規近似の導出は、冒頭の論文に書かれてるけど、少し記号を変えて再構成すると以下のような議論。

正数 $N$ と自然数 $k \le N$ 、0〜1の実数 $p$ に対して、
$q(k) = \dfrac{N!}{k!(N-k)!} p^k(1-p)^{N-k}$
とすると、
$q( k+1 ) - q( k ) = \displaystyle \frac{N-k}{k+1} \displaystyle \frac{1-p}{p} q(k) - q(k)$
が成り立つ。 $M=pN$ 、 $k=M+x$ , $f(x)=q(k)=q(M+x)$ と記号を導入すると、
$f(x+1) - f(x) = -\dfrac{\frac{1}{1-p} x + 1}{M+x+1} f(x)$
と変形できる。 $f(x)$ を実数上定義された微分可能な関数であれば
$f(x+1)-f(x) = \displaystyle \int_{x}^{x+1} f'(t)dt$
なので、 $f'(t)$ が、 $(x,x+1)$ で、ほぼ一定であれば、微分方程式
$f'(x) = -\dfrac{\frac{1}{1-p} x + 1}{M+x+1} f(x)$
を得る。これは初期条件 $f(0)=Q$ の元で解くことができる。

関数方程式を微分方程式で近似するタイプの議論は、他にもいくつか存在する。例えば、
$F(x+1)-F(x) = x^k$
の解は、ベルヌーイ多項式だけど、
$F'(x) \sim x^k$
と置き換えると
$F(x) \sim \dfrac{x^{k+1}}{k+1} + C$
と解けて、ベルヌーイ多項式の主要項が得られる。

もうひとつの例として
$F(x+1)= x F(x)$
なら、対数を取って
$\log F(x+1) - \log F(x) = \log(x)$
なので、同様の近似によって
$\log F(x) \sim x \log(x) - x$
という形で、de Moivre-Stirlingの公式の主要項が得られる。

そうは言っても、この議論は、ちょっと荒っぽい。後者の議論に関連して、大学一年レベルの解析学の教科書には
$\displaystyle \int_{0}^{N} \log(x) dx \le \log(N!) \le \displaystyle \int_{1}^{N+1} dx \log(x)$
という不等式が書いてある。正規近似の導出でも、同じようなことができそうか考える。

関数方程式に立ち戻ると、解きたいのは
$f(x+1) - f(x) = -\dfrac{\frac{1}{1-p} x + 1}{M+x+1} f(x)$
だけど
$g(x) = 1 - \dfrac{\frac{1}{1-p} x + 1}{M+x+1} = \dfrac{M- \frac{p}{1-p}x}{M+x+1}$
と置くと
$f(x+1)=g(x) f(x)$
という形になる。

$x \lt \dfrac{1-p}{p} M = (1-p)N$ の時、 $g(x) \gt 0$ なので対数を取って、
$\log f(x+1) - \log f(x) = \log g(x)$
になる。 $x \gt -M-1=-(pN+1)$ で、 $g(x)$ は、単調減少関数なので
$\log g(k+1) \le \displaystyle \int_{k}^{k+1} \log g(x) dx \le \log g(k)$
が成り立つ。従って
$\displaystyle \sum_{k=1}^n \log g(k) \le \displaystyle \int_{0}^{n} \log g(x) dx \le \displaystyle \sum_{k=0}^{n-1} \log g(k)$
なので
$\log g(n) + \displaystyle \int_{0}^{n} \log g(x) dx \le \displaystyle \sum_{k=0}^{n} \log g(k) \le \log g(0) + \displaystyle \int_{0}^{n} \log g(x) dx$
となる。

更に
$\displaystyle \int_{n}^{n+1} \log g(x) dx \lt \log g(n)$
$\log g(0) \lt 0$
を使うと、
$\displaystyle \int_{0}^{n+1} \log g(x) dx \le \displaystyle \sum_{k=0}^{n} \log g(k) \le \displaystyle \int_{0}^{n} \log g(x) dx$
が示せる。

上側の評価は、似たような議論で、もう少し改善できて、
$\sum_{k=0}^{n} \log g(k) \le \displaystyle \int_{-1}^{n} \log g(x) dx \lt \displaystyle \int_{0}^{n} \log g(x) dx$
を示すことが出来る。

$x \lt 0$ の時も、同じような議論ができるが省略。以下、 $x \ge 0$ で考える。 $x$ は、実数値を取るとしてるが、例えば、 $f(x)$ が、 $x \gt 0$ で単調減少すると仮定すると

$\log f(x) - \log f(0) = \log f(x) - \log f(x-[x]) + \log f(x-[x]) - \log f(0)$

と分解( $[x]$ は、xを超えない最大の整数)して、 $|\log f(x-[x]) - \log f(0)| \lt | \log f(1) - \log f(0)|$ は、Nが大きくなると、0に近付く。ちゃんと計算すると、( $f(x)$ が、 $x \gt 0$ で単調減少するという仮定のもとで)xが自然数の時と同じように、xが任意の実数の時でも同様の不等式の成立を示すことができる。

数値例を一つ見ておく。ベータ関数 $B(z,w)$ を使うと、関数方程式の一つの厳密解を書くことができ、これに対して、
$E(x) = \log f(x) - \log f(0) = x \log(\dfrac{1-p}{p}) - \log B(M+x+1,\dfrac{1-p}{p}M - x + 1) + \log B(M+1,\dfrac{1-p}{p}M+1)$
という関数を定義する。ここで考えた近似解に対して同様の関数を計算すると、
$M(x) = \displaystyle \int_{0}^{x} \log g(x) dx = \dfrac{(w x -M) \log(\frac{M-wx}{M+x+1}) - (Mw+M+w)\log(\frac{M+x+1}{M+1})}{w} + \dfrac{M}{w} \log(\dfrac{M}{M+1})$
$w=\dfrac{p}{1-p}$
を得る。また、ベルヌーイによる微分方程式の解からは
$C(x) = -(1+w)x + (Mw+M+w)\log( \dfrac{M+x+1}{M+1} )$
が得られる。

これらの関数を $N=106,p=0.23$ の時、プロットすると、以下のようになる。

f:id:m-a-o:20211114225129p:plain

厳密解は、 $M(x)$ と $M(x-1)-M(-1)$ のほぼ中間を通り、 $\dfrac{M(x)+M(x-1)-M(-1)}{2}$ は、Nが比較的小さい時でも、厳密解とよく一致する。微分方程式解から得られる $C(x)$ は、xが0の近くにある時には、 $M(x)$ と近い値を取るが、xが大きくなると、乖離が目立つ。

正規近似を得るには、 $M(x)$ をTaylor展開して、二次の項まで取ればいい。大体の計算は、素直に展開してけばいい。
$x \log \left(\dfrac{M-wx}{M+x+1} \right) = x \log (\dfrac{M-wx}{M+1}) - x \log(\dfrac{M+x+1}{M+1})$
という分解を使うといい。このへんの議論は、冒頭の論説と本質的には変わらないので省略。

グラフ作成コード

import math
from scipy.special import betaln
import numpy as np
import matplotlib.pyplot as plt

N,p=106,0.23
M=p*N

def exact(x):
   return x * math.log(p/(1-p)) - betaln(M+x+1 , (1-p)*M/p - x + 1) + 0*(M*math.log(p) + ((1-p)*M/p) * math.log(1-p) - math.log(N+1)) + betaln(M+1 , (1-p)*M/p+1)


def approx(x):
    w = p/(1-p)
    C = (M*math.log(M/(M+1)) + (M*w+M+w)*math.log(M+1))/w
    return ((w*x-M)*math.log((M-w*x)/(M+x+1)) - (M*w+M+w)*math.log(M+x+1))/w + C


def C(x):
  w = p/(1-p)
  return -(1+w)*x + (M*w+M+w)*math.log((M+x+1)/(M+1))


if __name__=="__main__":
   x = (np.arange(4000)/100 - 20)
   y1 = np.vectorize(exact)(x)
   y2 = np.vectorize(approx)(x)
   y3 = np.vectorize(approx)(x-1)-approx(-1)
   y4 = np.vectorize(C)(x)
   y5 = (y2+y3)/2
   plt.plot(x,y1 ,label="E(x)" , linewidth=3 ,alpha=0.6)
   plt.plot(x,y2 , label="M(x)",linestyle="dashed")
   plt.plot(x,y3 , label="M(x-1)-M(-1)",linestyle="dashed")
   plt.plot(x,y4 , alpha=0.8 , label="C(x)" , linewidth=5 , linestyle="dotted")
   plt.plot(x,y5 , alpha=1 , label="(M(x)+M(x-1)-M(-1))/2" , linewidth=5 , linestyle="dotted")
   plt.legend()
   plt.show()