LiosK-free Blog

• 2010-06-03
• カテゴリ: Client Side
• タグ: JavaScript Tips アルゴリズム

前回の記事で予告したとおり、今回はJavaScriptのMath.random()で生成できる乱数の精度の話。

前回の記事で、JavaScriptでは2^53未満の正整数を扱うことができるということがわかったから、今回の記事では2^53未満のランダムな正整数を生成してみる。

具体的には↓のようなコード。

var ub = Math.pow(2, 53), list = [];
for (var i = 0; i < 16; i++) {
  list[i] = Math.floor(Math.random() * ub).toString(2);
  while (list[i].length < 53) { list[i] = "0" + list[i]; }  // padding
}
window.alert(list.join("\n"));

実行する

2^53未満の乱整数を2進法表現で表示する。手許にあった各ブラウザで実行してみた結果が↓。

• Internet Explorer 8
  
• Firefox 3.6
  
• Opera 10.5
  
• Safari 4
  
• Google Chrome 4
  

Webkit系の2つ (Safari, Chrome) がおかしいということがすぐわかると思う。Webkit系のブラウザでは乱数の精度が53ビット未満だから、二進法表現にしたときに末尾にゼロが並んでしまうのだ。

乱数の精度を推測する

Math.random()がどの程度の精度を持っているかを調べるためには、↓のようなコードを実行してみるとよい。

var list = [];
for (var i = 0; i < 16; i++) {
  list[i] = Math.random().toString(2);
}
window.alert(list.join("\n"));

実行する

小数点以下の桁数を数えてみて、その最大値が乱数の精度になる。

SafariやChromeで数えてみると最大でも30桁しかないことがわかる。つまり、Webkit系ブラウザの乱数は30ビットの精度しかない。

Safari, Chromeでも53ビット乱数を生成する

ということで、乱数精度が30ビットしかないSafariやChromeでも53ビットの乱整数を生成できるような関数を書いてみた。

/**
 * Returns an unsigned x-bit random integer.
 *
 * @param int x A positive integer ranging from 0 to 53, inclusive.
 * @return int An unsigned x-bit random integer (0 <= getRandomInt(x) < 2^x).
 */
function getRandomInt(x) {
  if (x <   0) return NaN;
  if (x <= 30) return (0 | Math.random() * (1 <<      x));
  if (x <= 53) return (0 | Math.random() * (1 <<     30))
                    + (0 | Math.random() * (1 << x - 30)) * (1 << 30);
  return NaN;
}

↓のように使えば53ビット精度の乱整数を得られる。

getRandomInt(53);

仕組みは簡単で、例えば引数として53を指定した場合、30ビットと23ビットの乱整数を生成し、後者を30ビット分左シフトした後に両者を合成しているだけ。

この関数を使ってSafariで53ビット乱数を生成してみる。

var list = [];
for (var i = 0; i < 16; i++) {
  list[i] = getRandomInt(53).toString(2);
  while (list[i].length < 53) { list[i] = "0" + list[i]; }  // padding
}
window.alert(list.join("\n"));

実行する

ぱっと見た感じではうまくいってそう。

まとめ

Math.random()で生成できる乱数の精度はブラウザ依存。最新バージョンでもWebkit系のブラウザは精度が低めなので、JavaScriptで扱える整数の範囲いっぱいいっぱいの乱数を使いたければ、↑で用意したような関数を使わなければいけない。

ちなみに、最新版のOperaでは問題なく53ビットの乱数を生成できているけれど、以前のバージョンのOperaでは精度が低かったような気がする。他にも乱数の精度が低いブラウザがあったりするんだろうか？

ってか、万が一乱数の精度が30ビット未満のブラウザがあったりすると、↑で用意した関数でも正しく乱数を生成できなくなるので、ちょっと困る。

蛇足1

関数を書いたので、とりあえずテストもしてみた。ここから先は専門的だし長いので興味ない人は読まなくてok。

