書本內容：見相冊

注：鑒于很多網站隨意爬取數據，可能導致內容殘缺以及引用失效等問題，影響閱讀，請認準原創網址：

http://www.rzrgm.cn/lv-anchoret/category/1368696.html

 preface

還記的我們上一篇說的Monte Carlo 維度詛咒嗎

上一篇算是二維的例子吧，大家看了之后是否想著寫一個一維的Monte Carlo模擬積分？（我想了，沒寫出來）

為什么要整這個嘞

光照渲染中多處涉及重積分，最終結果是要求取一個近似值，因而需要對其值進行數值估計，Monte Carlo方法就是一個較為理想的方案。

其實我們的光線追蹤器不僅用了很多向量運算，還用了很多數值分析的知識，比如之前的背景用的是最簡單的插值，柏林噪聲紋理用的是三線性插值等

 ready

你需要對以下知識有了解：

二重定積分

概率密度函數

反函數

 Chapter 2：One Dimensional MC Integration

積分是用來求取面積或者體積的東東，常見形式為：I=∫x² dx

假設積分區間為0->2,我們用計算機符號表示為I = area(x²,0,2)

那么我們來模擬以下這個積分

我們先來看一個普通的：

那么，我們就可以寫代碼了（程序2-1）

void MC_integration(double(*f)(double x), double low, double high)
{
    size_t all = 1000000;
    double sum{ 0 };
    for (int i = 0; i < all; ++i)
    {
        double x = (high - low) * rand01() + low;
        sum += x*x;
    }
    stds cout << "I = " << (high - low) * sum / all << stds endl;
}

int main()
{
    MC_integration([](double x) {return x*x; }, 0, 2);
}

正如我們所期望的，但上述方法也可以計算那些我們用解析法無法解決的積分，如被積函數為log（sin（x）。

在圖形學中，我們經常有一些我們可以評估但不能明確寫下來的函數，或者我們只能概率評估的函數。 事實上，前面兩本書的光線追蹤代碼中的lerp函數，其實我們不知道在每個方向上看到的是什么顏色，但我們可以在任何給定的維度上統計估算它。

我 又想起來上本書剛剛寫過的區域光照渲染圖片，里面的黑色噪點真的把好精美的Cornell box 盒子給糟蹋了，那時因為我們對所有的東西做了統一的采樣，而沒有對光源進行足夠的光源采樣，所以我們可以對光源進行更多的隨機采樣，但是我們又需要控制這個比重，以防止過采樣的發生，所以，如何控制這個比例，已達到更好的效果，這就需要引入一個非常重要的概念——重要性采樣

但是在講重要性采樣之前，我們需要了解一些關于概率密度函數的東東

先看書上這張圖 

如果我們添加更多的tree，那么這張圖的長條將會增高（變長），因為他們的數量更多了，如果Height軸細分為更多的小塊，那么，每個長條的高度會降低（變短）。

離散密度函數與直方圖的不同之處在于它將頻率y軸歸一化為分數或百分比。

而連續直方圖中，我們將Height軸細分為無數個小塊，那么，縱軸將不會是一個分數，因為所有小格對應的長條基本沒有高度（每個小格代表的樹高區間基本只有一棵樹或者沒有樹，出現頻率幾乎為0）

所以對于密度函數，我們調整小格對應的影響因素，使得我們添加更多的小格時，長條不會太低

對于上述條形圖，我們采取如下：

小格的高 = 樹高介于H~H'的比例 / (H - H')

它將是很有用的！ 我們可以將其解釋為樹高的統計預測器：

即：一個介于H~H'的隨機樹高的概率 = 小格的高 * (H - H')

如果我們需要知道多個小格對應的概率，那么加和即可

其實上述公式就構成了概率密度函數，即，面積代表著對應橫坐標區間的概率，也就是上面的第二個公式

而概率密度函數就是一種連續的分數直方圖，簡稱為pdf(probability density function)

讓我們來整一個pdf，或許會有更深入的理解。如果我想要一個0~2的隨機數r它的出現概率和r自己成比例，我們將期望pdf為如下圖像：

r位于（x0,x1）的概率為C*area(p(r),x0,x1),其中C為常量,為了簡便，C一般取1

而area(p(r),0,2) = 1（概率密度函數的總面積就是總概率）

因為p(r)和r成比例，所以，設另一個常量C',則p(r) = C'r

故我們解積分：

1 = ∫_0->2 C'r dr = C'r²/2|_0->2 = 2C' =>C' = 1/2  公式2-1

所以 p(r) = r/2

我們如何用pdf p(r)生成一個隨機數？為此我們需要更多的機器。不要擔心這不會永遠持續下去！（不會像上一篇最后那個程序，永無止境。。）

給一個隨機數d = rand01(),我們知道我們寫的rand01()產生的隨機數對于0~1的每個浮點數的概率都是均勻一致的，我們需要找個函數f(d)得到我們想要的。假定e = f(d) = d²,這次將不再是一個均勻一致的pdf了，以前是的，之后我們講（把此處記為pos1遺留問題）

我們為了觀察這個函數還需要累積概率分布函數Cpdf(Cumulative probability distribution function)

