ぎるばーとのノート

もっともっと遠くへ行きたい　空が広く見える場所まで

三角乱数

数学・統計学

　マニアックな乱数記事を徐々に増やしていくキャンペーン。

　前回・前々回と、なぜ数学関数の使用を回避したいかを書いてなかったので、その補足から、、
　プログラムの数学関数、具体的にはサイン、コサイン、対数、平方根などですが、処理が重いことが多いです。除法も加・減・乗法と比べれば若干重いとはいえ、これの回避は難しい……。
　数学関数を使わないかわりに、単位円内の一様乱数を作るための繰返しが入りました。その繰返し回数を見積もってみます。
　円と外接正方形の面積の比は\pi/4 \approx 0.785 \ldotsで、円内に点が入ったときに脱出。その逆数が繰返し回数の平均になります。幾何分布なので数回程度かかる場合もありますが、平均は1回ちょっとです。
　ちなみに、実際にサインやコサインを使うより速いかは環境によりけりかも……？

　ここからが今日の本文。
　三角分布をPERTで使うという話は本当なのかどうか知りません。

Triangular distribution - Wikipedia
　三角乱数、平方根を使う方法がWikipediaに書いてあります。が、もっといいやり方（感動的！）があるので紹介します。

A new method to simulate the triangular distribution - ScienceDirect
u \sim \text{Uniform}(0,1),\ v \sim \text{Uniform}(0,1)として、
(1-c)\text{MIN}(u,v) + c\text{MAX}(u,v)
　これで、三角乱数（a=0, b=1, cは最頻値で0以上1以下）になります。MINとMAXは比較だけなので、平方根の計算よりずっと軽い！

double randomTriangular(double c) {
    double u = Math.random(), v = Math.random();
    
    return (1.0 - c) * Math.min(u, v) + c * Math.max(u, v);
}

　この方法が紹介されているのをあまり見ませんが、美しいアルゴリズムの一つだと思います。

コーシー乱数

数学・統計学

　連日の更新なんていつぶり……？

　コーシー分布について。ちなみにコーシーと孔子は別人です。（豆知識
　コーシー分布は自由度1のt分布、平均も分散も発散するたちが悪い分布として名を知られてますね。
　和に対して安定（和の分布もコーシー分布）なのは正規分布と同じです。が、正規分布はスケールの二乗（分散）が加法的で、コーシー分布はスケールが加法的になります。

　分散1の独立な正規乱数の和は分散2の正規乱数になって、スケール（標準偏差）で見ると√2倍。さらに標本平均の分散は、
\text{Var}(\bar X) = \text{Var}(\frac{\sum X_i}{n}) = \frac{1}{n^2}\text{Var}(\sum X_i) = \frac{n\sigma^2}{n^2} = \frac{\sigma^2}{n}
　標本サイズが大きくなることで標本平均の分散が小さくなり、母平均の推定の信頼性が高まるのでした。

　片や、スケール1の独立なコーシー乱数の和はスケール2のコーシー乱数になります。n個の和をとって、スケールnのコーシー乱数。nで割るとスケール1のコーシー乱数。と、標本平均で元の分布に戻ってしまいました。母確率分布がコーシー分布のとき、標本サイズを大きくしても母中心の推定の信頼性は高まりません。

Cauchy distribution - Wikipedia
　コーシー乱数の生成には、逆関数法とか標準正規乱数の比をとるとかいろいろ見つかりますが、ここでは数学関数を使わない方法を、、

double randomCauchy() {
    double u, v;
    do {
        u = 2.0 * Math.random() - 1.0;
        v = Math.random();
    } while (u * u + v * v >= 1.0 || v == 0.0);
    
    return u / v;
}

　単位円内（上半分）からランダムに選んだ点で、X座標をY座標で割るとコーシー分布(0,1)の乱数が得られます。説明は割愛、Box–Mullerとかで検索するといいんじゃないかな……。

逆正弦乱数

数学・統計学

physnotes.jp

　逆正弦分布（アークサイン分布）という確率分布があるのを、高校物理の備忘録さんの記事で知りました。記事の後半の逆正弦法則の話もとても興味深いです。

　振り子の速度は最低点で最大、両端では0になります。なので、最低点付近にいる時間は短く、両端付近にいる時間が長くなり、ランダムなタイミングでチラッと見たときの位置は逆正弦分布に従います。