まずはJavaScriptの関数として、仕様通りの値を返すかどうかテストする。とりあえずいくつか乱数を生成してみて、整数が生成されているか、値が期待した範囲に収まっているか、全てのビットが使用されているか、などをチェックする。

function getRandomIntTestA(put) {
  var N = 1024;

  if (!isNaN(getRandomInt(-1))) {
    put("NG\tgetRandomInt(-1) didn't return NaN");
  }
  if (!isNaN(getRandomInt(54))) {
    put("NG\tgetRandomInt(54) didn't return NaN");
  }

  var xs = [16, 32, 48, 53];
  while (xs.length) {
    var x = xs.shift(), ub = Math.pow(2, x);
    var bitsAppeared = new Array(x);

    for (var i = 0; i < N; i++) {
      var y = getRandomInt(x);
      if (!isFinite(y)) {
        put("NG\t" + "getRandomInt(" + x + ") didn't return a finite number: " + y);
      } else if (y < 0 || ub <= y) {
        put("NG\t" + y + " is out of the range of getRandomInt(" + x + ")");
      }
      if (y != parseInt(y)) {
        put("NG\t" + "getRandomInt(" + x + ") didn't return an integer: " + y);
      }

      var b = y.toString(2), c = x - b.length, z = "0";
      for (; c > 0; c >>>= 1, z += z) { if (c & 1) { b = z + b; } }

      for (var j = 0; j < x; j++) {
        if (b.charAt(j) == "1") {
          bitsAppeared[j] = true;
        }
      }
    }

    for (var i = 0; i < x; i++) {
      if (!bitsAppeared[i]) {
        put("NG\tBit " + i + " was never set to 1 during " + N + " times of getRandomInt(" + x + ") call.");
      }
    }
  }

  put("End of tests.");
}

実行する

NGが出ずにEnd of tests.まで行けばとりあえずは一安心。

蛇足2

↑のテストでは乱数の一様性などは全くテストされていないので、Rを使ってそれらを簡単にテストしてみる。

こちらのサイトなどを見てみると擬似乱数の検定方法なんかが載っているので、お絵描きしながら適当に検定などしてみる。

等確率性

乱数は各数値の出現確率が概ね等しくて、指定した範囲内で一様に分布するのが望ましい。カイ2乗検定を使って生成した乱数が一様に分布しているか確かめてみる。

↓のコードでは2000個ほどの乱数を16個のブロックに分割して一様性を検定する。さらに、これを100回ほど繰り返して、有意確率5%の範囲に入らない回数を数えている。一様に分布していればこの回数が95回(95%)ぐらいになるはず。

function getRandomIntTestB(put) {
  var N = 2048, ntrials = 100, nblocks = 16;

  put("ps <- c()");
  put("breaks <- seq(0, 2^53, length.out=" + (1 + nblocks) + ")");
  while (ntrials--) {
    var xs = new Array(N);
    for (var i = 0; i < N; i++) {
      xs[i] = getRandomInt(53);
    }
    put("xs <- c(\n" + xs.join(",\n") + "\n)");
    put("tab <- table(cut(xs, breaks))");
    put("ps <- c(ps, chisq.test(tab)$p.value)");
  }
  put("print(length(ps[0.05 < ps]))");
  put("print(chisq.test(tab))");
  put("barplot(tab)");
}

↓Safariでテストしてみた結果の例。最後の1回分はカイ2乗検定の結果とヒストグラムを出力している。

[1] 97

        Chi-squared test for given probabilities

data:  tab
X-squared = 12.2969, df = 15, p-value = 0.6564

この例では100回実行して97回は有意確率5%の範囲外にあった。そんなに悪くない数字。ヒストグラムを見てもそこそこ一様に分布している様子。

ちなみに、ブロックサイズを16ではなく32や64にしても概ね同じ結果がでる。

無規則性

隣あった乱数が何らかの相関を持っていたりするのは好ましくない。ということで、元の乱数列とそれをわずかにずらした乱数列の相関係数を検定し、相関がないことを確かめてみる。

