「統計学関連なんでもあり」の過去ログ--- 043

No.11404　Euler の sum 関数　　【青木繁伸】　2009/12/03(Thu) 18:19

数値計算ネタなので，統計学とは直接関係ないけど，R の計算は信頼できると思われるのでここに書いておきます。

中澤さんのページ http://phi.med.gunma-u.ac.jp/index.html で，Euler という，R に似た言語で，ζ(2) = sum(1/(1:10000000)^2) と円周率の2乗を6で割った値=pi^2/6 が等しい結果になることを示していたので，私も Euler をダウンロードしてやってみました。Mac 版がないのが残念。で，ちょっと，おもしろい結果を見つけたので報告。

Euler での結果

>longestformat
20 12
>pi^2/6
1.6449340668482264
>sum(1/(1:10000000)^2)
1.6449339668472596
>a=0; for i=1 to 10000000; a=a+1/i^2; end; a # 大きい方から足す
1.6449339668472596
>a=0; for i=10000000 to 1 step -1; a=a+1/i^2; end; a # 小さい方から足す
1.6449339668482315

sum(1/(1:10000000)^2) は，大きい方から（指定された方から）足した結果と同じ。つまり，指定順に足していっている。
R での結果

> options(digits=20)
> pi^2/6
[1] 1.644934066848226
> sum(1/(1:10000000)^2)
[1] 1.64493396684823
> sum(1/(10000000:1)^2)
[1] 1.64493396684823
> a <- 0; for (i in 1/(1:10000000)^2) a <- a+i; print(a) # 大きい方から足す
[1] 1.644933966847318
> a <- 0; for (i in 1/(10000000:1)^2) a <- a+i; print(a) # 小さい方から足す
[1] 1.64493396684823
>a <- 0; for (i in sample(1/(10000000:1)^2)) a <- a+i; print(a) # ランダムに足す
[1] 1.644933966847353

R の sum(1/(1:10000000)^2) は，引数の指定順にかかわらず，常に小さい方から足した結果と同じになるようだ。積み残しの危険があるので，足し算は常に小さい順に足すのがよい。もっと厳密に言えば，小さい2項を足した結果と次に小さいものを足すのではなく，毎回最小の2項を足し算するのがよい。
R は，そのあたりまで考えて計算を行っているのかなと。

No.11408　Re: Euler の sum 関数　　【青木繁伸】　2009/12/03(Thu) 21:59

多倍長演算プログラムを書いていて，どうしても合わないと思って，よく見たら，
> pi^2/6
[1] 1.644934066848226
> sum(1/(1:10000000)^2)
[1] 1.64493396684823
で，小数点以下6,7桁があっていない。それにしても，ここだけが違うとは？？
> sum(1/(1:100000000)^2)
[1] 1.644934056848227
なので，そうとう収束が遅いということだろう。