« 2017年4月 | トップページ | 2017年6月 »

2017年5月

2017年5月 3日 (水)

Stirling の公式(階乗の近似公式)の計算(PARI/GPにより)

■Stirling の公式(階乗の近似公式)の計算
計算を表示させてみましたが,わからないことが多く,
数値の精度は確認できていません.

●PARI/GP home
「Pari32-2-9-2.exe (83472 KB), Mar 22 2017 」をDL
http://pari.math.u-bordeaux.fr/

・近似式1
n!=n^n*e^(-n)
log(n!)=n*logn-n

・近似式2
n!=n^(n)*e^(-n)*sqrt(2*pi*n)
log(n!)=n*log(n)-n+(1/2)*log(2*Pi*n)

・近似式3
n!=n^(n)*e^(-n)*sqrt(2*pi*n)*(1+1/(12*n))
log(n!)=n*log(n)-n+(1/2)*log(2*Pi*n)+log(1+1/(12*n))

※n!=gamma(n+1)も表示

●gpのプロンプトで対話的に実行する場合
ファイルを読み込んで実行(画面出力)
ファイルをフルパスで指定

(14:57) gp > \r C:\GP\my_stirling.txt⏎

Stirling の公式(階乗の近似公式)の計算

n         log(n!)                                              lngamma(n+1)
1000   5912.128178488163348878130887   5912.128178488163348878130887
5000   37591.14350887676656917322088   37591.14350887676656917322088
10000   82108.92783681435345538503006   82108.92783681435345538503006
15000   129242.8090480940091717829996   129242.8090480940091717829996
20000   178075.6217371987003128679282   178075.6217371987003128679282
25000   228171.7598536769082061852894   228171.7598536769082061852894
30000   279274.6532369700774206467427   279274.6532369700774206467427
40000   383871.6065818260202389663666   383871.6065818260202389663666
50000   490995.2430498562320145778187   490995.2430498562320145778187
60000   600132.4104620969598935528420   600132.4104620969598935528420
70000   710944.0335371889194936427603   710944.0335371889194936427603
80000   823189.1169230131924946281581   823189.1169230131924946281581
90000   936687.4681600499774973443934   936687.4681600499774973443934
100000   1051299.221899121865129278108   1051299.221899121865129278108

n         log(n!)                                              n*log(n)-n
1000   5912.128178488163348878130887   5907.755278982137052053974364
5000   37591.14350887676656917322088   37585.96595708118713327366849
10000   82108.92783681435345538503006   82103.40371976182736071965819
15000   129242.8090480940091717829996   129237.0822012652067707496840
20000   178075.6217371987003128679282   178069.7510507225609097839588
25000   228171.7598536769082061852894   228165.7775962584450313873258
30000   279274.6532369700774206467427   279268.5798193287728240163317
40000   383871.6065818260202389663666   383865.3893238429341962572025
50000   490995.2430498562320145778187   490988.9142205141555336362576
60000   600132.4104620969598935528420   600125.9904722542642130665906
70000   710944.0335371889194936427603   710937.5364722047228824122994
80000   823189.1169230131924946281581   823182.5530924814931458929747
90000   936687.4681600499774973443934   936680.8454381161906976210663
100000   1051299.221899121865129278108   1051292.546497022842008995727

n         log(n!)                                              n*log(n)-n+(1/2)*log(2*Pi*n)
1000   5912.128178488163348878130887   5912.128095154832793321781681
5000   37591.14350887676656917322088   37591.14349221009992472877618
10000   82108.92783681435345538503006   82108.92782848102012482947450
15000   129242.8090480940091717829996   129242.8090425384536170504893
20000   178075.6217371987003128679282   178075.6217330320336465484837
25000   228171.7598536769082061852894   228171.7598503435748730297338
30000   279274.6532369700774206467427   279274.6532341922996429718456
40000   383871.6065818260202389663666   383871.6065797426869056764360
50000   490995.2430498562320145778187   490995.2430481895653479333743
60000   600132.4104620969598935528420   600132.4104607080710046768132
70000   710944.0335371889194936427603   710944.0335359984433031746683
80000   823189.1169230131924946281581   823189.1169219715258279669168
90000   936687.4681600499774973443934   936687.4681591240515714222779
100000   1051299.221899121865129278108   1051299.221898288531795947553

n         log(n!)                                    n*log(n)-n+(1/2)*log(2*Pi*n)+log(1+1/(12*n))
1000   5912.128178488163348878130887   5912.128178484694097322071836
5000   37591.14350887676656917322088   37591.14350887662770404974455
10000   82108.92783681435345538503006   82108.92783681431873613348564
15000   129242.8090480940091717829996   129242.8090480939937405644351
20000   178075.6217371987003128679282   178075.6217371986916326837074
25000   228171.7598536769082061852894   228171.7598536769026508198573
30000   279274.6532369700774206467427   279274.6532369700735627320765
40000   383871.6065818260202389663666   383871.6065818260180688738945
50000   490995.2430498562320145778187   490995.2430498562306257126953
60000   600132.4104620969598935528420   600132.4104620969589290604223
70000   710944.0335371889194936427603   710944.0335371889187850346411
80000   823189.1169230131924946281581   823189.1169230131919520992380
90000   936687.4681600499774973443934   936687.4681600499770686790583
100000   1051299.221899121865129278108   1051299.221899121864782058857
(14:58) gp >

● windowsのプロンプトを起動して,リダイレクションを使って実行する場合
gp.exeへパス(path)を設定しておく.

a) my_stirling.txtはC:\GPにある.画面出力.
C:\GP>gp -q < my_stirling.txt⏎

b) my_stirling.txtはC:\GPにある.出力ファイル(_res_my_stirling.dat)はC:\GPに保存される
C:\GP>gp -q < my_stirling.txt > res_my_stirling.dat⏎

C:\GP>type res_my_stirling.dat⏎     【res_my_stirling.datの中身を出力】

●使用したmy_stirling.txt
わからないことが多く,メモ用にかいたものです.
確認して,適宜修正してください.
/*      */のコメントに簡単な説明を記述しました.
「my_stirling.txt」 (SHIFT-JIS)

●参考資料
・Stirlingの公式
https://oku.edu.mie-u.ac.jp/~okumura/stat/stirling.html
・ガンマ分布の中心極限定理とStirlingの公式
http://www.math.tohoku.ac.jp/~kuroki/LaTeX/20160501StirlingFormula.pdf

« 2017年4月 | トップページ | 2017年6月 »

無料ブログはココログ
2017年7月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31