↓のコードでは長さ512の乱数列とその数列を右に2つずらした乱数列を用意して、両者の相関係数を検定する。↑のテストと同様に100回実行してみて、5%有意の範囲外に収まる回数を数える。

function getRandomIntTestC(put) {
  var N = 512, ntrials = 100;

  put("gap <- 2");
  put("ps <- c()");
  put("indice <- c(seq(1 + gap, " + N + "), seq(1, gap))");
  while (ntrials--) {
    var xs = new Array(N);
    for (var i = 0; i < N; i++) {
      xs[i] = getRandomInt(53);
    }
    put("xs <- c(\n" + xs.join(",\n") + "\n)");
    put("ps <- c(ps, cor.test(xs, xs[indice])$p.value)")
  }
  put("print(length(ps[0.05 < ps]))");
  put("print(cor.test(xs, xs[indice]))");
  put("plot(data.frame(xs=xs, ys=xs[indice]))");
}

↓IE8で実行してみた結果。↑のテストと同様に、最後の1回分に関しては検定結果と散布図を表示している。

[1] 95

        Pearson's product-moment correlation

data:  xs and xs[indice]
t = 0.4967, df = 510, p-value = 0.6196
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.06478947  0.10843932
sample estimates:
       cor
0.02198998

この例では100回中95回が5%有意の範囲外にあったので問題なし。散布図をみても散らばっている感じ。

ちなみに、ずらす量を2ではなくて1や4などにしてもほぼ同じ結果になる。

中心極限定理

中心極限定理は、任意の分布からランダムに選んだサンプルの平均値の分布は正規分布に従うという確率論の定理。乱数列は一様分布からランダムに選んだサンプルにあたるので、当然その平均値の分布は正規分布に従うはず。ということで、確かめてみる。

↓のコードは32個の乱数の平均値を2000個ほど集め、その分布が正規分布に従っているかをシャピロ-ウィルク検定で検定する。例によって100回実行する。

function getRandomIntTestD(put) {
  var N = 2048, ntrials = 100, blocksize = 32;

  put("ps <- c()");
  while (ntrials--) {
    var xs = new Array(N);
    for (var i = 0; i < N; i++) {
      xs[i] = 0;
      for (var j = 0; j < blocksize; j++) {
        xs[i] += getRandomInt(53) / blocksize;
      }
    }
    put("xs <- c(\n" + xs.join(",\n") + "\n)");
    put("ps <- c(ps, shapiro.test(xs)$p.value)")
  }
  put("print(length(ps[0.05 < ps]))");
  put("print(shapiro.test(xs))");
  put("print(summary(xs))");
  put("hist(xs)")
}

↓Operaで実行してみた結果の例。

[1] 94

        Shapiro-Wilk normality test

data:  xs
W = 0.9994, p-value = 0.7381

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
3.012e+15 4.180e+15 4.502e+15 4.504e+15 4.807e+15 6.032e+15

この例では100回中94回が5%水準の外にあったので問題なさそう。ヒストグラムを見てみてもきれいに正規分布している。

結論

↑のテストを見た感じだとちゃんと乱数として機能しているみたい。

スピードもそれなりに出るので、精度の高い乱数が必要なときはぜひ。

関連記事

• [メモ] Xubuntuでおひとり様サーバーを立てるまで 1
• 乱数の精度を調べてたらOpera Miniのバグを見つけた
• JavaScriptの乱数の精度の話
• JavaScript (ECMAScript) の整数の精度と乱数の精度について
• [JavaScript] prototypeに直接代入しちゃうのってダメじゃなかったっけ？

• コメント (0)
• トラックバック (1)
• 
• Tweet

http://liosk.blog103.fc2.com/tb.php/198-4a4ab930

乱数の精度を調べてたらOpera Miniのバグを見つけた

前回の記事に引き続いてJavaScriptの乱数の精度を調べてたら、iPhone用のOpera Miniのバグを見つけた。 NumberオブジェクトのtoStringメソッドがうまく働...

• 2010-06-05
• 発信元: LiosK-free Blog