y = sin(x)^2함수 하나가 있다고 해봅시다. 그리고 당신은 0과 pi 사이의 커브 아래쪽에 있는 면적을 알고 싶습니다.
우리는 분석적으로 정확한 답을 얻기 위해서 대수와 미적분을 사용하여 이 특정 문제를 해결할 수 있을 것입니다(답은 pi/2), 우리가 그럴 수 없거나 이런 방법으로 해결하고 싶지 않다고 해봅시다.
이 문제를 해결하는 다른 방법은 몬테 카를로 적분을 사용하는 것입니다, 이것은 당신이 문제를 수치적으로 해결하고 근사된 결과를 얻을 수 있게 해줍니다.
어떻게 하면 이렇게 할 수 있을까요:
0에서 pi 사이의 난수를 선택합니다.
선택한 값을 y = sin(x)^2 함수의 x에 넣어주고 y 값을 얻습니다.
이 함수의 y 값의 평균을 얻기 위해서 이 것을 여러 번 수행하고 평균을 얻습니다.
함수가 직사각형(rectangle)이라고 해봅시다, y의 평균은 직사각형의 높이로 사용할 수 있습니다, 그리고 우리가 0과 pi 사이를 알아보는 중이기 때문에 pi를 폭으로 사용할 수 있습니다.
직사각형의 면적을 얻기 위해서 폭과 높이를 곱합니다, 이것은 곡선 아래쪽의 추정된 면적입니다.
이것이 당신이 해야할 모든 것입니다.
몬테카를로 적분은 단순하다는 점과 매우 높은 차원에서도 잘 작동하는 점에서 꽤 강력합니다.
당신이 상상하듯, y값의 평균을 얻기 위해서 더 많은 샘플을 얻을 수록, 더 나은 추정치가 됩니다. 불행히도, 에러를 절반으로 줄이기 위해서 4배의 샘플 수를 가져야 합니다, 그래서 높은 수준의 정확도가 필요하다면 올바른 답(수렴)을 얻기 위해서 시간이 걸릴 수 있습니다.
두 가지 모두 샘플링 공간에 대해서 더 고른 적용 범위를 제공합니다. 이것은 당신의 샘플로부터 누락된 정보간에 큰 틈은 없을 것이라는 것입니다.
또 다른 방법은 stratified sampling 입니다, 여기서는 샘플링 공간을 몇개의 섹션 수로 나눕니다, 그리고 각 섹션 범위에서 난수를 선택합니다, 각 섹션안에서 확실히 샘플을 가질 수 있도록 합니다. 이것은 무작위성을 유지하면서 샘플링 공간에 대해서 더 고른 적용 범위를 제공합니다.
당신은 "만약 내가 100개의 샘플이 있으면, 그 공간에서 1/100 마다 고르게 샘플링 할 것입니다."라고 말하고 싶은 유혹을 받을 수 있습니다. 균일/정규 샘플링은 앨리어싱을 포함한 몇가지 문제 뿐만아니라 난수가 줄 수 있는 좋은 수학적인 특성을 잃기도합니다. (유리수가 아닌 위치에서 샘플링 할 수 있는 것과 같은)
Stratified sampling의 변경은 Pixar에 의해 발명된 "jittered grid"라 불리는 기술입니다. 이 기술은 샘플링을 고르게 하되 각 샘플에 작은 랜덤 값을 추가합니다.
정말 많은 다른 기술들이 있지만 이 글을 길어지게 할 것 같아서, 여기까지만 할 것입니다.
More General Monte Carlo Integration
마지막 섹션은 몬테카를로 적분의 단순화된 버젼이었습니다. 이것은 균일한 난수를 사용했기 때문에 가능했습니다.
몬테카를로 적분은 균일한 난수가 아닌 임의의 분포를 가진 난수들을 사용할 수 있습니다.
이 과정은 대부분 동일하지만 몇가지 차이점이 있습니다.
이전 섹션에서, 그것이 직사각형인 것 처럼, 평균 높이를 얻었고 폭을 곱하여 커브의 아래쪽 면적의 추정치를 얻었습니다.
첫번째 변경점은 폭이 곱해지는 부분이 반복문 내로 들어가는 것입니다. 평균 높이를 계산하는 대신에, 우리는 평균 직사각형의 면적을 계산합니다.
수학적으로 당신은 동일한 답을 얻습니다, 그래서 이상한 점은 없습니다.
두번째 변경점은 폭을 곱하는 대신에, 방정식에 연결시킨 그 수가 선택될 확률로 나누어줍니다.
0과 pi 사이의 샘플을 얻는 우리의 함수의 경우, 그 범위에서 어떤 단일 수의 선택 확률은 1/pi 와 같습니다. 우리가 그것으로 나누는 것은 그냥 pi 를 곱하는 것과 같은 의미입니다, 그래서 이것은 수학적으로 우리가 이전에 했던 것과 동일합니다.
여기에 더 일반화된 몬테카를로 적분이 있습니다.
원하는 특정 난수 분포를 사용하여 0과 pi 사이의 난수를 선택합니다.
선택된 값을 y = sin(x)^2 함수에 x에 넣고 y를 얻습니다.
함수의 추정된 면적을 얻기위해서 y 값을 그 수가 선택될 확률(PDF(x)라고도 함)로 나눕니다.
이것을 여러번하고 평균 내어 결과를 얻습니다.
여기에 여전히 균일하게 분포된 난수를 사용하는 일반적인 몬테카를로 적분 코드가 있습니다.
double GeneralMonteCarlo()
{
size_t numSamples = 10000;
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_real_distribution<double> dist(0.0f, 1.0f);
auto InverseCDF = [](double x) -> double
{
return x * c_pi;
};
auto PDF = [](double x) -> double
{
return 1.0f / c_pi;
};
double estimateSum = 0.0;
for (size_t i = 1; i <= numSamples; ++i)
{
double rnd = dist(mt);
double x = InverseCDF(rnd);
double y = sin(x)*sin(x);
double pdf = PDF(x);
double estimate = y / pdf;
estimateSum += estimate;
}
double estimateAverage = estimateSum / double(numSamples);
return estimateAverage;
}
흥미롭게도, PDF로 나누는 것은 마지막 섹션에서 폭을 곱해주는 것과 수학적으로 동일합니다 - 말 그대로 pi(폭)로 곱해주는 것이 됩니다. 유일한 차이점은 곱셈을 끝에 하지 않고 반복문 내로 가져온 것입니다.
최적화를 위해서, 당신은 확실히 나누기를 다시 바깥으로 옮길 수 있습니다(그리고 그것은 다시 곱셈으로 변함), 그러나 나는 가능하면 코드가 핵심 원리에 가까운 형태로 남겨두고 싶습니다.
Non Uniform Random Number Distributions
다른 난수 분포로 부터 샘플링 해봅시다. y = sin(x) 식에서 난수를 생성해봅시다. 당신은 아래 그래프에서 이것을 우리가 적분하는 y = sin(x)^2과 비교해볼 수 있습니다. 두 식은 꽤 비슷한 모양입니다!
몬테카를로 적분에 y = sin(x)를 난수 분포로 사용하기 위해서, 우리는 정규화된 PDF와 InverseCDF를 계산할 필요가 있습니다.
PDF로부터 숫자를 생성하기 위해서, 0에서 1 사이의 난수를 얻습니다. 그리고 inverse CDF 함수에 연결합니다: InverseCDF(x) = 2*asin(sqrt(x))
여기에 이 PDF와 inverseCDF를 사용하여 몬테카를로 적분을 수행하는 코드 조각이 있습니다.
double ImportanceSampledMonteCarlo()
{
size_t numSamples = 10000;
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_real_distribution<double> dist(0.0, 1.0);
auto InverseCDF = [](double x) -> double
{
return 2.0 * asin(sqrt(x));
};
auto PDF = [](double x) -> double
{
return sin(x) / 2.0f;
};
double estimateSum = 0.0;
for (size_t i = 1; i <= numSamples; ++i)
{
double rng = dist(mt);
double x = InverseCDF(rng);
double y = sin(x)*sin(x);
double pdf = PDF(x);
double estimate = y / pdf;
estimateSum += estimate;
}
double estimateAverage = estimateSum / double(numSamples);
return estimateAverage;
}
이것과 균일 랜덤 샘플링을 비교하기 위해서, 먼저 50,000,000 개 이상의 샘플을 만드는 과정을 보여줄 것입니다, 그리고나서 y = sin(x) 모양의 PDF를 사용할 것입니다.
Uniform 즉, 1/pi:
sin(x):
이전에 우리가 이야기 한 것 처럼, 샘플이 4배 늘어날 때마다, 표준편차(분산의 제곱근)는 절반으로 줄어듭니다. 이것이 패스트레이싱이 오래걸리는 이유입니다. 만약 당신이 패스트레이싱을 모른다면, 이것이 현대 애니매이션 영화가 렌더링 하는데 오랜시간이 걸리는 이유입니다.
그 결과, 당신은 적분하려는 함수의 유사한 모양의 PDF를 사용하면 추정치의 분산이 많이 낮은 것을 볼 수 있습니다. 우리는 더 좋고, 신뢰할 수 있는 답을 더 적은 샘플로 부터 얻을 수 있습니다.
Bad Random Number Distributions
만약 당신이 적분하려는 함수와 아주 다른 모양의 PDF를 사용하면, 당신은 더 큰 분산을 얻을 것이고 수렴하는데 더 오래걸릴 것입니다. 이것은 완전 실망스러운 일입니다.
y = (x/pi)^5 을 해봅시다, 이것은 우리가 적분하려는 함수와 비슷하게 보이지 않습니다.
여기에 PDF와 inverseCDF 가 있습니다:
PDF(x) = (6/pi) * (x/pi)^5
inverseCDF(x) = (x * pi^6)^(1/6)
여기에 50,000,000 샘플이 있습니다:
그리고 비교를 위해서 여기 다시 균일 샘플이 있습니다:
보시다시피, 정답에 근접하고 있습니다, 그러나 같은 결과(분산의 양)를 위해서 균일 샘플링과 비교해 10배 더 오래 걸립니다. 아얏!
Perfect Random Number Distributions
우리가 정말 운이 좋아서 우리가 적북하려는 함수와 완벽하게 일치하는 PDF와 inverseCDF를 얻었다고 해봅시다. 무슨 일이 일어날까요?
y = sin(x) 형태를 가진 난수 분포를 사용하여 y = sin(x) 함수를 적분하는 것을 확인해봅시다.
우리는 이미 PDF와 inverseCDF을 이전에 계산했습니다:
PDF(x) = sin(x)/2
inverseCDF(x) = 2 * asin(sqrt(x))
여기에 50,000,000 샘플로 그것을 수행합니다:
와! 보는 것과 같이, 첫번째 샘플부터 올바른 답을 얻었습니다, 분산(불균일성)이 0이고요. 그리고 50,000,000 샘플에 대해서 답을 계속해서 유지했습니다.
이것은 꽤 깔끔한 개념입니다, 그리고 당신이 "cosine weighted hemisphere sampling"에 대해서 안다면, 정확히 이런 일을 합니다.
Cosine weighted hemisphere samples는 cos(theta)를 라이팅 계산에서 제거될 수 있게 가중되어집니다, 왜냐하면 난수 분포가 당신을 위해 그것을 처리하기 때문입니다.
그것은 기본적으로 방정식에서 난수성을 제거합니다.
불행하게도 패스 트레이싱에서는 해당 텀 보다 더 많은 변수와 난수성이 있지만 도음이 됩니다.
이 외에도, multiple importance sampling을 포함한 다른 분산 감소 기법에도 관심이 있으면 살펴보시기 바랍니다.
Closing
이 블로그의 글에 들어가면서 나는 "에이 별것 없어, 나는 간단한 함수 몇 개 만들고, 그들의 PDF와 inverse CDF를 계산하고 내 갈길 가야지"라고 생각했습니다.
나는 내가 시도한 거의 모든 간단한 함수들이 그 과정을 거치는 것이 불가능하다는 점을 믿을 수 없었습니다.
예를들어, 당신은 x = sin(y)를 얻을 수 있고 y = asin(x)를 얻기 위해서 y에 대해 풀 수 있습니다, 그러나 만약 당신이 x = sin(y)*y를 y에 대해 풀고자 하면, 그날은 나쁜 날이 될 것입니다.
내 생각에 미래에 이런 것이 필요하면, 나는 (x, y) 데이터 점을 (y, x) 데이터 점으로 fitting a curve를 시도할 것입니다, 그러나 이런 일을 하는 다른 방법들도 많이 있습니다.
그나저나 적분하는 동안 만약 표준편차(즉, 분산의 제곱근)를 어떻게 계산했는지 궁금하다면, 분산은 "평균으로 부터 구한 차이의 제곱의 평균" 입니다. 즉, 만약 적분하려고 하는 항목의 정답을 알고있다면, 다음과 같이 표준편차를 계산할 수 있습니다:
// Variance is "The average of the squared differences from the mean"
double difference = integration - actualAnswer;
double differenceSquared = difference * difference;
averageDifferenceSquared = Lerp(averageDifferenceSquared, differenceSquared, 1.0 / double(i));
double stdDev = sqrt(averageDifferenceSquared);
적분은 현재의 평균 추정치입니다. (만약 100 개의 샘플이 있으면, 그것은 100개의 샘플의 평균입니다)
Anders Lindqvist (@anders_breakin)는 만약 당신이 이 글을 좋아하였다면 좋아할만한 몬테카를로, importance sampling 그리고 multiple importance sampling으 설명하는 블로그 글을 쓰고 있습니다. 그를 팔로우하면, 그 글이 곧 나올 것입니다.
