一様乱数から正規乱数を作る方法?


ソフトウェア関係の雑誌を読んでいたら、次のようなことが書かれていました。

「0以上1未満の一様乱数を12回足して6を引くと正規乱数になる」というテクニックを覚えておけば、一様乱数しか提供されていない開発環境で大いに役立つはずだ

これは中心極限定理のいい例ですね。実際に乱数を大量(下の絵では100万個)に発生させてそのヒストグラムを描いてみると、正規分布になっているように見えます(太線は正規分布)。これで十分な場合も多いでしょう。

しかし、正規性の検定をすると、正規分布であるという仮説を棄却したくなるようなp値が得られます。Mathematicaで試します(参考:DistributionFitTest)。

In[1]:= randomGaussian[] := Sum[RandomReal[], {12}] - 6

In[2]:= data = Table[randomGaussian[], {1000000}];

In[3]:= DistributionFitTest[data, NormalDistribution[0, 1]]

                  -8
Out[3]= 8.33717 10

この方法で生成される乱数の分布は、足す数(上の例では12)を増やせば正規分布に収束しますが、正規分布ではありません。ですから、「正規乱数になる」という表現はいけません。(例えば「randomGaussian[] := (Sum[RandomReal[], {24}] - 12)/Sqrt[2]」ならもっと正規分布に近づきますが、依然として正規分布ではありません。)

一様乱数しか提供されていない環境では、ボックス=ミュラー法などを使うのがいいと思います。Excelでは「=NORM.INV(RAND(),0,1)」が簡単です。

4874085601Pressほか『ニューメリカルレシピ・イン・シー』(技術評論社, 1993)にこんな記述がありました。

悪いrand()のため結果が疑わしいすべての科学論文を図書館の棚から消せば、各棚にこぶしが入るほどの隙間ができるであろう。(p.203)

「悪いrand()」ではなく「悪い乱数」を対象にすると、隙間はどのくらいになるのでしょう。

Excelを使っていると採点ミスで謝罪することになる?


事例が見つかるたびに騒ぎになっているようです。今回は、「8.2 – 7.2」の結果が1より小さくなってしまうというものです。最初に非難されたのは例によってExcelのようですが、他の環境でもこういう問題は起こります(C言語の例Rubyの例)。

#include <stdio.h>

int main() {
  double x = 8.2 - 7.2;
  if (x >= 1) printf("x >= 1");
  else printf("x < 1");
  return 0;
}
&#91;/c&#93;

<p>「コンピュータは数値を2進数としてしか表せなくて・・・」という説明をよく見ますが、この説明は間違いです。これでは、「コンピュータは正しい計算ができない」という誤解を生んでしまいます。</p>

<p>10進数の1, 2, 3, ...を2進数の1, 10, 11, ...で表すという具合に、数を2進数で直接表現しようとするとこういう問題が起きますが、別の表現方法を用いれば、このような問題は起きません。実際、Mathematicaは正しく計算します。<ins datetime="2011-08-16T15:48:38+00:00">(<a href="http://www.unfindable.net/umm/">UMM</a>でも確かめられます。)</ins></p>

<pre>
With[{x = 8.2 - 7.2}, If[x &gt;= 1, "x &gt;= 1", "x &lt; 1"]]

x &gt;= 1
</pre>

<p>コンパイルしても大丈夫です。<ins datetime="2011-08-16T15:48:38+00:00">(UMMでは動かないようです。)</ins></p>

<pre>
f = Compile[{a, b}, With[{x = a - b}, Print[If[x &gt;= 1, "x &gt;= 1", "x &lt; 1"]]]];
f[8.2, 7.2]

x &gt;= 1
</pre>

<p>C言語やRubyではどうすればいいのかというと、許容可能な誤差を決めておいて、その範囲内はよしとするのです。</p>

<p>さらに、この許容可能な誤差は、扱う数によって変えなければならないことに注意してください。以下の例では、許容可能な誤差を1より大きくしなければなりません(<a href="http://codepad.org/Fg4pP5H7">実行結果</a>)。</p>

[c]
#include <stdio.h>

int main() {
  double a = 9007199254740992;
  double b = a + 1;
  if (a == b) printf("a==b");
  return 0;
}

4891006269こういうことは、『Microsoft Visual C++入門』の3.2.11項「浮動小数点についての注意」にちゃんと書いてあったりします。ですから、入門を終えたプログラマはこういうことは知っています。

