ベンフォードの法則

世の中にあふれるさまざまな数値の先頭一桁（0は除く）の分布を調べると，
1～9まで，等確率ではありません．1が最も多く，数字が大きくなるほど出現しにくくなります．
その確率は，先頭一桁の数字をdとすると，

P(d)=log10(1+1/d)

となります．
下記サンプルプログラムは，適当な数値に適当な数値を掛け続けることで得られる数値の
先頭一桁に，１が一番出やすいことを数値実験で示しています．

**************************************************************
benford.sas
090815 翔
ベンフォードの法則
**************************************************************;

data benford;
  num=123.45;
  do i=1 to 10000;
    r=(int(99*ranuni(0))+1); /*1～99までの一様整数乱数を生成*/
    num=num*r;                      /*元の数字に乱数を掛けて*/
    num=num/10**int(log10(num));    /*桁を落として10未満にする*/
    d=int(num);                     /*最初の一桁の数字*/
    output;
  end;
run;

proc freq;tables d;run;

/*
                                 累積         累積
d      度数      パーセント      度数      パーセント
-----------------------------------------------------
1        3055       30.55          3055       30.55
2        1801       18.01          4856       48.56
3        1272       12.72          6128       61.28
4         951        9.51          7079       70.79
5         754        7.54          7833       78.33
6         662        6.62          8495       84.95
7         538        5.38          9033       90.33
8         517        5.17          9550       95.50
9         450        4.50         10000      100.00
 */

data p;
  do d=1 to 9;
    p=log10(1+1/d);
    output;
  end;
run;

proc print;run;

/*
理論値

OBS    d       p

 1     1    0.30103
 2     2    0.17609
 3     3    0.12494
 4     4    0.09691
 5     5    0.07918
 6     6    0.06695
 7     7    0.05799
 8     8    0.05115
 9     9    0.04576
 */

直感的に，「等確率でないこと」が不自然に感じるかもしれません．
でも，1が出やすい傾向は，九九表にもみることができます．

九九表
1	2	3	4	5	6	7	8	9
2	4	6	8	10	12	14	16	18
3	6	9	12	15	18	21	24	27
4	8	12	16	20	24	28	32	36
5	10	15	20	25	30	35	40	45
6	12	18	24	30	36	42	48	54
7	14	21	28	35	42	49	56	63
8	16	24	32	40	48	56	64	72
9	18	27	36	45	54	63	72	81
九九の先頭一桁表
1	2	3	4	5	6	7	8	9
2	4	6	8	1	1	1	1	1
3	6	9	1	1	1	2	2	2
4	8	1	1	2	2	2	3	3
5	1	1	2	2	3	3	4	4
6	1	1	2	3	3	4	4	5
7	1	2	2	3	4	4	5	6
8	1	2	3	4	4	5	6	7
9	1	2	3	4	5	6	7	8

data kuku;
  do i=1 to 9;do j=1 to 9;
    d=i*j;
    if d>9 then d=int(d/10);
    output;
  end;end;
run;

proc freq;tables d;run;

/*
                                 累積         累積
d      度数      パーセント      度数      パーセント
-----------------------------------------------------
1          18       22.22            18       22.22
2          15       18.52            33       40.74
3          11       13.58            44       54.32
4          12       14.81            56       69.14
5           6        7.41            62       76.54
6           7        8.64            69       85.19
7           4        4.94            73       90.12
8           5        6.17            78       96.30
9           3        3.70            81      100.00
 */