階乗の先頭桁の数字の分布

では次に，n の階乗 n!, n=1, 2, ... 100000 の先頭桁の数字の分布はどのようになるか？

ああ，簡単ですね。前のプログラムをまねて，っと...

import scipy as sp

n = 100000
tbl = sp.zeros(10, dtype=int)
f = 1
for i in range(1, n+1):
  f *= i
  tbl[int(str(f)[0])] += 1
print(tbl[1:])

..........

5時間回してもまだ答えが出ません（泣く）

だろうね。ずいぶん余分な処理をやっているからねえ。
それにしても，いきなり n = 100000 でやったんじゃないだろうね。

もちろんです。 n=1000, はあっという間に終わったので，n = 10000 で 50 秒くらいだったので，大丈夫と思ったんですけど。

概算だけど，n が 2 倍になると，計算時間は 10 倍くらいになりそうだから，100000 だと 20 時間以上かかるんじゃないか？

ひえー。何とかなりませんか。

なるね。まず，先頭桁の数字だけが必要なんだから，いくら Python なら正確な整数が扱えるといっても，無駄なことをしないこと。だって，100000 の階乗は，28242 から始まる 456574 桁の数になる。これを全桁バラバラの数字にしてその最初の要素（数字）だけを使う str(f)[0] のは，いかにもにも無駄。

ある大きな数 x の先頭桁はどうやって求める？

print(x) とすれば，わかります。

いやいや，そんな（笑）
log10(x) の「整数部+1」が桁数，「10^小数部」の整数部が先頭桁の数字。
x = 3141592653589793
とすると，
sp.log10(x) は 15.497149872694134 だから，桁数は 15+1 つまり 16 桁
10**0.497149872694134 = 3.1415926535897944 だから先頭の数字は 3 だとわかる。

ああ，高校の数学で対数の応用で習いましたね。思い出しました。

もう一つ。n の階乗は Γ(n＋1) に等しい。

Γ って何ですか？

ガンマ関数。
import math
math.factorial(4)
math.gamma(5)
math.factorial(100)
math.gamma(101)
をやってみればわかるね。
100 の階乗も Python だと math.factorial(100) で正確に全桁分かるけど，何度もいうように全桁は必要ないので，math.gamma(101) で十分。

なるほど....
でも，math.gamma(100001) だとエラーが出ました。OverflowError ですって。

とても大きい数になるので，計算できませんということだね。そんなこともあるので，結果の対数を返してくれる math.lgamma というのがある。

math.lgamma(100001) をやってみましたが，1051299.221899122 です。さっき，「100000 の階乗は，28242 から始まる 456574 桁の数になる」といわれたのですが，1051299+1 つまり1051300 桁で，先頭は 10**0.221899122 = 1.6668599890439895 だから，先頭数字は 1 じゃないですか？

気が早いね。「結果の対数を返してくれる」とはいったが，「結果の常用対数」を返してくれるとは言わなかっただろう？返してくれるのは，「結果の自然対数」だ。理数系では対数といえば自然対数を意味するのでちょっとまぎらわしかったね。

「結果の常用対数」を返してくれる関数はないのですか？

ないんだね。でも，「対数の底の変換公式」を覚えていないかい？

ああ，そうそう。log2(x) = log10(x)/log10(2) みたいなやつですね。
常用対数は log10(x) のように書けるとして，自然対数はどう書けばよいですか。そもそも，自然対数の底は幾つです？

ln(x) とかいたり，log(x) と書いたりするよ。自然対数の底は，2.718281828459.... だけど，e で表す。コンピュータ上では Python なら math.e として使えるよ。

ずいぶん半端な数ですが，それは置いておくことにします。
じゃあ，log(x) = log10(x)/log10(e) だから，結果の常用対数は log10(x) = log10(e) * log(x) となるので，math.log10(math.e)* math.lgamma(100001) が欲しい結果ですね。
456573.45089997095 となりましたから，桁数は 456573+1 つまり 456574 です（合ってます）
10**0.45089997095 = 2.8242294082311297 ですから，先頭数字は 2 で，そのあと 8, 2, 4, 2, 2, 9 と続くのですね（合ってます）

