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' 카테고리의 다른 글

Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
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

 

주어진 둘레 길이를 갖는 최대 면적 폐곡선

 

kipl.tistory.com

 

728x90

'Mathematics' 카테고리의 다른 글

Rejection Sampling  (0) 2023.02.20
Brachistochrone inside the Earth  (0) 2023.01.25
Snell's law and Brachistochrone Curve  (0) 2023.01.25
Tautochrone Curve(등시곡선)  (0) 2023.01.16
Fourier transform of the Heviside step function  (0) 2023.01.12
Posted by helloktk

댓글을 달아 주세요

같은 질량의 두 카트가 3개의 동일한 용수철로 그림과 같이 연결되어 있다. 두 카트를 움직이게 하면 다양한 양상의 진동을 보이는데 그중에 가장이 기본이 되는 두 가지 진동형태(normal mode)가 있다. 하나는 두 물체가 동일한 진동을 하는 것이고(각 카트의 평형위치에서 변위로 표현하면 $x_A(t) = x_B(t) = x_0 \cos(\omega_1 t)$) 다른 하나는 서로 반대방향 진동을($x_A(t) = -x_B(t) = x_0 \cos (\omega_2 t)$) 하는 것이다.  두 진동의 각진동수 비 $\omega_2 / \omega_1$=?

  1. $\sqrt{3}$
  2. $2$
  3. $2\sqrt{2}$
  4. $3$
  5. $5$

동영상: https://youtu.be/cu4TvUwk17g

728x90

'Physics > 역학' 카테고리의 다른 글

두 물체는 얼마나 가까워질 수 있는가?  (0) 2023.01.31
진동의 주기는?  (0) 2023.01.31
얼마나 빨리 회전해야 되는가?  (0) 2023.01.30
반구의 진동 주기는?  (0) 2023.01.28
반구의 회전관성은?  (0) 2023.01.27
Posted by helloktk

댓글을 달아 주세요

같은 질량의 두 물체가 일정한 거리 $d$ 만큼 떨어진 상태에서 서로 반대 방향을 향해 $v_0\ne 0$로 출발한다. 두 물체 사이에는 중력만 작용한다.

1. 언젠가 두 물체는 충돌할까?

2. 가장 가까워졌을 때 두 물체의 사이거리는?

 

728x90

'Physics > 역학' 카테고리의 다른 글

각진동수의 비는?  (0) 2023.01.31
진동의 주기는?  (0) 2023.01.31
얼마나 빨리 회전해야 되는가?  (0) 2023.01.30
반구의 진동 주기는?  (0) 2023.01.28
반구의 회전관성은?  (0) 2023.01.27
Posted by helloktk

댓글을 달아 주세요