Java で球面上・球体内の一様分布を生成する (2)：単位球面上の一様分布

今回は「単位球面上の一様分布」を生成するコードを考えます（一覧）。

setPointOnSphere(double x) メソッド：仕様

前回、宣言だけを示したメソッド setPointOnSphere(double x) をもう少し説明しましょう。

引数の double の配列（以下 x とする）は、PolarRandomGenerator オブジェクトによって「単位球面上に一様分布」する点の座標がセットされるオブジェクトです：

double[] x = new double[3];
PolarRandomGenerator generator = new PolarRandomGenerator(new Random(1L));

generator.setPointOnSphere(x);
// x に「単位球上に一様分布」する点の座標がセットされる

ただし、x の長さを n としたとき(n = x.length)、ここでの「単位球面の次元」は n-1 になるようにします。　上記の例では、x は2次元球面上に一様分布する点の座標がセットされます。　通常使う際には、この様にしておいた方が都合がよいと思います。

setPointOnSphere(double x) メソッド：実装

3次元以上の次元で球面上に一様分布する点の座標を計算するためには、次元に関して再帰的に計算する必要がありました（→）。　そのため、球面の次元（+1）を指定して座標を計算する別のメソッド setPointOnSphere(double x, int n) を定義して、switch 文で・・・　Java コード書いた方が早いかな？

class PolarRandomGenerator{
 
    private double nextDouble(){...}
    private double nextPhi(){...}
 
    /**
     * @param x (n-1) 次元球面 (n >= 2)上において、
     * 一様分布するランダムな点の n 次元座標をセットする配列
     */
    public void setPointOnSphere(double[] x){
        setPointOnSphere(x, x.length);
    }

    private void setPointOnSphere(double[] x, int n){
        switch(n){  // n は(球面の次元 + 1) に注意。　(n-1) 次元球面上の一様分布を生成
            case 2:
                setPointOnCircumference(x);  // 円周上（1次元）上の一様分布
                break;
            
            case 3:
                setPointOnSphere2D(x);  // 普通の球面（2次元）上の一様分布
                break;
            
            default:
                setPointOnSphereND(x, n);  // n >= 4 のとき
        }
    }

残りのメソッドを以降で見ていきます。　「一般次元での単位球体内・球面上の一様分布」を計算する公式はこちら。

1次元球面（円周）上の一様分布

円周 (circumference) 上の一様分布。　引数の配列の長さは2です。

    private void setPointOnCircumference(double[] x){
        assert x.length == 2;
     
        double phi = nextPhi();
     
        x[0] = sin(phi);
        x[1] = cos(phi);
    }

簡単、簡単。

2次元球面上の一様分布

通常の球面上の一様分布。　引数の配列の長さは3です。

    private void setPointOnSphere2D(double[] x){
        assert x.length == 3;
     
        double phi = nextPhi();
        double sinTheta = nextDouble() * 2.0 - 1.0;  // 区間 [-1, 1] での一様分布
        double cosTheta = sqrt(1 - sinTheta * sinTheta);
     
        x[0] = sinTheta;
        x[1] = cosTheta * sin(phi);
        x[2] = cosTheta * cos(phi);
    }

どうってことなし。

一般次元の球面上の一様分布

(n-1) 次元球面 (n >= 4) 上の一様分布。　引数の配列の長さは n です。　内部で (n-3) 次元球面上の一様分布を生成（setPointOnSphere(x, n-2)）しています。

    private void setPointOnSphereND(double[] x, int n){
        assert n >= 4;
     
        setPointOnSphere(x, n-2);  // (n-3) 次元球面上の一様分布を生成

        double phi = nextPhi();
        double theta = nextDouble();
        double sinTheta = pow(theta, 1.0/(n-2.0));
        double cosTheta = sqrt(1 - sinTheta * sinTheta);
     
        for(int i = 0; i < n-2; i++)
            x[i] *= sinTheta;
     
        x[n-2] = cosTheta * sin(phi);
        x[n-1] = cosTheta * cos(phi);
    }
}