単位球体内の一様分布： 4次元

単位球体内の一様分布を生成する方法を見ていくシリーズ（目次）& 倭式極座標を見ていくシリーズ（目次）。　今回は4次元。　3次元までの単位球体内（通常の円盤内、球体内）の一様分布はその生成方法はよく知られていますが、4次元以上では通常使われる極座標（『極座標のヤコビ行列とヤコビアン： n次元』を参照）では3次元までのように簡単には計算できません。　その理由は、この4次元極座標では体積要素が

　　 ${ \displaystyle\begin{align*} r^3\sin^2\theta_1\sin\theta_2 dr d\theta_1 d\theta_2 d\theta_3 \end{align*}}$

となりますが、この体積要素が、うまく変数を導入して

　　 ${ \displaystyle\begin{align*} dRd\Theta_1 d\Theta_2 d\Theta_3 \end{align*}}$

のように書けないことが原因です（今の場合、 ${ \Theta_1 }$ に当たるものがない*1）。　で、通常は4次元以上の球体内に一様分布する点を生成するには Box-Muller 法（wikipedia:ボックス＝ミュラー法）を使うんですが、以前の記事に書いた倭式極座標というのを使うと、2, 3次元の場合と同じように4次元以上の球体内でも一様に分布する点を生成することができます。　以下でその方法を見ていきましょう。

4次元倭式極座標

4次元の倭式極座標は以下のように与えられるのでした（『4次元倭式極座標』を参照）：

　　 ${ \displaystyle\begin{align*} &\begin{cases} x_1 = r \sin\theta \sin\varphi_1 \\ x_2 = r \sin\theta \cos\varphi_1 \\ x_3 = r \cos\theta \sin\varphi_2 \\ x_4 = r \cos\theta \sin\varphi_2 \\ \end{cases} & \begin{pmatrix} 0 \le r < \infty \\ 0 \le \theta \le \tfrac{\pi}{2} \\ 0 \le \varphi_* < 2\pi \end{pmatrix} \end{align*}}$

単位球内での積分は以下のようになります

　　 ${ \displaystyle\begin{align*} \int_0^1 r^3dr \int_0^{\pi/2} \sin\theta\cos\theta d\theta \int_0^{2\pi} d\varphi_1 \int_0^{2\pi} d\varphi_2 \end{align*}}$

一様分布する変数

一様分布する変数への変数変換は以下の通り：

　　 ${ \displaystyle\begin{align*} (r,\, \theta,\, \varphi_1,\, \varphi_2) \rightarrow (R,\, \Theta,\, \varphi_1,\, \varphi_2 ) = (\tfrac{1}{4} r^4,\, \tfrac{1}{2}\sin^2\theta,\, \varphi_1,\, \varphi_2) \end{align*}}$

このとき、単位球体内の積分は以下のようになります：

　　 ${ \displaystyle\begin{align*} \int_0^{1/4} dR \int_0^{1/2} d\Theta \int_0^{2\pi} d\varphi_1 \int_0^{2\pi} d\varphi_2 \end{align*}}$

各変数について被積分関数を1にすることができれば、その変数について一様分布を生成し、そこからの代数計算で望みの4次元球体内に一様分布する点の座標が計算できます。　次節は生成した一様分布から4次元球体内の点の座標を計算するための関係式を導きます。

直交座標の計算

極座標の変数を一様分布する変数で表すと

　　 ${ \displaystyle\begin{align*} (r,\, \theta,\, \varphi_1,\, \varphi_2) = (\sqrt[4]{4R},\, \sin^{-1}\sqrt{2\Theta},\, \varphi_1,\, \varphi_2) \end{align*}}$

となるので、直交座標は以下のようになります：

　　 ${ \displaystyle\begin{align*} \begin{cases} x_1 = \sqrt[4]{4R} \sqrt{2\Theta} \sin\varphi_1 \\ x_2 = \sqrt[4]{4R} \sqrt{2\Theta} \cos\varphi_1 \\ x_3 = \sqrt[4]{4R} \sqrt{1-2\Theta} \sin\varphi_2 \\ x_4 =\sqrt[4]{4R} \sqrt{1-2\Theta} \cos\varphi_2 \end{cases} \end{align*}}$

実際には、 ${ 4R,\, 2\Theta }$ をそれぞれ改めて ${ R,\,\Theta}$ とおき、以下の形で使う方が効率的：

　　 ${ \displaystyle\begin{align*} &\begin{cases} x_1 = \sqrt[4]{\tilde{R}} \sqrt{\tilde{\Theta}} \sin\tilde{\varphi}_1 \\ x_2 = \sqrt[4]{\tilde{R}} \sqrt{\tilde{\Theta}} \cos\tilde{\varphi}_1 \\ x_3 = \sqrt[4]{\tilde{R}} \sqrt{1-\tilde{\Theta}} \sin\tilde{\varphi}_2 \\ x_4 =\sqrt[4]{\tilde{R}} \sqrt{1-\tilde{\Theta}} \cos\tilde{\varphi}_2 \end{cases} & \begin{pmatrix} 0 \le r \le 1 \\ 0 \le \Theta \le 1 \\ 0 \le \varphi_* < 2\pi \end{pmatrix} \end{align*}}$

変数 ${ X \in [a,\,b] }$ （もしくは区間 ${ [a,\,b) }$ ）に対して、 ${ \tilde{X} }$ を区間 ${ [a,\,b] }$ 間に一様分布する変数とします。

Java コード

上式を使って4次元単位球体内に一様分布する点の座標を生成するコードは以下のようになります：

import static java.lang.Math.*;

public class UniformDistributionB4{

    public static void main(String... args){

        for(int i = 0; i < 10000; i++){
            double r = pow(random(), 0.25);
            double theta = random();
            double s = sqrt(theta);
            double c = sqrt(1- theta);
            double phi1 = random() * 2 * PI;
            double phi2 = random() * 2 * PI;
            
            double x1 = r * s * sin(phi1);
            double x2 = r * s * cos(phi1);
            double x3 = r * c * sin(phi2);
            double x4 = r * c * cos(phi2);
            
            System.out.println(x1+" "+x2+" "+x3+" "+x4);
        }
    }
}

UniformDistributionB4 は ${ B_4 }$ （4次元球 4-dimensional ball）中の一様分布 (uniform distribution) の意。　標準出力に書き出さずに配列やリストとして返したり、4次元座標を生成する部分をメソッドとして抽出したりといろいろ改善の余地はありますが、そのあたりは後日に n 次元でやれればと。　重要なのは、Box-Muller 法のように条件に合わない場合に再試行したりしなくていいってところです。　これは高次元になればなるほどパフォーマンスに差が出てきます。　倭式極座標スゲーよ。