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

単位球体内の一様分布を生成する方法を見ていくシリーズ（目次）。　今回は3次元。　やることは2次元の場合と同じ。　ちなみに3次元の単位球体は普通の「球」です。

極座標

3次元の極座標は以下のように与えられるのでした：

　　 ${ \displaystyle\begin{align*} &\begin{cases} x = r \sin\theta \cos\varphi \\ y = r \sin\theta \sin\varphi \\ z = r\cos\theta \end{cases} & \begin{pmatrix} 0 \le r < \infty \\ 0 \le \theta \le \pi \\ 0 \le \varphi < 2\pi \end{pmatrix} \end{align*}}$

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

　　 ${ \displaystyle\begin{align*} \int_0^1 r^2dr \int_0^\pi \sin\theta d\theta \int_0^{2\pi} d\varphi \end{align*}}$

一様分布する変数

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

　　 ${ \displaystyle\begin{align*} (r,\, \theta,\, \varphi ) \rightarrow (R,\, \Theta,\, \varphi) = (\tfrac{1}{3} r^3,\, \cos\theta,\, \varphi ) \in [0,\, \tfrac{1}{3}] \times [-1,\, 1] \times [0,\,2\pi ) \end{align*}}$

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

　　 ${ \displaystyle\begin{align*} \int_0^{1/3} dR \int_{-1}^1 d\Theta \int_0^{2\pi} d\varphi \end{align*}}$

直交座標の計算

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

　　 ${ \displaystyle\begin{align*} (r,\, \theta,\, \varphi) = (\sqrt[3]{3R},\, \cos^{-1}\Theta,\, \varphi) \end{align*}}$

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

　　 ${ \displaystyle\begin{align*} \begin{cases} x = \sqrt[3]{3R} \sqrt{1-\Theta^2} \cos\varphi \\ y = \sqrt[3]{3R} \sqrt{1-\Theta^2} \sin\varphi \\ z = \sqrt[3]{3R} \Theta \end{cases} \end{align*}}$

これを用いて3次元球体内の一様乱数を生成するには、効率化のために ${ R }$ に3を吸収させて（定義域にも）、以下のようにします：

　　 ${ \displaystyle\begin{align*} \begin{cases} x = \sqrt[3]{\tilde{R}} \sqrt{1-\tilde{\Theta}^2} \cos\tilde{\varphi} \\ y = \sqrt[3]{\tilde{R}} \sqrt{1-\tilde{\Theta}^2} \sin\tilde{\varphi} \\ z = \sqrt[3]{\tilde{R}} \tilde{\Theta} \end{cases} \begin{pmatrix} 0 \le R \le 1 \\ -1 \le \Theta \le 1 \\ 0 \le \varphi < 2\pi \end{pmatrix} \end{align*}}$

ただし、変数 ${ X \in [a,\,b] }$ に対して ${ \tilde{X} }$ を区間 ${ [a,\,b] }$ 間に一様分布する変数とします。

Java コード

実際に Java で3次元球体内の一様乱数を生成してみましょう。　Java コードは以下のようになります：

import static java.lang.Math.*;

public class UniformDistributionB3{

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

        for(int i = 0; i < 30000; i++){
            double r = cbrt(random());
            double cosTheta = random() * 2 - 1;
            double sinTheta = sqrt(1- cosTheta*cosTheta);
            double phi = random() * 2 * PI;
            
            double x = r * sinTheta * cos(phi);
            double y = r * sinTheta * sin(phi);
            double z = r * cosTheta;
            
            System.out.println(x+" "+y+" "+z);
        }
    }
}