記為P(x),則P(x) = area(p,-inf,x),幾何意義就是從負無窮到x的累積概率和，也就是x左邊部分的積分形成的面積

而且我們規定在沒有定義的定義域內p(x) = 0,若采用pdf函數p(r) = r/2，那么，P(x)如下：

P(x) = 0,　　x < 0

P(x) = x²/4, 0 < x < 2

P(x) = 1,　　x > 2

當x∈（0,2）,P(x) = ∫_0->x p(r) dr = ∫_0->x r/2 dr = x²/4

那么，x和r是什么關系？其實就是兩個虛擬變量，類似于程序中的參數，如果我們把x調整為有效區間的中點，那么

P(1.0) = 1/4

這就是說基于我們的pdf概率密度函數p(r) = r/2的設定下產生一個小于1的隨機數的概率為25%

這產生了一種巧妙的觀察結果，這種觀察結果是產生非均勻隨機數的許多方法的基礎。

我們想要一個函數f（），它滿足：

當我們調用f（drand48（））時，我們得到一個基于設定的Cpdf（x²/4）的一個返回值。 我們不知道那是什么，但我們知道當x為25％時應該小于1，而參數為75％應該高于1。 如果f（）為增函數，那么我們期望f（0.25）= 1.0。 

上面說了一堆，意思就是x和r無非就是函數參數，即P(1.0) = 0.25的時候，我們同樣想得到一個函數f()，它滿足f(0.25)=1.0

即：根據Cpdf得知隨機數小于1.0的概率為0.25，而存在一個f()，使得以0.25為參數的時候能夠得到那個1.0（即累積概率分布函數右端的邊界線）

這可以推廣到每個可能的輸入

即：f(P(x)) = x

也就是說f(x) = P^-1(x)

也就是P(x)的反函數。對于我們的目的而言，上述的意思就是：如果我們有了pdf函數p()和對應的Cpdf函數P(),我們為f()輸入一個隨機數，那么它就會給出我們想要的：

e = P^-1(rand01())

對于我們的pdf函數p(x) = x/2，以及對應的P(x),我們需要計算P^-1

即：如果我們有 y = x²/4

那么，x = sqrt(4y)

因此基于pdf,我們輸入隨機數，得到的就是

e = sqrt(4*rand01())

那么e∈（0,2），正如我們所期望的，我們輸入一個0.25，它就返回一個1.0

到了這里，其實上面重復了很多廢話，上面一大堆，表面上干了一個什么事情呢，用一個0~1的隨機數映射0~2這個積分區間

我們想一下我們要做的是什么，模擬x的二次方曲線的積分，我們要用隨機數模擬隨機選取積分區間的值，而我們的積分區間就是0~2，

但是沒必要這么麻煩啊，0~1隨機數隨機模擬0~2取值，直接2 * rand01()不就完啦，也就是程序2-1

不不不，它們有著天差地別，我們上述一頓操作是為了引入一個概念——Important Sample

我們通過設置pdf使得，對于積分區間的自變量取值有了權重比例，也就是說有些地方的采樣要多一點，有些地方要少一點。

也就是：I = 1/n * Σ(F(x_i)/pdf(x_i))

由此可知，每個積分采樣 I_i = 采樣高度（函數值）/該采樣點占的權重，所有的積分采樣取均值，實現了不同采樣點為整個積分做的貢獻是不同的

constexpr double pdf(const double x) {return 0.5 * x; }

void pdf1()
{
    size_t all{ 100000 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = sqrt(4 * rand01());
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

我們來解釋pos1處的遺留問題

之前的程序2-1所指的積分方案中，是均勻一致采樣，即對于任意采樣點同等對待

即,設p(x) = C

則公式2-1 寫為 

1 = ∫_0->2 C dr = Cr|_0->2 = 2C =>C = 1/2

即，開篇的方法中的pdf函數為p(x) = 1/2

所以程序2-1可以改寫為：

constexpr double pdf(const double x) {return 0.5; }

void pdf1()
{
    size_t all{ 100000 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = 2 * rand01();
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

在important sample中我們覺得函數峰高的部分對于整個積分的貢獻更多，低峰對于積分貢獻不大，所以我們對于x²積分函數采取的pdf函數為x/2

其實，pdf應該采取的就是積分函數本身，在知道答案的情況下，我們采取如下pdf函數最佳：

p(x) = (3/8)x²

則P(x) = x³/8

P^-1(x) = pow(8x,1/3)

則函數為：

constexpr double pdf(const double x) { return 3 * x * x / 8; }

void pdf1()
{
    size_t all{ 1 };
    double sum{ 0. };
    for (int i = 0; i < all; ++i)
    {
        double x = pow(8 * rand01(), 1. / 3);
        sum += x*x / pdf(x);
    }
    stds cout << "I = " << sum / all << stds endl;
}

至此，所有的內容都已經闡述完畢

我們回顧一下，因為這個是光線追蹤蒙特卡羅中的大部分概念

1. 你有一個被積函數f(x)的積分域【a,b】

2. 在域【a,b】中選擇一個非零的pdf函數p()

3. 對所有的積分采樣I_i = f(r)/p(r)取平均，pdf函數p()，隨機數r

這總是會收斂到正確答案，p越接近f，收斂越快

感謝您的閱讀，生活愉快~