では，この考え方でプログラムするとどうなるかな？まずは n = 1000 でやってみよう。

はい。

import scipy as sp
import math

const = math.log10(math.e)
tbl = sp.zeros(10, dtype=int)
n = 1000
for i in range(1, n+1):
  a = const*math.lgamma(i+1)
  a = a-math.floor(a)
  b = int(str(math.floor(10**a))[0])
  tbl[b] = tbl[b]+1
print(tbl[1:])

結果は，
[294 175 124 102  69  87  51  51  47]
となりました。

前の，とてつもなく遅いプログラムでは n = 1000 のときの結果はどうだった？

はい，もう一度やってみます。
[293 176 124 102  69  87  51  51  47]
あれ，'1' と '2' の頻度が一つずつ違います。
おかしいなぁ。 n = 2000 でも，n = 4000 でも同じように '1' と '2' の頻度が一つずつ違います。

そうだね。for ループの中に i が 10 以下のときに 10**a がどのような値になるかを書き出してご覧。

for i in range(1, n+1):
  a = const*math.lgamma(i+1)
  a = a-math.floor(a)
  if i <= 10:
      print(i, 10**a)

のようにですね。
実行してみると...
1 1.0
2 1.9999999999999993
3 6.000000000000002
4 2.399999999999998
5 1.2000000000000008
6 7.200000000000008
7 5.040000000000002
8 4.03199999999999
9 3.628799999999987
10 3.6287999999999943
あれ，n = 2 のとき，本当は 2 になるところが 1.9999999999999993 となってしまって，先頭桁が '2' ではなくて '1' だとされています。

そうだね。ここがコンピュータで数値計算をするときに気をつけないといけないところだ。コンピュータは 2 の羃乗数の和で表せる数 0.625 = 0.5 + 0.125 = 2^-1 + 2^-3, 13 = 1 + 4 + 8 = 2**0 + 2**2 + 2**3 は正確に保存できるが，そうでない一般の数はあくまでも近似値が保存されているに過ぎない。
対処法は場合場合によるが，今回のように先頭桁の数字が必要ならば 10**a に極々小さい数を加えてやれば良い。たとえば 1e-14 を加えると，
1 1.00000000000001
2 2.0000000000000093
3 6.0000000000000115
4 2.4000000000000083
5 1.2000000000000108
6 7.200000000000018
7 5.040000000000012
8 4.032
9 3.6287999999999974
10 3.6288000000000045
のようになり，n = 2 の場合にも，先頭桁は正しく '2' になる。他の場合には 1e-14 を加えて先頭桁の数字が変わるということはない。

では，最終的なプログラムは

import scipy as sp
import math

const = math.log10(math.e)
tbl = sp.zeros(10, dtype=int)
n = 100000
for i in range(1, n+1):
  a = const*math.lgamma(i+1)
  a = a-math.floor(a)
  b = int(str(math.floor(10**a+1e-14))[0])
  tbl[b] = tbl[b]+1
print(tbl[1:])

で，結果は，[29925 17691 12655  9760  7888  6683  5781  5083  4534]
となりました。
あっという間に計算できましたね。
あれれ，ベンフォードの法則とは微妙にずれてますね。
フィボナッチ数列の場合は，ほぼ理論通り
[30103 17610 12494  9690  7918  6695  5798  5117  4575]
でしたから。

本当だね〜。なんでかな？

 

なお，R の場合だと，極々小さい数を足したりしなくてもよい。

> const = log10(exp(1))
> tbl = numeric(9)
> n = 100000
> for (i in 1:n) {
+   a = const*lfactorial(i)
+   a = a-floor(a)
+   b = as.integer(unlist(strsplit(as.character(floor(10^a)), "")))
+   tbl[b] = tbl[b]+1
+ }
> print(tbl)
[1] 29925 17691 12655  9760  7888  6683  5781  5083  4534

また，そもそも for ループなどは必要なく，ベクトル演算を用いて，より短く書くことができる。

> n = 100000
> a = const*lfactorial(1:n)
> a = a-floor(a)
> table(as.integer(unlist(strsplit(as.character(floor(10^a)), ""))))
b
    1     2     3     4     5     6     7     8     9
29925 17691 12655  9760  7888  6683  5781  5083  4534

