Geometry & Recognition

이전 포스팅(https://kipl.tistory.com/731)에서 단위구에서 두 점 사이거리의 평균을 구하는 과정에서 사이거리에 대한 확률밀도함수를 소개하였고, 또 단위원에서 사이거리에 대한 확률밀도함수를 구체적으로 구했다. 이제 단위원에서의 기법을 이용해서 단위구에서 사이거리 분포에 대한 확률밀도함수를 구하자. 두 점의 사이거리 $s$의 확률밀도함수는

$$P(s)=\frac{1}{B_1^3}{\iint }_{\text{unit ball}}{d}^{3}r{d}^3{r}^{\prime }\delta (|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|-s),B_1=\frac{4\pi }{3}$$

로 쓰인다. 주어진 단위구 내부의 한 지점 $\overrightarrow{r}$에 대해서 적분은 $$S(\overrightarrow{r},s)={\int }_{\text{unit ball}}{d}^3{r}^{\prime }\delta (|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|-s)$$

는 $\overrightarrow{r}$에서 반지름 $s$인 구면이 단위구에 포함된 면적을 나타낸다. 그리고 이 값은 $\overrightarrow{r}$의 방향에 무관하게 단위구 중심에서 거리에만 의존함을 쉽게 알 수 있다, 즉 

$$S(\overrightarrow{r},s)=S(r,s)$$

이다. 따라서 $r+s$가 1보다 작을 때와 클 때 두 경우를 별도로 고려해야 한다. 작은 경우는 반지름 $s$인 구면이 완전히 단위구 내부에 포함되므로 

$$r+s<1S(r,s)=4\pi s^2$$

이고, 큰 경우는 (단위원의 경우 그림을 참조하면: https://kipl.tistory.com/732) $4\pi s^2$에서 단위구 밖으로 나가는 구면캡(spherical cap)의 면적을 제외하면 된다. 구면캡의 경도각 범위가 $\cos \theta =\frac{1-r^2-s^2}{2rs}$이므로 구면캡의 입체각은 

$$\mathrm{\Delta }\mathrm{\Omega }=2\pi {\int }_0^{\theta }\sin \theta d\theta =2\pi (1-\frac{1-r^2-s^2}{2rs})=2\pi \frac{(r+s{)}^2-1}{2rs}$$

$$\to r+s>1S(r,s)=(4\pi -\mathrm{\Delta }\mathrm{\Omega })s^2=4\pi s^2\frac{1-(r-s{)}^2}{4rs}$$

따라서 사이거리에 대한 확률밀도함수는

$$P(s)=\frac{1}{B_1^2}4\pi {\int }_0^1{r}^{2}drS(r,s)$$

$$=\frac{1}{B_1^2}4\pi {\int }_0^{1-s}{r}^{2}dr(4\pi s^2)+\frac{1}{B_1^2}4\pi {\int }_{1-s}^1{r}^{2}dr\left(4\pi s^2\frac{1-(r-s{)}^2}{4rs}\right)$$

$$=\frac{3}{16}{s}^2(2-s{)}^2(4+s)$$

거리의 평균:

$$<s>={\int }_0^{2}P(s)sds=\frac{36}{35}$$

거리제곱의 평균:

$$<s^2>={\int }_0^{2}P(s)s^2=\frac{6}{5}$$

거리역수의 평균:

$$<\frac{1}{s}>={\int }_0^{2}P(s)\frac{1}{s}ds=\frac{6}{5}$$

단위원에서 선택한 임의의 두 지점의 사이거리 $s$가 $0\le s\le 2$이 확률밀도함수는 

$$P(s)=\frac{1}{D_1^2}{\iint }_{\text{unit disk}}{d}^{2}r{d}^2{r}^{\prime }\delta (|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|-s),D_1=\pi$$

먼저

$$L(\overrightarrow{r},s)={\int }_{\text{unit disk}}{d}^2{r}^{\prime }\delta (|{\overrightarrow{r}}^{\prime }-\overrightarrow{r}|-s)$$

은 델타함수 제한조건 때문에 $\overrightarrow{r}$을 중심으로 한 반지름 $s$인 원과 단위원의 겹치는 영역의 경계의 길이를 의미함을 알 수 있다. 이 $L(\overrightarrow{r},s)$가 구해지면 $P(s)$는 이 값을 단위원에 대해 적분을 하면 되는데, 그 값은 단위원의 중심엣 $\overrightarrow{r}$가 거리에만 의존함을 쉽게 알 수 있다. 따라서 $\overrightarrow{r}$가, 예를 들면  $x$축 위에 있는 경우만 고려하면 충분하다. 또, $L(\overrightarrow{r},s)=L(|\overrightarrow{r}|=r,s)$이므로

$${\int }_{\text{unit disk}}{d}^{2}rL(r,s)=2\pi {\int }_0^{1}rdrL(r,s)$$

이제 $\overrightarrow{x}$을 중심으로 한 반지름 $s$인 원이 완전히 단위원에 포함이 되는 경우는 원과 단위원의 겹치는 영역이 반지름 $s$인 원이므로 

$$r+s<1\to L(r,s)=2\pi s$$

이고, 원호 일부가 단위원 밖으로 나가는 경우는 두 가지 경우(그림의 $\theta$가 예각 또는 둔각)를 생각할 수 있는데,

두 경우 모두 단위원과 겹치는 영역의 원호의 길이는 그림에서 $2(\pi -\theta )s$인데, $$\sin \phi =s\sin \theta ,\cos \phi -s\cos \theta =r$$

$$\to \cos \theta =\frac{1-r^2-s^2}{2rs}$$

$$r+s>1\to L(r,s)=2(\pi -\theta )s=2s(\pi -{\cos }^{-1}\frac{1-r^2-s^2}{2rs})$$

로 표현된다. 따라서 사이거리에 대한 확률밀돟함수는

$$P(s)=\frac{2\pi }{{\pi }^2}{\int }_0^{1}rdrL(r,s)$$$$=4s{\int }_0^{1-s}rdr+\frac{4s}{\pi }{\int }_{1-s}^{1}rdr(\pi -{\cos }^{-1}\frac{1-r^2-s^2}{2rs})$$

$$=\frac{4s}{\pi }{\cos }^{-1}\frac{s}{2}-\frac{2s^2}{\pi }\sqrt{1-\frac{s^2}{4}}$$

이 분포를 이용하면 단위원 내부에서 선택된 두 점간의 평균거리는

$$⟨|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|⟩={\int }_0^{2}P(s)sds=\frac{128}{45\pi }=0.9054$$

$$⟨|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }{|}^2⟩={\int }_0^{2}P(s)s^{2}ds=1$$

$$⟨\frac{1}{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}⟩={\int }_0^{2}P(s)\frac{1}{s}ds=\frac{16}{3\pi }$$

단위구 내의 두 지점을 ${\overrightarrow{r}}_1$, ${\overrightarrow{r}}_2$라면 두 지점이 선택될 확률이 

$$P({\overrightarrow{r}}_1,{\overrightarrow{r}}_2)=\frac{d{r}_1^3}{B_1}\times \frac{d{r}_2^3}{B_1},B_1=\text{volume of unit ball}=\frac{4\pi }{3}$$이므로 두 지점의 사이거리의 평균은

$$⟨|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2|⟩=\frac{1}{B_1^2}{\iint }_{\text{unit ball}}|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2|d{r}_1^{3}d{r}_2^3$$

${\overrightarrow{r}}_1$과 $\overrightarrow{r_2}$의 사이각이 $\theta$일 때 고정된 ${\overrightarrow{r}}_1$에 대해서 ${\overrightarrow{r}}_2$의 극좌표 고도각으로 선택하면 $r_2$에 대한 적분은

$$I_2={\int }_{\text{unit ball}}|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2|d{r}_2^3$$

$$=2\pi {\int }_0^1{r}_2^{2}d{r}_2{\int }_{-1}^{1}d\cos \theta \sqrt{r_1^2+r_2^2-2r_1{r}_2\cos \theta }$$

$$=\frac{2\pi }{3}{\int }_0^1{r}_2^{2}d{r}_2\frac{(r_1+r_2{)}^3-|r_1-r_2{|}^3}{r_1{r}_2}$$

$$=\pi (1+\frac{2}{3}{r}_1^2-\frac{r_1^4}{15})$$

따라서 사이거리의 평균은

$$⟨|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2|⟩=\frac{4\pi }{B_1^2}{\int }_0^1{r}_1^{2}d{r}_1\pi (1+\frac{2}{3}{r}_1^2-\frac{r_1^4}{15})$$

$$=\frac{1}{B_1^2}\frac{64{\pi }^2}{35}=\frac{36}{35}$$

이전 포스팅에서 거리 역수의 평균이 $\frac{6}{5}$임을 보였는데, 평균거리가 이 값의 역수로 주어지지는 않는 점이 의외다.

아래 그림은 거리에 대한 확률밀도함수를 Monte Carlo 시뮬레이션으로 구한 결과이다. 실제로 거리에 따른 확률밀도함수를 임의의 차원에서 구할 수 있는데, 3차원의 단위구는(ref: J. Math. Phys., Vol. 41, No. 4, April 2000, https://kipl.tistory.com/733)

$$P_3(s)=\frac{3}{16}{s}^2(2-s{)}^2(4+s),0\le s\le 2$$

로 계산이 된다. 이 확률밀도함수의 최대값은 $s=1$이 아니라 $s=\frac{1}{5}\sqrt{105}-1=1.04939$ 에서 나온다. $P_3(s)$를 이용하면$$<s>={\int }_0^2{P}_3(s)sds=\frac{36}{35}$$를 얻고, 거리의 역수의 평균은$$<\frac{1}{s}>={\int }_0^2{P}_3(s)\frac{1}{s}ds=\frac{6}{5}$$도 얻을 수 있다. 

#define TWO_PI 6.28318530717958647692
// picking a random vector inside an unit ball;
void randomVector(double &x, double &y, double &z) {
    double a = double(rand()) / RAND_MAX;		//[0,1] = radius;
    double b = double(rand()) / RAND_MAX * 2 - 1; //[-1,1] = cos(theta);
    double c = double(rand()) / RAND_MAX * TWO_PI; //[0, 2PI] = azimuthal; 
    a = pow(a, 1.0 / 3.0);
    x = a * sqrt(1 - b * b) * cos(c);
    y = a * sqrt(1 - b * b) * sin(c);
    z = a * b;
}
// Calculate the distance between two randowm points selected 
// inside an unit ball;
double randomDistance() {
    double x1, y1, z1, x2, y2, z2;
    randomVector(x1, y1, z1);
    randomVector(x2, y2, z2);
    x1 -= x2;
    y1 -= y2;
    z1 -= z2;
    return sqrt(x1 * x1 + y1 * y1 + z1 * z1);
}
void meanDistSphere(CRaster& raster) {
    double R = 100;
    int trials = 5000000;
    srand(unsigned(time(NULL)));
    std::vector<int> hist(2 * int(R + .5) + 1, 0);
    for (int i = trials; i-->0;) {
        double d = R * randomDistance();
        hist[int(d)]++;
    }
    double s = 0, sx = 0, sxx = 0;
    double maxh = 0;
    for (int i = hist.size(); i--> 0; ) {
        double h = hist[i];
        maxh = max(h, maxh);
        s += h;
        sx += h * i;
        sxx += h * i * i;
    }
    double mean = sx / s;
    double var = (sxx - sx / s) / s;
    double sdev = sqrt(var);
    // CRaster raster;
    raster.SetDimensions(hist.size(), hist.size(), 24);
    CSize sz = raster.GetSize();
    memset(raster.GetDataPtr(), 0xFF, raster.GetDataSize());
    for (int x = sz.cx; x--> 0; ) {
        int h = int(double(sz.cy) * hist[x] / maxh);
        for (int y = h; y--> 0; ) 
            raster.SetPixel(x, sz.cy-1-y, 0);
    }
    // SaveRaster(raster, "hist.bmp");
}