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