Arcsine distribution - Wikipedia, the free encyclopedia
　逆正弦分布(0,1)は、ベータ分布の特別な場合（α=β=0.5）の別名です。この分布の乱数を生成する方法（用途はまったく不明、、）を書き記しておきます。

　Wikipediaにある通り、一様乱数U \sim \text{Uniform}(-\pi,\pi)から、(\sin U)^2とすればいいわけですが、サイン関数を使わずに計算する方法もあります。正規乱数でおなじみのMarsaglia法を応用します。

double randomArcsine() {
    double u, v, s;
    do {
        u = Math.random();
        v = Math.random();
        s = u * u + v * v;
    } while (s >= 1.0 || s == 0.0);
    
    return (u * u) / s;
}

　第1象限で単位円内の点をランダムに選んで、原点からの距離で割ると円周上に分布します。座標を一つとって二乗すると目的の乱数が得られますね。正規乱数の場合のように二つ独立に作れるわけではないのに注意！　(\sin x)^2 + (\cos x)^2 = 1なので……。

Javaの知られざる構文〜幻のfor-while文

コンピュータ Java

【この記事はジョーク記事です！】

　この件に関連して……。

　あまり知られていないのだけど、初期のJavaに、for文とdo-while文を組み合わせたfor-while文というのがありました。いつ頃からか、言語規定に載らなくなりましたが、最新版でもコンパイル・実行が可能です。

// for-while文の例
int i;
for (i = 0; i < 5; i++) {
    System.out.println(i);
} while (i < 5);

　iはメモリに記憶されますが、ループ判定のタイミングでCPU内のレジスタに読み込まれます。

　宇宙空間や過酷な放射線環境下で使用する機器では、メモリは放射線対策していることが多いですが、CPUは特に何でもないものを使ったりするそうです。そこで、CPU内のビット反転に対する耐性を高めたい場合に使用できる構文ということだったのですが……。
　うーん、可読性や保守性が悪すぎたのか……。幻の構文となってしまいました。





　……信じましたか？(^^;

　「種明かし」しましょう。上のコードは、for文の後ろに「中身が空っぽ」のwhile文が続いているだけです。つまり、

// for-while文の例
int i;
// for文
for (i = 0; i < 5; i++) {
    System.out.println(i);
}
// while文
while (i < 5)
    ;

　まあ、いくら何でもそんな変な構文あるわけないですね！（信じてしまった方は、悪意の人間に騙されないよう気をつけて……。）

Stirlingの近似式

数学・統計学

　Stirlingの近似式は、大きな数の階乗を見積もることができます。
\ln n! \simeq n \ln n - n \tag{1}
　この式は「両辺の比が1に収束する」という意味です。もっと強く、「両辺の差が0に収束する」近似式は次の通り。
\ln n! \simeq n \ln n - n + \frac{1}{2} \ln 2\pi n \tag{2}

　ところで、階乗の定義域を実数に広げたものはガンマ関数で、
\Gamma(x) = \int_0^\infty t^{x-1}e^{-t} dt
\Gamma(n + 1) = n! = 1 \times 2 \times 3 \times \cdots \times n
　ガンマ関数の対数をとって微分するとディガンマ関数になります。
\psi(x) = \frac{d}{dx} \ln \Gamma(x) = \frac{\Gamma'(x)}{\Gamma(x)}
　ディガンマ関数は調和級数の拡張版といったところ。nが非負整数なら、
\psi(n + 1) = -\gamma + H_n = -\gamma +1 + \frac{1}{2} + \frac{1}{3} + \cdots + \frac{1}{n}
　式中の\gammaは有名なEulerの定数です。これも語りたいことはあるけれど、また別の機会に……。（たぶん

　話を元に戻しまして、Stirlingの式(1)の両辺をnで微分してみます。
　左辺は、
\frac{d}{dn} \ln n! = \frac{d}{dn} \ln \Gamma(n + 1) = \psi(n + 1)
　右辺は、
\frac{d}{dn} (n \ln n - n) = \ln n
　ディガンマ関数と対数関数が出てきました！　この二つは興味深い関係がたくさんありますが、発散する速さに注目すると、
\lim_{n \to \infty} \{\psi(n + 1) - \ln n\} = 0
　つまり同じ速さで発散します。（1を任意の実数\alphaに変えても同じ。）

　下の図はディガンマ関数と対数関数のグラフです。
f:id:zwxadz:20151001214330p:plain

　まとめると、Stirlingの式が示しているのは「階乗の対数（ディガンマ関数の積分）を対数関数の積分で近似できるよ！」ということになります。