3 次元空間の単位ベクトルを表す乱数

新しい記事「3 次元空間の単位ベクトルを表す乱数（その4）」を参照

http://sssslide.com/www.slideshare.net/T2C_/kadai9-report-forpublic
ピアソン問題（課題9）であるが...

私がやってみると，結果がだいぶ違う。
どうも，作者が作ったシミュレーションデータに問題があるのではないかと思う。
データの要点は，3次元空間におけるランダムな方向を向いている単位ベクトルの x, y, z 座標を作るところにある。
作者のデータの作成法は，以下のようにまとめられる。
1. 3 個の一様乱数 [0, 1) ξ1, ξ2, ξ3 を生成する
2. u = 2*ξ1-1, v = 2*ξ2-1, w = ξ3 とする
3. d = u^2 + v^2 + w^2 とする
4. d ≦ 1 の場合に
   x = 2*u*w/d
   y = 2*v*w/d
   z = (w^2-u^2-v^2)/d とする
となっている。実際には 4. は，ベクトルの標準化をしているのだから，「d ≦ 1 の場合に」という条件は不要
このように作られたベクトルの座標を表す x, y, z は，歩数 n ごとに決まり，それぞれの合計が最終地点の X, Y, Z 座標になる。X, Y, Z は互いに無相関である。しかし，ベクトル X, Y の平均値はほぼ 0 であるが，Z の平均値は「ほぼ -1」になってしまう。

trial = 100000
X = Y = Z = numeric(trial)
for (i in 1:trial) {
    x = y = z = 0
    for (j in 1:3) {
        repeat {
            xi1 = runif(1)
            xi2 = runif(1)
            xi3 = runif(1)
            u = 2 * xi1 - 1
            v = 2 * xi2 - 1
            w = xi3
            d = u ^ 2 + v ^ 2 + w ^ 2
            if (d <= 1) break
        }
        x = x + 2 * u * w / d
        y = y + 2 * v * w / d
        z = z + (w ^ 2 - u ^ 2 - v ^ 2) / d
    }
    X[i] = x
    Y[i] = y
    Z[i] = z
}
> mean(X)
[1] -8.390827e-05
> mean(Y)
[1] 0.002259335
> mean(Z)
[1] -0.9993304 # **** !!!!!!!!!!!!!

実はこの X, Y, Z は 3 次元の（連続）ランダムウォークの結果なので，X, Y, Z の平均値が 0 でないと，結果が歪んでしまう。
のではないかなと思う。
作者に質問したいけど，コメントを付ける手立てがない。

私の書いたプログラムは，3 次元角座標の角度 2 つをランダムに生成し，角座標から直交座標 x, y, z に戻す方法を使っている。

PearsonProblem = function(n = 3, dR = 0.1, trial = 1e+05) {
    angle = matrix(runif(2 * n * trial, min = 0, max = 2 * pi), 2 * n, trial)
    ans = apply(angle, 2, function(ang) {
        ph = ang[1:n]
        th = ang[-(1:n)]
        X = sum(cos(ph) * cos(th))
        Y = sum(cos(ph) * sin(th))
        Z = sum(sin(ph))
        sqrt(X^2 + Y^2 + Z^2)
    })
}
ans = PearsonProblem(n = 3, trial = 1e+06)
mean(ans)
hist(ans, nclass = 50)

3 歩，歩いたときの原点からの距離の分布は，下図のようになった。平均値は 1.617797