フィボナッチ数列の各項の先頭の数字の度数分布はどのようになると思うか?
そんなの簡単。どの数字も同じ頻度。
ブッブー。違います。1が一番多く,順に少なくなり,9が一番少ない。
ベンフォードの法則という。
Python だと,int 変数は自動的に精度が確保されるので,特別な仕掛なしでフィボナッチ数列の正確な値が簡単に計算できる。
import scipy as sp
n = 100000
tbl = sp.zeros(10, dtype=int)
a = b = 1
tbl[1] = 2
for i in range(2, n):
a, b = b, a+b
tbl[int(str(b)[0])] += 1
print(tbl[1:])
これによって,先頭桁が 1 ~ 9 である項数が以下のようであることが分かる。
[30103 17610 12494 9690 7918 6695 5798 5117 4575]
理論確率 p は log10(i+1) - log10(i), i = 1, ..., 9
期待値は n * p
Observed log10 p Expected
1 30103 0.301030 0.301030 30102.999566
2 17610 0.477121 0.176091 17609.125906
3 12494 0.602060 0.124939 12493.873661
4 9690 0.698970 0.096910 9691.001301
5 7918 0.778151 0.079181 7918.124605
6 6695 0.845098 0.066947 6694.678963
7 5798 0.903090 0.057992 5799.194698
8 5117 0.954243 0.051153 5115.252245
9 4575 1.000000 0.045757 4575.749056
実に良く一致している。
ちなみに,第 100000 項は 上のプログラムを実行した後 print(b) とすれば表示されるが,259740 から始まる 20899 桁の数である。