Rejection Sampling

Mathematics 2023. 2. 20. 22:08

절대온도 $T$일 때 질량 $m$인 이상기체의 속력은 Maxwell-Boltzmann 분포를 따른다.

$$P(v)dv = 4\pi \left(\frac{m}{2\pi kT}\right)^{3/2}v^2 e^{-mv^2/2kT}dv$$

속력 $v$는 0부터 $\infty$까지의 값을 갖는다. 그럼 이 속력분포를 가지는 $N$개의 속력 샘플을 만들고 싶은 경우 어떻게 할까? 우선 무한대 범위까지 모두 처리할 수 없으므로 적당한 속력 범위 내의 샘플만 만들자. 계산하기 쉽도록 하기 위해서 $x^2 = mv^2 /kT$인 차원없는 변수 $x$를 사용한다. 확률이 보존되어야 하므로 $P(v)dv= f(x)dx$에서 $x$에 대한 확률밀도함수는

$$ f(x) = \sqrt{\frac{2}{\pi}} x^2 e^{-x^2/2}$$

로 주어진다. 가능한 샘플 $x$는 다음의 단계로 얻어진다.

1. 주어진 $x$ 범위 $[x_\text{min}, x_\text{max}]$ 사이에서 랜덤하게 한 값 $x$를 취하고,

2. $[0,1]$ 사이에서 랜덤한 값 $p$을 취한 후 

3. $ p> f(x)$이면 $x$는 버리고 다시 과정 1로 돌아간다. 

 

다음 예는 $x$를 $[0,4]$사이에서 500 구간으로 나누어서 만든 50000개의 샘플의 분포와 Maxwell-Boltzmann 분포를 비교한 것이다.

그림의 붉은점은 accepted된 값이고 푸른점은 rejected된 값이다.

    const CSize sz = CSize(500, 500);
    raster.SetDimensions(sz, 24);
    srand(unsigned(time(0)));
    int hist[500] = {0};
    int required = 50000;
    const double xmin  = 0, xmax = 4;
    while (required) {
        double x = double(rand()) / RAND_MAX; //[0,1]
        double v = (xmax - xmin) * x + xmin;  //[xmin,xmax]
        double p = double(rand()) / RAND_MAX; //[0,1]
        double maxwell = MaxwellDist(v);
        CPoint q = CPoint(int(x * (sz.cx - 1)), int(p * (sz.cy - 1)));
        if (p <= maxwell) {
            raster.SetPixel(q.x, sz.cy - 1 - q.y, 0xFF0000);
            hist[q.x]++;
            required--;
        }
        else
            raster.SetPixel(q.x, sz.cy - 1 - q.y, 0xFF);	
    }
    int s = 0;
    for (int i = 0; i < sz.cx; i++) 
        s += hist[i];
    FILE *fp = fopen("maxwell.txt", "w");
    for (int i = 0; i < sz.cx; i++)
        fprintf(fp, "%f, %f\n", (xmax - xmin) * double(i) /(sz.cx - 1) + xmin, \
                                double(hist[i]) / s);
    fclose(fp);
728x90
Posted by helloktk
,

p=0.2, imax=255, imin=0

int AddImpulseNoise(CRaster& raster, double p, int imax, int imin, CRaster& noised) {
    srand(unsigned(time(NULL)));
    CSize sz = raster.GetSize();
    p = max(min(p, 1),0); //noise probability;
    int np = int(p * sz.cx * sz.cy);
    // clone;
    noised = raster;
    imax = max(min(255, imax), 0); //maximum impulse value;
    imin = min(max(0, imin), 255); //minimum impulse value;
    for (int count = 0, turn = 0; count < np; count++) {
        //x in [0, sz.cx - 1];
        int x = int((sz.cx - 1) * double(rand()) / RAND_MAX); 
        //y in [0, sz.cy - 1];
        int y = int((sz.cy - 1) * double(rand()) / RAND_MAX); 
        if (turn == 0) {
            turn = 1;
            noised.SetPixel(x, y, imax);
        } else {
            turn = 0;
            noised.SetPixel(x, y, imin);
        }
    }
    return 1;
}
728x90

'Image Recognition > Fundamental' 카테고리의 다른 글

Ellipse Fitting  (0) 2024.03.02
Bilateral Filter  (0) 2024.02.18
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Posted by helloktk
,

