スターリングの近似

提供: testwiki
ナビゲーションに移動 検索に移動

テンプレート:No footnotes

log n! と n log nnn → ∞ のとき漸近する

スターリングの近似テンプレート:Lang-en-short)またはスターリングの公式テンプレート:Lang-en-short)は、階乗、あるいはその拡張の一つであるガンマ関数漸近近似である。名称は数学者テンプレート:仮リンクにちなむ。

概要

スターリングの近似は精度に応じていくつかの形がある。応用上よく使われる形の公式は、ランダウの記号を用いて、

logn!=nlognn+O(logn)

である。テンプレート:Math における次の項は テンプレート:Math である。故に、次によい近似のテンプレート:仮リンク

n!2πn(ne)n

である[1]。(ここで記号 は両辺の比が(n → ∞ のとき) 1 に収束することを意味する。)

テンプレート:Math の漸近近似よりもむしろ上下からの評価が必要なことがある。そのような評価として、任意の正の整数 テンプレート:Math に対して

2πnn+1/2enn!enn+1/2en

が成り立ち[2]、従って任意の テンプレート:Math に対して比 テンプレート:Mathテンプレート:Mathテンプレート:Math の間にある。

スターリングの近似は階乗の複素引数への拡張の一つであるガンマ関数 テンプレート:Math(正の整数 テンプレート:Mvar に対し テンプレート:Math が成り立つ;ボーア・モレルップの定理も参照)に拡張することができ、

Γ(z)2πz(ze)z(|argz|πε,|z|)