Excelのような、こうことは知らないであろう人が使うソフトウェアでは、こういう工夫をしなくても正しく計算できるようにしておいてほしいものです。フェルマーの最終定理の「反例」のような冗談を言っているうちはいいのですが、「x >= 0.6 なら合格」というようなことって、よくやられているのではないでしょうか(ボーダーの学生が間違って不合格になる危険があります)。間違いが公になるときには、「採点ミス」になるのでしょうか。

Excelのようなプログラミング


Mathematicaはバージョン6でDynamicというおもしろい機能が追加されました(動的言語とは関係ありません。もとから動的ですし)。いろんな解釈があると思いますが、Dynamicで何ができるかというと、Excelのようなプログラミングです。ここで考えているExcelの機能とは、複数のセルが関連しているところで、あるセルの値が変更されると、関連するセルの値が自動的に再計算される、というものです。

Dynamicによって、Mathematicaのプログラミングスタイルは大きく変わると思うのですが、残念ながら、今のところは「GUIを簡単に作れる」という応用が主になっています。

たとえばこんな感じです。

center[p_, q_, r_] :=
 With[{np = Normalize@p, nq = Normalize@q},
  r Normalize[np + nq]/Sin[ArcCos[np.nq]/2]]

Manipulate[
 Graphics[{
   Line@{o, p},
   Line@{o, q},
   Circle[o + center[p - o, q - o, r], r]
   }, PlotRange -> 2],
 {{o, {-1, 0}}, Locator},
 {{p, {0, 1}}, Locator},
 {{q, {0, -1}}, Locator},
 {{r, 1}, 0, 2},
 SaveDefinitions -> True]

Wolfram CDF Playerがインストールされていれば実際に動かすことができます(UMMでも動きます)。


「こういうことがCAD上でできたらおもしろい」と思ったのですが、AutoCADなんかですでに実現されていますね。

だから、Mathematicaを使う人たちは、もっと別の応用を考えなければいけないのです。

フェルマーの最終定理の「反例」(表計算ソフト編)


「A1^3+B1^3-C1^3」の結果が0なので、フェルマーの最終定理の「反例」になっています。Excel 2007, 2008だけでなく、OpenOffice Calc 3.0RC1, gnumeric 1.8.2, KPsread 1.6.3, Google Docs BETA等でもこうなります。

表計算ソフトの利用者の大部分は、コンピュータが初頭的な計算において間違いを犯すとは思ってはいないのではないでしょうか。

参考:Excel 2007に乗算のバグが発見される(解決済み)

Excelの計算ミスのようなことはMathematicaにはない?


Excelで「850*77.1」を計算すると100000になるように見えるというバグが話題になっている。「なる」ではなく「なるように見える」と書いたのは、結果をグラフなどで2次利用する際には問題ないから(つまり内部では正しい値を持っている)。もちろん、結果をそのまま使う場合には問題になる。

この話題がWolfram Blogで取り上げられていて(上の画像もこのブログへの直リンク)、似たような話として、40000.223と入力したつもりが40000.2229999999になってしまう(逆ではない)例が紹介されていた。

「テストで見つけるには相当な知識が必要だよ」とか、「ふつうの処理系はハードウェアに組み込まれた精度での計算しかしないが、Mathematicaは任意の精度で計算できるから大丈夫」とのこと。

「テスト」なんてそんなにありがたいものではないと私は思っているし、計算結果の信頼性はExcelよりMathematicaのほうが上だというのもそのとおりかもしれない。でも、それは程度の問題であって、Mathematicaだって完璧ではない。たとえば先日の問題で、「10^-4」と書くところを「0.0001」にしたりすると、Mathematicaだって正しい結果は返してくれない(UMMで確認できる)。

Select[Range@199,
 With[{f = Exp[Pi Sqrt@#]},
  Abs[f - Round@f] <= 0.0001 &]]

{37, 58, 96, 100, 103, 104, 117, 118, 119, 123, 124, 126, 128, 
132, 134, 138, 140, 143, 144, 146, 147, 148, 150, 152, 153, 
156, 160, 162, 164, 166, 167, 168, 172, 175, 178, 179, 180, 
181, 182, 184, 186, 187, 188, 189, 190, 192, 194, 195, 197, 
198, 199}