フィボナッチ数列

フィボナッチ数列の各項の先頭の数字の度数分布はどのようになると思うか？

そんなの簡単。どの数字も同じ頻度。

ブッブー。違います。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 桁の数である。

フィボナッチ数列の一般項（その2）

R で，実数の内部表現をプリントさせてみると，4項目から既に，末尾にゴミが入る。

これは，多倍長演算をしても，例えビット数を無限にしても最下位にゴミが入る（だって，2の羃乗数でない限り2進数では正確に表せないのだから，あたりまエ）

第 4 項
> a
[1] 3.0000000000000004441
> print(sprintf("%a", a))
[1] "0x1.8000000000001p+1"
> a==3
[1] FALSE

第 5 項
> a
[1] 5.0000000000000008882
> print(sprintf("%a", a))
[1] "0x1.4000000000001p+2"
> a==5
[1] FALSE
：
第 10 項
> a
[1] 55.000000000000014211
> print(sprintf("%a", a))
[1] "0x1.b800000000002p+5"
> a==55
[1] FALSE

フィボナッチ数列の一般項

「フィボナッチ数列の一般項の実数演算結果がExcelとLibreOffice Calcで違う」ということだ。

まあ，そうなんだろうけど

議論の展開とは別だが，実際問題として，結果は整数になって欲しいんだから，演算誤差があろうと，整数に丸めてやれば何の問題もない。というか，整数に丸めない方が異常だ。

ちゃんとやっていれば，実数演算結果が違うってことに気づかなかったはずだワン。

n = 70
A = 1:n
B = round((((1+sqrt(5))/2)^A-(((1-sqrt(5))/2)^A))/sqrt(5))
C = numeric(n)
C[1] = C[2] = 1
for (i in 3:n) {
    C[i] = C[i-2] + C[i-1]
}
B==C
B[n]

n = 70 まで正しい。

防衛省の調査ミス（その3）

件の本山だけど，標高も間違えていた

敵性の Google なんか使うからだ。わが国には国土地理院があるぞ。本山の標高は 715 メートルだ。

防衛省の調査ミス（その2）

> 遮蔽（しゃへい）となる尾根などがあった場合はそこの標高、ない場合は山頂の標高、さらに国有地までの水平距離を紙の上で定規を使って計測。高さと水平距離の縮尺の違いに気付かないまま三角関数を用いて計算し、誤った仰角を算出した

ああ，一応は三角関数を使ったんだ。

でも，高さと距離は定規で測った。

ちなみに，定規とは直線を引くためのもの，長さを測るのは物差し。

まあ，測った対象が間違っていたということ。でもね，測った山の高さと水平距離の数値を見れば，？？？となるのが普通でしょ。

そもそも，標高も水平距離も数値で出ているのだから，それをわざわざ定規（物差し）で測ったということが納得いかねぇ。

防衛省の調査ミス

調査ミスではなく，初歩的な計算ミス
指摘する新聞なども，基本的なことが分かってるのか？という懸念がある。

標高 a メートルの A 地点から，標高 b メートルの B 山の山頂を見るとき，その仰角はいくつか？ただし，A，B 間の水平距離は x メートルとする。という問題は小学校でもあるんじゃないか？まあ，高校生までは三角比を習わないので，正確な答えを得られなくても，図に描いて分度器で測ればよい。

防衛省の担当お役人は，この図を描くことを自分でやらずに Google map かなんかを使って，分度器で測ったらしい。よくあることではあるが，山の高さを目で見てわかるようにするとき，水平距離と垂直の高さの縮尺は後者が誇張されている。

常識人は，このようなとき分度器なんか使わない（そもそも，分度器なんて持ってますか？）。

tan(仰角) = (b-a) / x なのだから，仰角 = atan( (b-a) / x )。角度を度で表すなら，仰角 = atan( (b-a) / x ) * 180 / pi

最近，データの集計ミスの事案も何件も出ているけど，担当のお役人があまりにも文系過ぎるのではないか？それとも，いわゆるゆとり世代がお役人になっているのか？

きっと，この様子では，現地で見てみると，実際には標高が低い山（あるいは高いビル）が手前にあって，目的の山が見えません！ということもあるんじゃないか？（笑い事じゃない）。

