書本內容：見相冊
注：鑒於很多網站隨意爬取數據，可能導致內容殘缺以及引用失效等問題，影響閱讀，請認准原創網址：
https://www.cnblogs.com/lv-anchoret/category/1368696.html
preface
還記的我們上一篇說的Monte Carlo 維度詛咒嗎
上一篇算是二維的例子吧，大家看了之后是否想着寫一個一維的Monte Carlo模擬積分？（我想了，沒寫出來）
為什么要整這個嘞
光照渲染中多處涉及重積分，最終結果是要求取一個近似值，因而需要對其值進行數值估計，Monte Carlo方法就是一個較為理想的方案。
其實我們的光線追蹤器不僅用了很多向量運算，還用了很多數值分析的知識，比如之前的背景用的是最簡單的插值，柏林噪聲紋理用的是三線性插值等
ready
你需要對以下知識有了解：
二重定積分
概率密度函數
反函數
Chapter 2：One Dimensional MC Integration
積分是用來求取面積或者體積的東東，常見形式為：I=∫x2 dx
假設積分區間為0->2,我們用計算機符號表示為I = area(x2,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'r2/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) = d2,這次將不再是一個均勻一致的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) = x2/4, 0 < x < 2
P(x) = 1, x > 2
當x∈（0,2）,P(x) = ∫0->x p(r) dr = ∫0->x r/2 dr = x2/4
那么，x和r是什么關系？其實就是兩個虛擬變量，類似於程序中的參數，如果我們把x調整為有效區間的中點，那么
P(1.0) = 1/4
這就是說基於我們的pdf概率密度函數p(r) = r/2的設定下產生一個小於1的隨機數的概率為25%
這產生了一種巧妙的觀察結果，這種觀察結果是產生非均勻隨機數的許多方法的基礎。
我們想要一個函數f（），它滿足：
當我們調用f（drand48（））時，我們得到一個基於設定的Cpdf（x2/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 = x2/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(xi)/pdf(xi))
由此可知，每個積分采樣 Ii = 采樣高度（函數值）/該采樣點占的權重，所有的積分采樣取均值，實現了不同采樣點為整個積分做的貢獻是不同的
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中我們覺得函數峰高的部分對於整個積分的貢獻更多，低峰對於積分貢獻不大，所以我們對於x2積分函數采取的pdf函數為x/2
其實，pdf應該采取的就是積分函數本身，在知道答案的情況下，我們采取如下pdf函數最佳：
p(x) = (3/8)x2
則P(x) = x3/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. 對所有的積分采樣Ii = f(r)/p(r)取平均，pdf函數p()，隨機數r
這總是會收斂到正確答案，p越接近f，收斂越快
感謝您的閱讀，生活愉快~