)
做光學(xué)仿真和隨機模擬這些年我發(fā)現(xiàn)自己繞不開一個坎所有編程語言和仿真軟件能直接生成的隨機數(shù)幾乎都是均勻分布。可現(xiàn)實世界里真正常用的是高斯分布——也就是正態(tài)分布。測量噪聲是高斯分布光斑的能量分布接近高斯分布人群身高的統(tǒng)計也是高斯分布。于是“均勻分布產(chǎn)生高斯分布”就成了一個高頻問題網(wǎng)上搜一下相關(guān)討論特別多連LightTools這種光學(xué)仿真軟件里怎么設(shè)置高斯分布都被反復(fù)問。這篇文章我打算把這件事徹底講透從數(shù)學(xué)原理講到代碼實現(xiàn)再落到LightTools這類工程工具里的實際操作把我踩過的坑和驗證過的方法都整理出來。1. 均勻分布和高斯分布先搞清楚我們要干什么1.1 兩個分布到底差在哪里均勻分布的概率密度函數(shù)是一條水平直線在定義區(qū)間內(nèi)每個點出現(xiàn)的概率一樣。打個比方均勻分布就像抽簽箱子里十個球抽中任何一個的概率都是十分之一。而高斯分布是一條鐘形曲線中間高、兩邊低絕大多數(shù)樣本落在均值附近極端值幾乎不會出現(xiàn)。它的概率密度函數(shù)長這樣[ f(x) \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{(x-\mu)^2}{2\sigma^2}} ]這里的μ是均值σ是標準差σ2是方差。μ決定了鐘形曲線在x軸上的位置σ決定了曲線是“胖”還是“瘦”。σ越小曲線越尖數(shù)據(jù)越集中σ越大曲線越平數(shù)據(jù)越分散。高中數(shù)學(xué)里大家可能背過這個公式但沒有多少人認真想過它背后的幾何含義。高斯分布之所以無處不在本質(zhì)上是因為自然界里大多數(shù)“誤差”和“波動”都是大量微小獨立因素疊加的結(jié)果——一個人的身高受幾百個基因位點影響光學(xué)系統(tǒng)的噪聲來自熱漲落、散粒噪聲、讀出噪聲等多種源頭這些獨立因素的求和效應(yīng)會自發(fā)收斂到高斯分布。這就是中心極限定理的基本思想后面我會專門講。1.2 為什么計算機偏偏只給均勻分布你可能會問既然高斯分布這么重要為什么所有編程語言的隨機數(shù)接口不直接生成高斯分布這里有個歷史原因也有實現(xiàn)層面的原因。最底層的原因是計算機產(chǎn)生的是偽隨機數(shù)序列。無論用哪種算法本質(zhì)上都是從一個種子出發(fā)經(jīng)過一系列確定性數(shù)學(xué)運算生成一個在[0,1)區(qū)間內(nèi)均勻分布的序列。生成均勻分布本來就只需要讓這些序列“盡量均勻地鋪滿區(qū)間”判定標準很清晰。而高斯分布是無界的、形狀復(fù)雜沒法用簡單的線性同余之類的操作直接生成。所以標準做法是先用底層引擎生成均勻分布隨機數(shù)再做數(shù)學(xué)變換得到高斯分布。這個思路貫穿所有領(lǐng)域——Python里調(diào)用numpy.random.standard_normal底層用的也是這個邏輯C里std::normal_distribution也是。這樣做有個好處隨機數(shù)引擎和高斯變換是兩個獨立的模塊。引擎負責保證均勻隨機數(shù)的質(zhì)量和周期變換方法負責保證從均勻到高斯的映射正確。哪一邊出了問題都能單獨替換整個架構(gòu)非常干凈。我在做蒙特卡洛光線追跡時也習(xí)慣沿用這個分層思想先產(chǎn)生高質(zhì)量的均勻隨機數(shù)再根據(jù)物理模型做各種分布采樣絕不混在一起。2. 核心方法拆解Box-Muller變換的原理和證明2.1 Box-Muller變換一句公式解決大問題1958年Box和Muller發(fā)表了一篇簡短但影響深遠的論文給出了一個非常優(yōu)雅的結(jié)論如果U1和U2是相互獨立的均勻分布隨機數(shù)都滿足U(0,1)那么定義[ Z_0 \sqrt{-2\ln U_1}\cos(2\pi U_2) ] [ Z_1 \sqrt{-2\ln U_1}\sin(2\pi U_2) ]得到的Z0和Z1就是相互獨立的標準正態(tài)分布隨機數(shù)均值0、方差1。需要任意均值和標準差時再用公式Z μ σ * Z0做線性變換就行。這套公式第一次看到會覺得莫名其妙憑什么開個根號、乘個三角函數(shù)就變成高斯了我當時也困惑了好久直到我把推導(dǎo)過程完整走了一遍才真正理解。核心思路是把二維標準正態(tài)分布的聯(lián)合密度函數(shù)放到極坐標里看。二維標準正態(tài)分布的聯(lián)合密度是[ \frac{1}{2\pi} e^{-\frac{x^2y^2}{2}} ]這個函數(shù)只依賴x2y2也就是只依賴到原點的距離r。在極坐標下做變換x rcosθy rsinθ雅可比行列式給出了面積元從dxdy變成rdrdθ。于是分布可以拆成兩個獨立部分角度θ在[0, 2π)上均勻分布半徑平方R的定義要小心處理。具體來說令R X2 Y2。X和Y獨立且各服從標準正態(tài)分布時R服從自由度為2的卡方分布也就是參數(shù)為1/2的指數(shù)分布。而指數(shù)分布可以用逆變換采樣直接從均勻分布生成——如果U是U(0,1)均勻隨機數(shù)那么-2lnU就是參數(shù)為1/2的指數(shù)分布。這一下就把均勻隨機數(shù)U1和半徑R連起來了。角度θ本來就均勻分布直接取2πU2即可。再把極坐標換回直角坐標就有了上面的公式。理解了這個推導(dǎo)過程你就不會再“背公式背到懷疑人生”了。無非是高斯分布從極坐標看半徑服從指數(shù)分布角度均勻分布而指數(shù)分布恰好能用均勻分布逆變換生成。三個環(huán)節(jié)環(huán)環(huán)相扣。2.2 另一條路中心極限定理近似法除了Box-Muller變換還有一個流傳很廣的方法就是利用中心極限定理把12個獨立的U(0,1)均勻隨機數(shù)相加再減去6結(jié)果近似服從標準正態(tài)分布。為什么偏偏是12個因為單個U(0,1)均勻分布的均值為0.5、方差為1/12。12個獨立均勻分布之和均值是12×0.56方差是12×(1/12)1。這樣減6之后均值歸零、方差正好是1不需要額外的縮放系數(shù)。這個方法實現(xiàn)起來極其簡單我最早在單片機項目里生成高斯噪聲時就用的這個辦法因為MCU上跑浮點三角函數(shù)開銷不小加法卻很快。但是這個方法的缺點是尾巴很“禿”。12個[0,1)區(qū)間的數(shù)加起來最大就是12最小是0減6之后輸出的取值范圍嚴格落在[-6, 6]之間。而真正的標準正態(tài)分布理論上可以取到任意大的值雖然|Z|6的概率非常小約為十億分之一但在蒙特卡洛仿真里如果樣本量過億尾部事件就會開始影響結(jié)果。用中心極限定理生成的近似正態(tài)分布尾部是截斷的這對風險評估、極端情況分析這類場景是致命的。下表把兩種方法放在一起對比對比維度Box-Muller變換中心極限定理12個均勻相加精度精確服從正態(tài)分布近似尾部截斷計算開銷需要ln、cos、sin只需要12次加法和1次減法單次輸出數(shù)量每次生成2個獨立樣本每次生成1個樣本適合場景仿真精度要求高快速原型、嵌入式低算力環(huán)境易實現(xiàn)程度中等有邊界條件要處理非常簡單我個人的經(jīng)驗是除非是嵌入式環(huán)境實在不方便調(diào)用數(shù)學(xué)庫否則默認用Box-Muller或者它的改進版本。工程上求穩(wěn)精度不夠后面排查問題非常痛苦。3. 手寫代碼從Python到C的完整落地3.1 一段干凈的Box-Muller實現(xiàn)理論說了一堆代碼才是硬道理。下面是我用了很多年的Python實現(xiàn)注釋寫得比較詳細import math import random def box_muller_sample(): 用Box-Muller變換生成兩個獨立的標準正態(tài)分布隨機數(shù)。 返回: (z0, z1)均服從N(0, 1)。 # random.random() 返回 (0, 1] 區(qū)間有些實現(xiàn)是[0,1) # 注意必須嚴格大于0否則ln(0)會得到負無窮 u1 random.random() while u1 0.0: u1 random.random() u2 random.random() # 核心變換公式 mag math.sqrt(-2.0 * math.log(u1)) z0 mag * math.cos(2.0 * math.pi * u2) z1 mag * math.sin(2.0 * math.pi * u2) return z0, z1 def gaussian_sample(mu0.0, sigma1.0): 生成一個服從 N(mu, sigma^2) 的隨機數(shù)。 z0, _ box_muller_sample() return mu sigma * z0 # 驗證一下 if __name__ __main__: samples [gaussian_sample() for _ in range(100000)] mean sum(samples) / len(samples) var sum((x - mean) ** 2 for x in samples) / (len(samples) - 1) print(f均值: {mean:.4f}) print(f標準差: {math.sqrt(var):.4f})跑一下這段代碼輸出大致是這樣的均值: -0.0012 標準差: 0.9996在十萬個樣本量下均值和標準差都非常接近理論值0和1。偏差在0.01以內(nèi)是正常的畢竟是隨機抽樣存在天然的統(tǒng)計波動。如果你看到均值明顯偏離0比如達到0.05以上那就要懷疑隨機數(shù)質(zhì)量或者實現(xiàn)有沒有問題了。3.2 避免三角函數(shù)的極坐標法Marsaglia Polar MethodBox-Muller原始版本需要計算cos和sin這兩個函數(shù)在循環(huán)里調(diào)幾百萬次性能會很不好看。George Marsaglia在1962年提出一個改進版本用拒絕采樣繞開三角函數(shù)這就是極坐標法。算法思路很巧妙先在單位正方形內(nèi)隨機生成一個點(u, v)如果它落在單位圓內(nèi)u2v2 1就接受否則拒絕重來。然后利用這個點的坐標和半徑直接把角度信息藏在了坐標里不需要再用atan2或cos/sin去重建角度。import math import random def marsaglia_polar(): Marsaglia極坐標法生成兩個獨立標準正態(tài)隨機數(shù)。 不需要三角函數(shù)但可能需要多次生成(u,v)對。 while True: u random.uniform(-1.0, 1.0) v random.uniform(-1.0, 1.0) s u * u v * v if 0.0 s 1.0: break factor math.sqrt(-2.0 * math.log(s) / s) z0 u * factor z1 v * factor return z0, z1這個算法的拒絕率是多少呢單位正方形的面積是4內(nèi)切單位圓的面積是π所以隨機點落在圓內(nèi)的概率是π/4約78.5%。也就是說每生成一對(u,v)平均有21.5%的概率被拒絕需要再來一次。這個開銷遠小于三角函數(shù)計算的開銷實測下來整體速度比基礎(chǔ)版快30%以上。我在C/C項目里基本都用這個極坐標版本因為C標準庫的sin/cos依賴FPU高頻調(diào)用時性能波動明顯。如果你在做實時信號處理建議直接抄這個版本。3.3 用NumPy批量生成和驗證實際工程中很少一次只生成一兩個隨機數(shù)更多是要一整個數(shù)組。NumPy里可以直接用但為了驗證我們的Box-Muller實現(xiàn)也可以自己向量化import numpy as np def box_muller_batch(n): 用Box-Muller批量生成n個標準正態(tài)隨機數(shù)。 n為偶數(shù)時效率最高因為一次生成兩個。 n_half n // 2 u1 np.random.random(n_half) u2 np.random.random(n_half) # 防止log(0) u1 np.maximum(u1, np.finfo(float).eps) mag np.sqrt(-2.0 * np.log(u1)) z0 mag * np.cos(2.0 * np.pi * u2) z1 mag * np.sin(2.0 * np.pi * u2) result np.concatenate([z0, z1]) return result[:n] # 驗證分布形狀 data box_muller_batch(1000000) import matplotlib.pyplot as plt plt.hist(data, bins200, densityTrue, alpha0.7) # 畫出理論高斯曲線 x np.linspace(-4, 4, 500) y 1 / np.sqrt(2 * np.pi) * np.exp(-x**2 / 2) plt.plot(x, y, r-, linewidth2) plt.show()畫出來的直方圖和紅色理論曲線應(yīng)該幾乎完全重合。這種可視化驗證是判斷隨機數(shù)生成器容不容易出錯的最直觀方法比只看均值和方差靠譜多了——分布形狀是否正確、尾部是否對稱、有沒有明顯缺口一眼就能看出來。4. 工程場景實戰(zhàn)LightTools里的高斯分布設(shè)置4.1 光學(xué)仿真里的高斯分布從哪來光學(xué)仿真軟件里高斯分布出現(xiàn)得非常頻繁。激光二極管發(fā)出的光束其橫截面上的光強分布通常用高斯函數(shù)來描述這就是所謂的高斯光束模型。LED的配光曲線也經(jīng)常用高斯型分布來近似。在LightTools里做雜散光分析或者照明設(shè)計很多時候都需要設(shè)置光線的出射位置或者出射方向服從高斯分布。LightTools這類基于蒙特卡洛光線追跡的軟件本質(zhì)上做了大量隨機采樣。每一條光線的起點位置、發(fā)射方向、波長甚至表面反射的方向偏移都是靠隨機數(shù)決定的。如果采樣分布搞錯了追跡幾百萬條光線的結(jié)果也會整體跑偏而且這種錯誤非常隱蔽因為你從最終的照度圖上很難直接看出是分布參數(shù)設(shè)錯了還是仿真本身收斂不夠。4.2 LightTools中設(shè)置高斯分布的具體路徑不同版本的LightTools菜單位置略有差異但核心邏輯一脈相承。我以常用的設(shè)置方式說明在LightTools里進入光源屬性設(shè)置光源的發(fā)光特性里通常有“出射角度分布”或“強度分布”這樣的下拉選項。在下拉列表中選擇高斯分布后最關(guān)鍵的是設(shè)置兩個參數(shù)一個是分布的均值位置在角度分布里通常對應(yīng)0°也就是光軸中心方向另一個是標準差σ它決定了光束的角寬度。需要特別強調(diào)的是LightTools中的高斯分布參數(shù)絕大多數(shù)場景指的是“角度分布”而不是光源面的空間能量分布。角度分布的意思是光線出射方向相對于光軸的夾角θ其概率密度呈高斯分布。如果你設(shè)置σ10°那么大約68.3%的光線會落在偏離光軸±10°的范圍內(nèi)大約95.4%的光線落在±20°范圍內(nèi)。這個規(guī)律和標準高斯分布完全對應(yīng)。還有一個常用設(shè)置是光源面的空間強度分布。比如當你模擬一個高斯光束照射在接收面上時接收面上的輻照度分布是高斯型。LightTools里這類分布有時也被叫做“高斯輪廓”或者“自定義高斯型分布”配置方式同理會讓你輸入峰值位置和半寬參數(shù)。注意有些版本用的是半高全寬有些版本用的是1/e2寬度這個定義差異最容易讓人翻車。我自己的習(xí)慣是設(shè)置完后先在接收面上放一個探測器看實測的照度分布剖面確認一下半寬數(shù)值到底是按哪種定義算的。4.3 從均勻隨機數(shù)到高斯采樣的內(nèi)部邏輯LightTools內(nèi)部怎么把均勻隨機數(shù)變成高斯分布光線理解這一點對排查問題非常有幫助。它的底層思路和前面講的代碼一樣先用偽隨機數(shù)引擎生成均勻分布的隨機數(shù)序列再通過變換把它們映射到期望的分布上。光線從光源表面發(fā)射首先要決定發(fā)射點坐標。如果光源面是矩形坐標通常從均勻分布采樣然后決定發(fā)射方向如果發(fā)射方向要求高斯分布就會用Box-Muller變換或等價的查表法生成角度偏差。每一個這樣的采樣點對應(yīng)一條光線幾百萬條光線疊加起來就能統(tǒng)計出一個平滑的照度分布。所以你在LightTools里看到“光線數(shù)量”這個參數(shù)背后其實是一組隨機采樣序列的長度。光線數(shù)量太小高斯分布的統(tǒng)計漲落就會很明顯照度圖看起來毛躁不平滑。實際項目中我通常會用至少20萬條光線做初步仿真到了出圖驗證階段再用100萬條以上確保分布穩(wěn)定。4.4 參數(shù)設(shè)置案例與驗證步驟舉個具體例子。我在做一個激光照明系統(tǒng)的勻光設(shè)計時需要把激光二極管的快軸發(fā)散角模擬成高斯分布。激光二極管的快軸半高全寬大約30°對應(yīng)的標準差大約是12.7°半高全寬除以2.3548。在LightTools里新建一個光源把出射角度分布改為高斯分布均值設(shè)0°標準差設(shè)12.7°光線數(shù)量臨時設(shè)10萬條。在距離光源100mm的位置放一個接收面接收面尺寸覆蓋±50°發(fā)散角對應(yīng)的范圍。追跡完成后查看接收面的輻照度分布沿著x軸切一刀得到的輪廓應(yīng)該近似高斯鐘形曲線。如果輪廓偏平頂或者明顯不對稱多半是角度分布選項選成了均勻分布或者標準差定義換算錯了。這個驗證步驟很值得養(yǎng)成習(xí)慣任何光源模型改動之后先花十分鐘做個簡單的正向驗證確認分布形態(tài)正確再跑完整的系統(tǒng)仿真。否則幾小時的追跡結(jié)果可能全部作廢。5. 實操中踩過的坑常見問題與排查技巧5.1 生成的序列“不那么高斯”是怎么回事表格整理我這些年遇到的高頻問題現(xiàn)象可能原因排查思路均值偏離目標值很大變換公式寫錯或邊界值沒處理用幾組U1、U2手算驗證或者畫直方圖看分布中心方差偏小采樣時用了有偏方法或隨機數(shù)序列周期太短檢查是否誤用了CLT近似加大樣本量直方圖左右不對稱隨機數(shù)引擎質(zhì)量差或變換中用了截斷換引擎測試檢查是否有while循環(huán)誤截斷生成速度太慢循環(huán)里反復(fù)調(diào)用三角函數(shù)或每次只生成一個改用Marsaglia極坐標法或批量生成出現(xiàn)NaN或inf輸入U1為0log(0)導(dǎo)致無窮在采樣函數(shù)里加邊界判斷確保U1 05.2 邊界條件一個零值引發(fā)的血案之前我在一個C語言模塊里實現(xiàn)Box-Muller測試時偶爾冒出NaN。追了半天發(fā)現(xiàn)是最底層的均勻隨機數(shù)生成器偶爾返回精確的0.0。log(0)等于負無窮sqrt(負無窮)直接得到NaN。解決方案有兩種。最簡單的是在采樣前做一個保護判斷如果U1等于0就重新采樣一次。因為連續(xù)均勻分布取到精確0的概率微乎其微重新采一次幾乎不可能再次為0。另一種方案是用U11-U1做變換把(0,1]區(qū)間的值映射到[0,1)區(qū)間再取一個極小值做上下限夾逼。我推薦第一種邏輯簡單、不引入額外偏差。5.3 隨機數(shù)質(zhì)量對結(jié)果的影響很多人沒意識到偽隨機數(shù)生成器的質(zhì)量會直接影響高斯樣本的質(zhì)量。早期C語言的rand()函數(shù)周期短、低位隨機性差用的時候會發(fā)現(xiàn)生成的高斯序列在高位和低位分布不均勻?，F(xiàn)代推薦用PCG或者Mersenne Twister這類經(jīng)過驗證的引擎。Python的random模塊底層是Mersenne Twister一般夠用NumPy從1.17版本開始默認用PCG64質(zhì)量更好。C里std::mt19937也是成熟選擇。有個判斷隨機數(shù)質(zhì)量的小技巧生成一批高斯樣本后算一下樣本的四分位數(shù)和理論標準正態(tài)分布的四分位數(shù)對比。如果偏差持續(xù)超過幾個百分比就要懷疑引擎了。另外可以做自相關(guān)檢查看看生成的序列里有沒有周期性規(guī)律——正規(guī)的高斯白噪聲自相關(guān)系數(shù)應(yīng)該幾乎為零。5.4 大規(guī)模生成時的性能優(yōu)化方向當隨機數(shù)需求膨脹到千萬甚至億級時Box-Muller就算不上最優(yōu)解了。這時業(yè)界常用的是Ziggurat算法它用拒絕采樣和預(yù)計算查找表的方式把生成成本壓縮到每次僅需一次比較和一次查表速度比Box-Muller快2到4倍。不過Ziggurat的實現(xiàn)復(fù)雜度更高需要精心預(yù)計算表格。如果沒有極端的性能要求我建議先用極坐標法畢竟維護起來省心得多。另外可以做的優(yōu)化是向量化。在Python里用NumPy一次性生成上百萬個u1和u2數(shù)組利用底層C實現(xiàn)的向量化運算比for循環(huán)逐個生成快幾個數(shù)量級。我之前把一個Python循環(huán)版本改成向量化版本十萬個樣本的生成時間從1.2秒左右降到了毫秒級別。5.5 一個小技巧直接用Box-Muller生成二維高斯光斑采樣最后分享一個工程上很實用的技巧。在做光學(xué)仿真前處理時經(jīng)常需要在一個圓形光斑內(nèi)生成服從高斯分布的采樣點坐標。這時可以直接利用Box-Muller生成的z0和z1兩個獨立標準正態(tài)隨機數(shù)把它們直接當作x、y坐標使用。因為二維標準正態(tài)分布的等概率密度線就是同心圓聯(lián)合分布天然是中心對稱的圓形高斯光斑。如果你想要半高全寬可控的光斑只需要做縮放x FWHM / 2.3548 * z0y FWHM / 2.3548 * z1。這樣生成的坐標點自然形成高斯圓形彌散斑。相比先均勻生成半徑和角度再變換的方法這個做法不需要計算反正切代碼更簡潔分布也精確。我把這個函數(shù)封裝在自己的工具庫里凡是需要模擬高斯光斑的地方都直接調(diào)它。我個人在實際項目中的體會是均勻分布到高斯分布的變換看起來只是一個公式的事但越深入就越發(fā)現(xiàn)它連接著概率論、數(shù)值計算、仿真工程好幾個層面的知識。這些原理性的東西一旦吃透了在LightTools、Zemax等軟件里遇到分布相關(guān)的設(shè)置時就不會再犯迷糊——因為你一眼就能看出來軟件底層在做什么數(shù)學(xué)操作。這大概就是“底層原理”和“工具使用”之間最有趣的關(guān)系工具會過時但原理不會。 均勻分布生成高斯分布：從Box-Muller到LightTools實戰(zhàn)
做光學(xué)仿真和隨機模擬這些年我發(fā)現(xiàn)自己繞不開一個坎所有編程語言和仿真軟件能直接生成的隨機數(shù)幾乎都是均勻分布。可現(xiàn)實世界里真正常用的是高斯分布——也就是正態(tài)分布。測量噪聲是高斯分布光斑的能量分布接近高斯分布人群身高的統(tǒng)計也是高斯分布。于是“均勻分布產(chǎn)生高斯分布”就成了一個高頻問題網(wǎng)上搜一下相關(guān)討論特別多連LightTools這種光學(xué)仿真軟件里怎么設(shè)置高斯分布都被反復(fù)問。這篇文章我打算把這件事徹底講透從數(shù)學(xué)原理講到代碼實現(xiàn)再落到LightTools這類工程工具里的實際操作把我踩過的坑和驗證過的方法都整理出來。1. 均勻分布和高斯分布先搞清楚我們要干什么1.1 兩個分布到底差在哪里均勻分布的概率密度函數(shù)是一條水平直線在定義區(qū)間內(nèi)每個點出現(xiàn)的概率一樣。打個比方均勻分布就像抽簽箱子里十個球抽中任何一個的概率都是十分之一。而高斯分布是一條鐘形曲線中間高、兩邊低絕大多數(shù)樣本落在均值附近極端值幾乎不會出現(xiàn)。它的概率密度函數(shù)長這樣[ f(x) \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{(x-\mu)^2}{2\sigma^2}} ]這里的μ是均值σ是標準差σ2是方差。μ決定了鐘形曲線在x軸上的位置σ決定了曲線是“胖”還是“瘦”。σ越小曲線越尖數(shù)據(jù)越集中σ越大曲線越平數(shù)據(jù)越分散。高中數(shù)學(xué)里大家可能背過這個公式但沒有多少人認真想過它背后的幾何含義。高斯分布之所以無處不在本質(zhì)上是因為自然界里大多數(shù)“誤差”和“波動”都是大量微小獨立因素疊加的結(jié)果——一個人的身高受幾百個基因位點影響光學(xué)系統(tǒng)的噪聲來自熱漲落、散粒噪聲、讀出噪聲等多種源頭這些獨立因素的求和效應(yīng)會自發(fā)收斂到高斯分布。這就是中心極限定理的基本思想后面我會專門講。1.2 為什么計算機偏偏只給均勻分布你可能會問既然高斯分布這么重要為什么所有編程語言的隨機數(shù)接口不直接生成高斯分布這里有個歷史原因也有實現(xiàn)層面的原因。最底層的原因是計算機產(chǎn)生的是偽隨機數(shù)序列。無論用哪種算法本質(zhì)上都是從一個種子出發(fā)經(jīng)過一系列確定性數(shù)學(xué)運算生成一個在[0,1)區(qū)間內(nèi)均勻分布的序列。生成均勻分布本來就只需要讓這些序列“盡量均勻地鋪滿區(qū)間”判定標準很清晰。而高斯分布是無界的、形狀復(fù)雜沒法用簡單的線性同余之類的操作直接生成。所以標準做法是先用底層引擎生成均勻分布隨機數(shù)再做數(shù)學(xué)變換得到高斯分布。這個思路貫穿所有領(lǐng)域——Python里調(diào)用numpy.random.standard_normal底層用的也是這個邏輯C里std::normal_distribution也是。這樣做有個好處隨機數(shù)引擎和高斯變換是兩個獨立的模塊。引擎負責保證均勻隨機數(shù)的質(zhì)量和周期變換方法負責保證從均勻到高斯的映射正確。哪一邊出了問題都能單獨替換整個架構(gòu)非常干凈。我在做蒙特卡洛光線追跡時也習(xí)慣沿用這個分層思想先產(chǎn)生高質(zhì)量的均勻隨機數(shù)再根據(jù)物理模型做各種分布采樣絕不混在一起。2. 核心方法拆解Box-Muller變換的原理和證明2.1 Box-Muller變換一句公式解決大問題1958年Box和Muller發(fā)表了一篇簡短但影響深遠的論文給出了一個非常優(yōu)雅的結(jié)論如果U1和U2是相互獨立的均勻分布隨機數(shù)都滿足U(0,1)那么定義[ Z_0 \sqrt{-2\ln U_1}\cos(2\pi U_2) ] [ Z_1 \sqrt{-2\ln U_1}\sin(2\pi U_2) ]得到的Z0和Z1就是相互獨立的標準正態(tài)分布隨機數(shù)均值0、方差1。需要任意均值和標準差時再用公式Z μ σ * Z0做線性變換就行。這套公式第一次看到會覺得莫名其妙憑什么開個根號、乘個三角函數(shù)就變成高斯了我當時也困惑了好久直到我把推導(dǎo)過程完整走了一遍才真正理解。核心思路是把二維標準正態(tài)分布的聯(lián)合密度函數(shù)放到極坐標里看。二維標準正態(tài)分布的聯(lián)合密度是[ \frac{1}{2\pi} e^{-\frac{x^2y^2}{2}} ]這個函數(shù)只依賴x2y2也就是只依賴到原點的距離r。在極坐標下做變換x rcosθy rsinθ雅可比行列式給出了面積元從dxdy變成rdrdθ。于是分布可以拆成兩個獨立部分角度θ在[0, 2π)上均勻分布半徑平方R的定義要小心處理。具體來說令R X2 Y2。X和Y獨立且各服從標準正態(tài)分布時R服從自由度為2的卡方分布也就是參數(shù)為1/2的指數(shù)分布。而指數(shù)分布可以用逆變換采樣直接從均勻分布生成——如果U是U(0,1)均勻隨機數(shù)那么-2lnU就是參數(shù)為1/2的指數(shù)分布。這一下就把均勻隨機數(shù)U1和半徑R連起來了。角度θ本來就均勻分布直接取2πU2即可。再把極坐標換回直角坐標就有了上面的公式。理解了這個推導(dǎo)過程你就不會再“背公式背到懷疑人生”了。無非是高斯分布從極坐標看半徑服從指數(shù)分布角度均勻分布而指數(shù)分布恰好能用均勻分布逆變換生成。三個環(huán)節(jié)環(huán)環(huán)相扣。2.2 另一條路中心極限定理近似法除了Box-Muller變換還有一個流傳很廣的方法就是利用中心極限定理把12個獨立的U(0,1)均勻隨機數(shù)相加再減去6結(jié)果近似服從標準正態(tài)分布。為什么偏偏是12個因為單個U(0,1)均勻分布的均值為0.5、方差為1/12。12個獨立均勻分布之和均值是12×0.56方差是12×(1/12)1。這樣減6之后均值歸零、方差正好是1不需要額外的縮放系數(shù)。這個方法實現(xiàn)起來極其簡單我最早在單片機項目里生成高斯噪聲時就用的這個辦法因為MCU上跑浮點三角函數(shù)開銷不小加法卻很快。但是這個方法的缺點是尾巴很“禿”。12個[0,1)區(qū)間的數(shù)加起來最大就是12最小是0減6之后輸出的取值范圍嚴格落在[-6, 6]之間。而真正的標準正態(tài)分布理論上可以取到任意大的值雖然|Z|6的概率非常小約為十億分之一但在蒙特卡洛仿真里如果樣本量過億尾部事件就會開始影響結(jié)果。用中心極限定理生成的近似正態(tài)分布尾部是截斷的這對風險評估、極端情況分析這類場景是致命的。下表把兩種方法放在一起對比對比維度Box-Muller變換中心極限定理12個均勻相加精度精確服從正態(tài)分布近似尾部截斷計算開銷需要ln、cos、sin只需要12次加法和1次減法單次輸出數(shù)量每次生成2個獨立樣本每次生成1個樣本適合場景仿真精度要求高快速原型、嵌入式低算力環(huán)境易實現(xiàn)程度中等有邊界條件要處理非常簡單我個人的經(jīng)驗是除非是嵌入式環(huán)境實在不方便調(diào)用數(shù)學(xué)庫否則默認用Box-Muller或者它的改進版本。工程上求穩(wěn)精度不夠后面排查問題非常痛苦。3. 手寫代碼從Python到C的完整落地3.1 一段干凈的Box-Muller實現(xiàn)理論說了一堆代碼才是硬道理。下面是我用了很多年的Python實現(xiàn)注釋寫得比較詳細import math import random def box_muller_sample(): 用Box-Muller變換生成兩個獨立的標準正態(tài)分布隨機數(shù)。 返回: (z0, z1)均服從N(0, 1)。 # random.random() 返回 (0, 1] 區(qū)間有些實現(xiàn)是[0,1) # 注意必須嚴格大于0否則ln(0)會得到負無窮 u1 random.random() while u1 0.0: u1 random.random() u2 random.random() # 核心變換公式 mag math.sqrt(-2.0 * math.log(u1)) z0 mag * math.cos(2.0 * math.pi * u2) z1 mag * math.sin(2.0 * math.pi * u2) return z0, z1 def gaussian_sample(mu0.0, sigma1.0): 生成一個服從 N(mu, sigma^2) 的隨機數(shù)。 z0, _ box_muller_sample() return mu sigma * z0 # 驗證一下 if __name__ __main__: samples [gaussian_sample() for _ in range(100000)] mean sum(samples) / len(samples) var sum((x - mean) ** 2 for x in samples) / (len(samples) - 1) print(f均值: {mean:.4f}) print(f標準差: {math.sqrt(var):.4f})跑一下這段代碼輸出大致是這樣的均值: -0.0012 標準差: 0.9996在十萬個樣本量下均值和標準差都非常接近理論值0和1。偏差在0.01以內(nèi)是正常的畢竟是隨機抽樣存在天然的統(tǒng)計波動。如果你看到均值明顯偏離0比如達到0.05以上那就要懷疑隨機數(shù)質(zhì)量或者實現(xiàn)有沒有問題了。3.2 避免三角函數(shù)的極坐標法Marsaglia Polar MethodBox-Muller原始版本需要計算cos和sin這兩個函數(shù)在循環(huán)里調(diào)幾百萬次性能會很不好看。George Marsaglia在1962年提出一個改進版本用拒絕采樣繞開三角函數(shù)這就是極坐標法。算法思路很巧妙先在單位正方形內(nèi)隨機生成一個點(u, v)如果它落在單位圓內(nèi)u2v2 1就接受否則拒絕重來。然后利用這個點的坐標和半徑直接把角度信息藏在了坐標里不需要再用atan2或cos/sin去重建角度。import math import random def marsaglia_polar(): Marsaglia極坐標法生成兩個獨立標準正態(tài)隨機數(shù)。 不需要三角函數(shù)但可能需要多次生成(u,v)對。 while True: u random.uniform(-1.0, 1.0) v random.uniform(-1.0, 1.0) s u * u v * v if 0.0 s 1.0: break factor math.sqrt(-2.0 * math.log(s) / s) z0 u * factor z1 v * factor return z0, z1這個算法的拒絕率是多少呢單位正方形的面積是4內(nèi)切單位圓的面積是π所以隨機點落在圓內(nèi)的概率是π/4約78.5%。也就是說每生成一對(u,v)平均有21.5%的概率被拒絕需要再來一次。這個開銷遠小于三角函數(shù)計算的開銷實測下來整體速度比基礎(chǔ)版快30%以上。我在C/C項目里基本都用這個極坐標版本因為C標準庫的sin/cos依賴FPU高頻調(diào)用時性能波動明顯。如果你在做實時信號處理建議直接抄這個版本。3.3 用NumPy批量生成和驗證實際工程中很少一次只生成一兩個隨機數(shù)更多是要一整個數(shù)組。NumPy里可以直接用但為了驗證我們的Box-Muller實現(xiàn)也可以自己向量化import numpy as np def box_muller_batch(n): 用Box-Muller批量生成n個標準正態(tài)隨機數(shù)。 n為偶數(shù)時效率最高因為一次生成兩個。 n_half n // 2 u1 np.random.random(n_half) u2 np.random.random(n_half) # 防止log(0) u1 np.maximum(u1, np.finfo(float).eps) mag np.sqrt(-2.0 * np.log(u1)) z0 mag * np.cos(2.0 * np.pi * u2) z1 mag * np.sin(2.0 * np.pi * u2) result np.concatenate([z0, z1]) return result[:n] # 驗證分布形狀 data box_muller_batch(1000000) import matplotlib.pyplot as plt plt.hist(data, bins200, densityTrue, alpha0.7) # 畫出理論高斯曲線 x np.linspace(-4, 4, 500) y 1 / np.sqrt(2 * np.pi) * np.exp(-x**2 / 2) plt.plot(x, y, r-, linewidth2) plt.show()畫出來的直方圖和紅色理論曲線應(yīng)該幾乎完全重合。這種可視化驗證是判斷隨機數(shù)生成器容不容易出錯的最直觀方法比只看均值和方差靠譜多了——分布形狀是否正確、尾部是否對稱、有沒有明顯缺口一眼就能看出來。4. 工程場景實戰(zhàn)LightTools里的高斯分布設(shè)置4.1 光學(xué)仿真里的高斯分布從哪來光學(xué)仿真軟件里高斯分布出現(xiàn)得非常頻繁。激光二極管發(fā)出的光束其橫截面上的光強分布通常用高斯函數(shù)來描述這就是所謂的高斯光束模型。LED的配光曲線也經(jīng)常用高斯型分布來近似。在LightTools里做雜散光分析或者照明設(shè)計很多時候都需要設(shè)置光線的出射位置或者出射方向服從高斯分布。LightTools這類基于蒙特卡洛光線追跡的軟件本質(zhì)上做了大量隨機采樣。每一條光線的起點位置、發(fā)射方向、波長甚至表面反射的方向偏移都是靠隨機數(shù)決定的。如果采樣分布搞錯了追跡幾百萬條光線的結(jié)果也會整體跑偏而且這種錯誤非常隱蔽因為你從最終的照度圖上很難直接看出是分布參數(shù)設(shè)錯了還是仿真本身收斂不夠。4.2 LightTools中設(shè)置高斯分布的具體路徑不同版本的LightTools菜單位置略有差異但核心邏輯一脈相承。我以常用的設(shè)置方式說明在LightTools里進入光源屬性設(shè)置光源的發(fā)光特性里通常有“出射角度分布”或“強度分布”這樣的下拉選項。在下拉列表中選擇高斯分布后最關(guān)鍵的是設(shè)置兩個參數(shù)一個是分布的均值位置在角度分布里通常對應(yīng)0°也就是光軸中心方向另一個是標準差σ它決定了光束的角寬度。需要特別強調(diào)的是LightTools中的高斯分布參數(shù)絕大多數(shù)場景指的是“角度分布”而不是光源面的空間能量分布。角度分布的意思是光線出射方向相對于光軸的夾角θ其概率密度呈高斯分布。如果你設(shè)置σ10°那么大約68.3%的光線會落在偏離光軸±10°的范圍內(nèi)大約95.4%的光線落在±20°范圍內(nèi)。這個規(guī)律和標準高斯分布完全對應(yīng)。還有一個常用設(shè)置是光源面的空間強度分布。比如當你模擬一個高斯光束照射在接收面上時接收面上的輻照度分布是高斯型。LightTools里這類分布有時也被叫做“高斯輪廓”或者“自定義高斯型分布”配置方式同理會讓你輸入峰值位置和半寬參數(shù)。注意有些版本用的是半高全寬有些版本用的是1/e2寬度這個定義差異最容易讓人翻車。我自己的習(xí)慣是設(shè)置完后先在接收面上放一個探測器看實測的照度分布剖面確認一下半寬數(shù)值到底是按哪種定義算的。4.3 從均勻隨機數(shù)到高斯采樣的內(nèi)部邏輯LightTools內(nèi)部怎么把均勻隨機數(shù)變成高斯分布光線理解這一點對排查問題非常有幫助。它的底層思路和前面講的代碼一樣先用偽隨機數(shù)引擎生成均勻分布的隨機數(shù)序列再通過變換把它們映射到期望的分布上。光線從光源表面發(fā)射首先要決定發(fā)射點坐標。如果光源面是矩形坐標通常從均勻分布采樣然后決定發(fā)射方向如果發(fā)射方向要求高斯分布就會用Box-Muller變換或等價的查表法生成角度偏差。每一個這樣的采樣點對應(yīng)一條光線幾百萬條光線疊加起來就能統(tǒng)計出一個平滑的照度分布。所以你在LightTools里看到“光線數(shù)量”這個參數(shù)背后其實是一組隨機采樣序列的長度。光線數(shù)量太小高斯分布的統(tǒng)計漲落就會很明顯照度圖看起來毛躁不平滑。實際項目中我通常會用至少20萬條光線做初步仿真到了出圖驗證階段再用100萬條以上確保分布穩(wěn)定。4.4 參數(shù)設(shè)置案例與驗證步驟舉個具體例子。我在做一個激光照明系統(tǒng)的勻光設(shè)計時需要把激光二極管的快軸發(fā)散角模擬成高斯分布。激光二極管的快軸半高全寬大約30°對應(yīng)的標準差大約是12.7°半高全寬除以2.3548。在LightTools里新建一個光源把出射角度分布改為高斯分布均值設(shè)0°標準差設(shè)12.7°光線數(shù)量臨時設(shè)10萬條。在距離光源100mm的位置放一個接收面接收面尺寸覆蓋±50°發(fā)散角對應(yīng)的范圍。追跡完成后查看接收面的輻照度分布沿著x軸切一刀得到的輪廓應(yīng)該近似高斯鐘形曲線。如果輪廓偏平頂或者明顯不對稱多半是角度分布選項選成了均勻分布或者標準差定義換算錯了。這個驗證步驟很值得養(yǎng)成習(xí)慣任何光源模型改動之后先花十分鐘做個簡單的正向驗證確認分布形態(tài)正確再跑完整的系統(tǒng)仿真。否則幾小時的追跡結(jié)果可能全部作廢。5. 實操中踩過的坑常見問題與排查技巧5.1 生成的序列“不那么高斯”是怎么回事表格整理我這些年遇到的高頻問題現(xiàn)象可能原因排查思路均值偏離目標值很大變換公式寫錯或邊界值沒處理用幾組U1、U2手算驗證或者畫直方圖看分布中心方差偏小采樣時用了有偏方法或隨機數(shù)序列周期太短檢查是否誤用了CLT近似加大樣本量直方圖左右不對稱隨機數(shù)引擎質(zhì)量差或變換中用了截斷換引擎測試檢查是否有while循環(huán)誤截斷生成速度太慢循環(huán)里反復(fù)調(diào)用三角函數(shù)或每次只生成一個改用Marsaglia極坐標法或批量生成出現(xiàn)NaN或inf輸入U1為0log(0)導(dǎo)致無窮在采樣函數(shù)里加邊界判斷確保U1 05.2 邊界條件一個零值引發(fā)的血案之前我在一個C語言模塊里實現(xiàn)Box-Muller測試時偶爾冒出NaN。追了半天發(fā)現(xiàn)是最底層的均勻隨機數(shù)生成器偶爾返回精確的0.0。log(0)等于負無窮sqrt(負無窮)直接得到NaN。解決方案有兩種。最簡單的是在采樣前做一個保護判斷如果U1等于0就重新采樣一次。因為連續(xù)均勻分布取到精確0的概率微乎其微重新采一次幾乎不可能再次為0。另一種方案是用U11-U1做變換把(0,1]區(qū)間的值映射到[0,1)區(qū)間再取一個極小值做上下限夾逼。我推薦第一種邏輯簡單、不引入額外偏差。5.3 隨機數(shù)質(zhì)量對結(jié)果的影響很多人沒意識到偽隨機數(shù)生成器的質(zhì)量會直接影響高斯樣本的質(zhì)量。早期C語言的rand()函數(shù)周期短、低位隨機性差用的時候會發(fā)現(xiàn)生成的高斯序列在高位和低位分布不均勻?，F(xiàn)代推薦用PCG或者Mersenne Twister這類經(jīng)過驗證的引擎。Python的random模塊底層是Mersenne Twister一般夠用NumPy從1.17版本開始默認用PCG64質(zhì)量更好。C里std::mt19937也是成熟選擇。有個判斷隨機數(shù)質(zhì)量的小技巧生成一批高斯樣本后算一下樣本的四分位數(shù)和理論標準正態(tài)分布的四分位數(shù)對比。如果偏差持續(xù)超過幾個百分比就要懷疑引擎了。另外可以做自相關(guān)檢查看看生成的序列里有沒有周期性規(guī)律——正規(guī)的高斯白噪聲自相關(guān)系數(shù)應(yīng)該幾乎為零。5.4 大規(guī)模生成時的性能優(yōu)化方向當隨機數(shù)需求膨脹到千萬甚至億級時Box-Muller就算不上最優(yōu)解了。這時業(yè)界常用的是Ziggurat算法它用拒絕采樣和預(yù)計算查找表的方式把生成成本壓縮到每次僅需一次比較和一次查表速度比Box-Muller快2到4倍。不過Ziggurat的實現(xiàn)復(fù)雜度更高需要精心預(yù)計算表格。如果沒有極端的性能要求我建議先用極坐標法畢竟維護起來省心得多。另外可以做的優(yōu)化是向量化。在Python里用NumPy一次性生成上百萬個u1和u2數(shù)組利用底層C實現(xiàn)的向量化運算比for循環(huán)逐個生成快幾個數(shù)量級。我之前把一個Python循環(huán)版本改成向量化版本十萬個樣本的生成時間從1.2秒左右降到了毫秒級別。5.5 一個小技巧直接用Box-Muller生成二維高斯光斑采樣最后分享一個工程上很實用的技巧。在做光學(xué)仿真前處理時經(jīng)常需要在一個圓形光斑內(nèi)生成服從高斯分布的采樣點坐標。這時可以直接利用Box-Muller生成的z0和z1兩個獨立標準正態(tài)隨機數(shù)把它們直接當作x、y坐標使用。因為二維標準正態(tài)分布的等概率密度線就是同心圓聯(lián)合分布天然是中心對稱的圓形高斯光斑。如果你想要半高全寬可控的光斑只需要做縮放x FWHM / 2.3548 * z0y FWHM / 2.3548 * z1。這樣生成的坐標點自然形成高斯圓形彌散斑。相比先均勻生成半徑和角度再變換的方法這個做法不需要計算反正切代碼更簡潔分布也精確。我把這個函數(shù)封裝在自己的工具庫里凡是需要模擬高斯光斑的地方都直接調(diào)它。我個人在實際項目中的體會是均勻分布到高斯分布的變換看起來只是一個公式的事但越深入就越發(fā)現(xiàn)它連接著概率論、數(shù)值計算、仿真工程好幾個層面的知識。這些原理性的東西一旦吃透了在LightTools、Zemax等軟件里遇到分布相關(guān)的設(shè)置時就不會再犯迷糊——因為你一眼就能看出來軟件底層在做什么數(shù)學(xué)操作。這大概就是“底層原理”和“工具使用”之間最有趣的關(guān)系工具會過時但原理不會。