Brunner-Munzel 検定（モンテカルロ法）

Brunner-Munzel 検定 ッテ，もはや，参照する価値もない記事です。パチもん検定です。

 

いくら，並べ替え検定が正確（？）といっても，サンプルサイズによる縛りがあるので，いつでも使えるわけではない。

そんなとき，サンプルサイズが大きい場合でも，任意の並べ替えに基づく検定を何回も行い，そのときの検定統計量と観察データにおける検定統計量の関係から完全な並べ替え検定の近似値を求めるのがモンテカルロ法（まあ，円周率を求めるときに，[0,1) の一様乱数を2組発生させて，原点からの距離が 1 未満になる確率を求めるようなことと同じ）

対象とする検定統計量は lawstat パッケージの brunner.munzel.test 関数の戻り値 estimate つまり P(X<Y)+.5*P(X=Y)

しつこいけど，t 値を対象としてはいけない。

bm.perm = function(x, y, trials=5000) {
  s = c(brunner.munzel.test(x, y)$estimate, brunner.munzel.test(y, x)$estimate)
  s.lo = min(s) + .Machine$double.eps
  s.hi = max(s) - .Machine$double.eps
  xy = c(x, y)
  nx = length(x)
  ny = length(y)
  nxy = nx + ny
  results = replicate(trials, {
    xy = sample(xy, nxy)
    brunner.munzel.test(xy[1:nx], xy[nx+1:ny])$estimate
  })
  mean(results <= s.lo | s.hi <= results)
}

これを使って，モンテカルロ法による p 値の推定値と，brunner.munzel.test の戻り値の p.value を比較する。

状況としては，第 2 群のサンプルサイズが 第 1 群のサンプルサイズの 2 倍 で 標準偏差も 2 倍

m = 200
p1 = p2 = numeric(m)
nx = 100; meanx = 50; sdx = 10
ny = 200; meany = 50; sdy = 20
for (i in 1:m) {
  cat(i, " ")
  x = rnorm(nx, meanx, sdx)
  y = rnorm(ny, meany, sdy)
  p1[i] = brunner.munzel.test(x, y)$p.value
  #t.test(x, y)$p.value
  p2[i] = bm.perm(x, y)
}
plot(p1, p2, xlab="brunner.munzel.test", ylab="permutation test", pch=19, cex=0.7, col="#66000020")
abline(0, 1, col=2)
abline(v=0.05, h=0.05, col=3, lty=2)

結果は以下の図のようになる。brunner.munzel.test()$p.value は，モンテカルロ法による p 値より小さめである。brunner.munzel.test()$p.value ＜ 0.05 となるのは 200 回中 9 回であるが，モンテカルロ法の場合には 200 回中 4  回だった。

trials=10000, m=1000 でやってみたら，同様な傾向であった。brunner.munzel.test()$p.value ＜ 0.05 となるのは 1000 回中 59 回であるが，モンテカルロ法の場合には 1000 回中 40  回だった。

cor.test の method="kendall" で

> set.seed(1234567)
> n = 200
> x = rnorm(n)
> y = rnorm(n)

n が50 以上のとき，exact を指定しない場合は正規近似が行われる

> cor.test(x, y, method="kendall")

    Kendall's rank correlation tau

data:  x and y
z = 1.2089, p-value = 0.2267
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau
0.05748744

exact = FALSE を指定しても問題はない。

> cor.test(x, y, exact=FALSE, method="kendall")

    Kendall's rank correlation tau

data:  x and y
z = 1.2089, p-value = 0.2267
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau
0.05748744

しかし，exact=TRUE を指定すると，暴走する

> cor.test(x, y, exact=TRUE, method="kendall")

    Kendall's rank correlation tau

data:  x and y
T = 10522, p-value = NA
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau
0.05748744

プログラムでは，

cor.test(x, y, ...., exact=NULL, ...)

  :

                if (is.null(exact))
                  exact <- (n < 50)

となっているので，exact が 50 以上であるにもかかわらず，exact=TRUE を指定すると歯止めが掛からず，プログラムが想定しない事態が生じる。

まあ，ほとんどの人はそんなことをしないのだけど，実装はまずいというしかない。