が成り立つ(ただし テンプレート:Mathテンプレート:Sfnテンプレート:Math のときは収束が遅くなるため、応用上は相補公式などを用いて テンプレート:Math 程度に制限することが多い。

導出

初等的な導出

スターリングの公式の厳密な証明にはオイラーの和公式、あるいは鞍点法といった複素解析の技法などを用いられることが多いが、初等的に導くことも可能である。まず階乗の対数を積分で近似する。logが凹関数であることから k-1<x<k (k=2,3,...) に対して

logk(kx){logklog(k1)}<logx<logk1k(kx)

これを k-1 から k まで積分して

logk12{logklog(k1)}<k1klogxdx<logk12kk1klogxdx+12k<logk<k1klogxdx+12{logklog(k1)}

k=m+1,m+2,...,n に対して足し合わせると

logn!m!=k=m+1nlogk>mnlogxdx+k=m+1n12k>mnlogxdx+12n12m+mn12xdx=(n+1/2)lognn+1/(2n)(m+1/2)logm+m1/(2m)logn!m!=k=m+1nlogk<mnlogxdx+12{lognlogm}=(n+1/2)lognn(m+1/2)logm+m
nn+1/2en+1/(2n)mm+1/2em+1/(2m)<n!m!<nn+1/2enmm+1/2em

ここで an=n!nn+1/2en と定めると

e1/(2n)e1/(2m)<anam<1

m,n→∞のとき最左辺は1に収束するから、特に n=2m のとき

limma2mam=1

これとウォリスの公式の系 limn4n(n!)2n(2n)!=π と比較すると、

4n(n!)2n(2n)!=4n(annn+1/2en)2na2n(2n)2n+1/2e2n=an22a2n
limnan=limn4n(n!)2n(2n)!2a2nan=2π

を得る。

精度の改善

精度を改善するために an を評価する。

logamam+1=logm!(m+12)logm+mlog(m+1)!+(m+32)log(m+1)m1=(m+12)log(1+1m)1=(m+12)k=1(1)k1k1mk1=k=0(1)kk+11mk+k=1(1)k12k1mk1=k=2(k1)(1)k2k(k+1)1mk
logan2π=m=nlogamam+1=k=2(k1)(1)k2k(k+1)m=n1mk=k=2(1)k2k(k+1){1nk1+O(1nk)}=112n+O(1n2)

従って

n!2πn(ne)nexp(112n)

オイラーの和公式による導出

オイラーの乗積表示によるガンマ関数の定義の対数をとり

logΓ(z)=log{(z1)Γ(z1)}=limN(z1)logN+n=1N{lognlog(n+z1)}

f(n)=lognlog(n+z1)オイラーの和公式を適用すれば

logΓ(z)=limN(z1)logN+n=1Nf(n)dn+12(f(N)+f(1))+k=1mB2k(2k)!(f(2k1)(N)f(2k1)(1))+n=1NB2m+1(nn)(2m+1)!f(2m+1)(n)dn=limN(z1)logN+[nlognn(n+z1)log(n+z1)+(n+z1)]n=1N+12{logNlog(N+z1)logz}+k=1mB2k(2k)(2k1){1N2k11(N+z1)2k11+1z2k1}+n=1NB2m+1(nn)2m+1{1n2m+11(n+z1)2m+1}dn=limN(N+z12){logNlog(N+z1)}+(z12)logz+k=1mB2k(2k)(2k1){1N2k11(N+z1)2k11+1z2k1}+n=1NB2m+1(nn)2m+1{1n2m+11(n+z1)2m+1}dn=z+1+(z12)logzk=1mB2k(2k)(2k1)(11z2k1)+n=1NB2m+1(nn)2m+1{1n2m+11(n+z1)2m+1}dn

となる。右辺の定数を集めて

C=1k=1mB2k(2k)(2k1)+n=1NB2m+1(nn)dn(2m+1)n2m+1

とすれば

logΓ(z)=Cz+(z12)logz+k=1mB2k(2k)(2k1)z2k1n=1NB2m+1(nn)dn(2m+1)(n+z1)2m+1

となり、この主要部をガンマ関数の相補公式に代入して zi とすれば

Γ(z)Γ(1z)=zΓ(z)Γ(z)=πsinπz
πi+logz+logΓ(z)+logΓ(z)logsinπzlogπ=0
πi+logz+Cz+(z12)logz+C+z+(z12)(πi+logz)logsinπzlogπ=0
2C+πi2πizlogsinπzlogπ=0

となるが

logsinπz=log(eπizeπiz)πi2log2πizπi2log2

であるから

C=log2π2

を得る。剰余項については

α=21+cosargz

として

|n=1NB2m+1(nn)dn(2m+1)(n+z1)2m+1||B2m|2mn=1Ndn|n+z1|2m+1|B2m|2mα2m+1n=1Ndn(n+|z|1)2m+1=|B2m|α2m+1|z|2m=O(z2m)

である。故に テンプレート:Indent を得る。最初の数項を書き下せば

logΓ(z)log2πz+(z12)logz+112z1360z3+11260z511680z7+11188z9
Γ(z)2πz(ze)zexp(112z1360z3+11260z511680z7+11188z9)

とやり、指数関数のテイラー展開により

Γ(z)2πz(ze)z(1+112z+1288z213951840z3)

となる。

鞍点法による導出

スターリングの公式は鞍点法の好適例とされることが多いが、実際に複素平面全体(負の実数を除く)で漸近近似が成立することを鞍点法によって示すのは困難であるから、ここでは z を正の実数に限定する。ガンマ関数t=z(1+u) の置換により

Γ(z+1)=0tzetdt=1zz(1+u)zezzuzdu=zz+1ez1ez{ulog(1+u)}du=zz+1ez(1εez{ulog(1+u)}du+εεez{ulog(1+u)}du+εez{ulog(1+u)}du)(ε1)

となるが、z が十分に大きければ u=0 の付近が支配的であるから

Γ(z+1)zz+1ezεεez{ulog(1+u)}duzz+1ezεεezu22duzz+1ezezu22du

という近似が許され、ガウス積分により

Γ(z+1)zz+1ez2πz=2πz(ze)z

を得る。ε=z13 として、近似の誤差は

|1εez{ulog(1+u)}du|1ε1ε|u(1+u)ez{ulog(1+u)}|du=1εz[ez{ulog(1+u)}]1ε=1εz{ez(ε22+O(ε3))}0z23ez1/32
|εez{ulog(1+u)}du|2εε|2u(1+u)εez{ulog(1+u)}|du=2εz[ez{ulog(1+u)}]ε=02εz{ez(ε22+O(ε3))}2z23ez1/32
εεez{ulog(1+u)}du=εεez(u22u33+O(u4))du=εεezu22ez(u33+O(u4))du=εεezu22(1+zu33+zO(u4))du=εεezu22(1+zO(u4))du
|εεez{ulog(1+u)}duεεezu22du||zO(ε4)εεezu22du|=|O(z13)εεezu22du|

であり

z23ez1/32z12(z)

であるから

Γ(z+1)=2πz(ze)z(1+O(z13))

を得る。これは

limzΓ(z+1)2πz(ze)z=1(|argz|<π)

を示すに十分である。ただし、実際の誤差は O(z1) であるが、それを鞍点法で示すのは困難である。

収束の速度と誤差見積もり

より正確に記すと、次のようになる[3]

n!=2πn(ne)neλn

ここで

112n+1<λn<112n.

スターリングの公式は以下の級数(スターリング級数)の近似(初項で打ち切ったもの)である。

n!2πn(ne)n(1+112n+1288n213951840n35712488320n4+)

n としたとき、省かれた級数はその最初の項とそれ以降が相殺するように漸近していく。これは漸近展開の一例である。

以下のような階乗の対数の漸近展開も「スターリング級数」と呼ぶ。

lnn!nlnnn+12ln2πn+112n1360n3+11260n511680n7+

この場合、誤差は打ち切った級数の初項と同じ符号で同程度の大きさであることが知られている。

ガンマ関数に対するスターリングの公式

すべての正の整数に対して、

n!=Π(n)=Γ(n+1)

が成り立つ。ここで Γ はガンマ関数を表す。

しかしながら、テンプレート:仮リンク Π(z)=Γ(z+1) は、階乗とは異なり、より広く、正でない整数を除いてすべての複素数に対して定義される。それにもかかわらず、スターリングの公式をなお適用することができる。Re z > 0 であれば

logΓ(z)=(z12)logzz+12log2π+20arctan(t/z)e2πt1dt

が成り立つテンプレート:Sfn。 部分積分を繰り返すことで次が得られる

logΓ(z)(z12)logzz+12log2π+n=1B2n2n(2n1)z2n1.

ここで Bテンプレート:Subn 番目のベルヌーイ数である。(無限和は収束しないので、この公式は漸近展開にすぎないことに注意する。)公式はεを正数として |arg z| < テンプレート:Π − ε であるときに絶対値の十分大きい z に対して成り立つ。公式の右辺に現れる級数はスターリング級数と呼ばれるテンプレート:Sfn。最初の m 項が使われるとき誤差項は O(z2m1) である。対応する近似は

Γ(z)=2πz(ze)z(1+O(1z))

のように書ける。この漸近展開のより進んだ応用は Re z が定数の複素変数 z に対してである。例えば直線 テンプレート:Sfrac + it 上でテンプレート:仮リンクの Im z において適用されたスターリングの公式を見よ。

ビネーの公式

スターリングの公式は収束しない級数を伴うので解析的に扱いづらいが、収束しない級数を収束する積分に換えたものとしてビネーの(第二)公式があるテンプレート:Sfn

Γ(z+1)=2πz(ze)zeμ(z),μ(z)=20arctan(t/z)e2πt1dt,(z>0)

ビネーの公式は、スターリングの級数を形式的に(収束条件を無視して)操作することによっても導かれるが、厳密には対数ガンマ関数の導関数にアーベル・プラナの和公式を適用して得られる。 テンプレート:Indent

ddzlogΓ(z)=limNlogNn=0N1n+z
d2dz2logΓ(z)=n=01(n+z)2

z>0 なら f=(n+z)2 は右半平面において正則であるからプラナの和公式により

d2dz2logΓ(z)=01(n+z)2dt+12z2+i0(z+it)2(zit)2e2πt1dt=1z+12z2+i0(z+it)2(zit)2e2πt1dt

積分して

ddzlogΓ(z)=C1+logz12z+i0(z+it)1+(zit)1e2πt1dt
logΓ(z)=C2+C1z+zlogzlogz2+i0log(z+it)+log(zit)e2πt1dt=C2+C1z+(z12)logz+20arctan(t/z)e2πt1dt

z>0 なら|tan(t/z)|<M は有界であるから

|0arctan(t/z)e2πt1dt||0|z|1/2arctan(t/z)e2πt1dt|+||z|1/2arctan(t/z)e2πt1dt|0|z|1/2k=1|t/z|kdt2πt+M|z|1/2e2πdt(e2π1)ete2π(|z|>1)0|z|1/2dt2π(|z|t)+Me2πe2π1t=|z|1/2e2πtdt
limz|0arctan(t/z)e2πt1dt|limz12π[log(|z|t)]0|z|1/2+limzMe2π2π(e2π1)[e2πt]|z|1/2=0

である。 スターリングの公式と比較して積分定数を求め

logΓ(z)=12log2πz+(z12)logz+20arctan(t/z)e2πt1dt

真数に直して

Γ(z+1)=zΓ(z)=2πz(ze)zeμ(z),μ(z)=20arctan(t/z)e2πt1dt,(z>0)

を得る。なお、ビネーの公式を元にして部分積分を繰り返すとスターリングの級数が得られる。

収束級数形式のスターリングの公式

トーマス・ベイズの John Canton への書簡が1763年王立協会により公表されている。それによると、スターリングの公式は収束級数ではないとされていた[4]

スターリングの公式の収束級数形式を得るには以下を評価する。

02arctan(t/z)exp2πt1dt=lnΓ(z)(z12)lnz+z12ln2π

一つの方法として、階乗冪の逆数の収束級数を使う方法がある。zn=z(z+1)(z+n1) としたとき、次のようになる。

02arctan(t/z)exp2πt1dt=n=1cn(z+1)n

ここで

cn=1n01xn(x12)dx

である。以上から次のようなスターリング級数が得られる。

lnΓ(z)=(z12)lnzz+12ln2π+112(z+1)+112(z+1)(z+2)+59360(z+1)(z+2)(z+3)+2960(z+1)(z+2)(z+3)(z+4)+

これは、z>0 のとき収束する。

計算機向けの変形

ガンマ関数の(関数電卓などの)計算機向けの近似として次の式がある。

Γ(z)2πz(zezsinh1z+1810z6)z

これは、次と同等である。

2lnΓ(z)ln2πlnz+z{2lnz+ln(zsinh1z+1810z6)2}

これらはスターリングの公式を組み替えて、その結果生じる冪級数と双曲線正弦関数のテイラー展開の間の合致を観察することで得られる。この近似は z の実数部が 8 以上のとき、小数点以下 8 桁を超える精度を持つ。2002年、Robert H. Windschitl がリソースの制限された計算機(電卓など)でのそれなりの正確性を持った近似としてこれを示した[5]

Gergő Nemes は 2007年にほぼ同程度の結果を与える近似式を提案した。こちらはより単純である。

Γ(z)2πz{1e(z+112z110z)}z

これは、次と同等である。

lnΓ(z)12(ln2πlnz)+z{ln(z+112z110z)1}

歴史

この公式は最初に次の形でアブラーム・ド・モアブルにより1730年に発見された[6]

n![constant]nn+12en

スターリングの貢献は定数が 2π であることを示したことである[6]。より正確な形式はジャック・ビネが見出した。

スターリングの近似の「一次」バージョン n!=nn は、マックス・プランクが1901年の黒体放射の論文で使用した。これは多量の光子や振動子についての黒体放射エネルギーの方程式にリンクしている。この近似は量子論でよく使われ、例えばピーター・デバイルイ・ド・ブロイも使っている。アルベルト・アインシュタインサティエンドラ・ボースは違う方式を採用した。非常に大きな n について確率分布をグラフに描画してみると、両者はほぼ平行になる。

脚注

テンプレート:Reflist

参考文献

関連項目

外部リンク