면적이 $A$인 단순 평면 도형이 있을 때 그 둘레의 길이 $L$과의 사이에는 다음 부등식이 성립한다:

$$ L^2 \ge 4\pi  A$$

정 $n$ 각형에 대해서 이 부등식을 체크해보면

$$ \frac{L^2}{4\pi A} = \frac{\tan(\pi/n)}{\pi/n}=\left\{\begin{array}{cc} 1.654 & \text{정삼각형}\\ 1.273&\text{정사각형} \\ 1.156&\text{정오각형}\\ 1.034&\text{정십각형}\end{array} \right.$$

증명: 평면에서 둘레의 길이가 $L$인 한 폐곡선을 곡선을 따라 움직인 길이를 매개변수로 사용해서 표현하자: $(x(s), y(s)),~s\in [0,L]$. 이 폐곡선을 복소평면으로 옮기면 $z(s)= x(s) + i y(s)$로 표시할 수 있다. $z(s)$는 $s$에 대해서 주기가 $L$ 복소함수로 볼 수 있다. Fourier series를 사용하기 편리하게 하기 위해서 매개변수를 바꿔 $z(t) = x(Lt/2\pi) + i y(Lt/2\pi)$로 표현하자. 그러면 $t$에 대한 주기가 $2\pi$가 되고 다음과 같이 Fourier 전개를 쓸 수 있다. (영상처리에서는 물체의 윤곽선에 대한 Fourier  계수를 이용해서 물체의 모양을 분류하는데 이용이 된다)

$$ z(t) = \sum_{k=-\infty}^\infty a_k e^{i k t}$$

먼저 $L^2$항을 뽑아내기 위해 다음의 적분을 고려하자. $s=Lt/2\pi$로 놓으면 $z(s)$은 곡선의 길이를 매개변수로 사용한 것이어서  $|dz/ds|=1$이므로

$$  \int_0^{2 \pi} \left|  \frac{dz}{dt} \right|^2 dt = \int_0^{2\pi} \left|\frac{dz}{ds} \right| ^2 \left( \frac{ds}{dt}\right)^2 dt = \frac{L^2}{4\pi^2} \int_0^{2 \pi} dt = \frac{L^2}{2\pi}$$ 

좌변에 Fourier 전개를 대입하면 (note: $ \int_0^{2\pi} e^{i (k-j)t} dt = 2\pi \delta _{kj}$)

\begin{align} \int_0^{2\pi}\left|  \frac{dz}{dt} \right|^2 dt= 2 \pi \sum_{k=-\infty}^\infty k^2 |a_k|^2 \end{align}

을 얻는다. 따라서

$$ \frac {L^2 }{4\pi^2 } = \sum_{k=-\infty}^\infty k^2 |a_k|^2 $$

다음으로 도형의 면적으로 표현하자.

\begin{align} A &= \frac {1}{2} \int_0^ {2\pi } \left(x \frac{dy}{dt}- y \frac{dx}{dt}  \right)dt \\ &= \frac{1}{2} \text {Re} \int_0^{2\pi} i z \overline{\frac{dz}{dt}}  dt= \pi \sum _{k=-\infty}^\infty k |a_k|^2 \end{align}

그런데 $$ \sum_k k^2 |a_k| ^2 \ge   \sum _k k |a_k|^2$$

이므로

$$ \frac{L^2}{4\pi} \ge A$$

을 얻는다. 그리고 등호 성립하려면 Fourier 계수가 $k=0, 1$이외에는 모두 0이어야 된다. 즉

$$ z(t) = a_0 + a_1 e^{i t}$$

인 경우인데, 이는 복소평면에서  중심이 $a_0$, 반지름이 $|a_1|$인 원을 나타낸다. 변분법을 이용해서도 원이 주어진 둘레길이를 가지는 평면도형 중에서 가장 넓은 면적의 도형임을 쉽게 보일 수 있다. 

https://kipl.tistory.com/187

 

주어진 길이의 폐곡선으로 가둘 수 있는 최대 면적의 도형은?

평면에서 주어진 길이의 폐곡선으로 가둘 수 있는 최대의 면적은 얼마이고 그 모양은 어떤 것일까? 이 문제는 변분법을 이용하면 쉽게 해결할 수 있다. 평면에서 폐곡선은 매개변수 $t$의 함수로

kipl.tistory.com

 

728x90
Posted by helloktk
,