Box-Muller算法：从均匀分布随机数到高斯分布的MATLAB实现

简介：面向统计计算与随机模拟场景，这份MATLAB资源聚焦如何将均匀分布随机数转换为高斯分布（正态分布）随机数，适合需要掌握随机数生成原理的科研人员、工程师及在校学生。资源围绕Box-Muller变换与12法则两种经典方法展开，对应实现代码可直接运行或二次修改，帮助理解正态分布模拟的底层逻辑。包内共4个文件，以3个.m脚本文件为核心，分别提供均匀分布转换、高斯随机数生成等函数示例，另含1个txt说明文档，便于对照查阅；压缩包整体仅1KB，轻量易用。目前已有399人学习浏览，适合快速入门随机数生成与MATLAB数值实验。通过学习本包，读者可以获得完整可用的转换函数源码、方法原理注释以及实际调用示例，能够直接嵌入自己的仿真或统计项目中，提升建模效率与结果可信度。

1. Box-Muller 生成高斯分布：把 uniformgauss 从 randn 的黑匣子里拆出来

自己做 matlab 仿真时，高斯分布随机数通常一句话就来了：randn。可一旦你需要把随机数生成过程写进论文、移植到硬件仿真器，或者复现教科书里的粒子滤波，randn就成了黑匣子。Box-Muller 是上世纪五十年代提出的经典算法，输入两个独立均匀分布随机数，输出两个独立高斯分布随机数。很多叫uniformgauss的资源包里装的就是这个东西的 matlab 封装，用来做高斯分布 MATLAB 仿真、图像加噪、通信系统误码率统计。这篇文章直接拆掉黑匣子：先立住原理，再给出能跑的向量化代码，最后把四个最容易翻车的边界坑和验证手段一起讲清楚，让你拿到手就能复现并敢在自己的工程里用。

2. 从均匀到高斯：Box-Muller 原理、两种实现和为什么不用逆变换

2.1 Box-Muller 公式在说什么：两次均匀采样换一对高斯样本

Box-Muller 的核心思想是坐标变换。我们盯着二维标准正态分布看：它的概率密度是

f(x,y) = (1/2π) * exp(-(x²+y²)/2)

这个式子只依赖r² = x²+y²，和角度无关。换到极坐标后，半径和角度可以拆成两个独立问题。设半径平方R = r²，做一次变量替换，能得到R服从参数为1/2的指数分布，而角度Θ服从[0,2π)上的均匀分布。指数分布恰好可以用逆变换法生成：取均匀随机数U1，令R = -2 ln(U1)，再取第二个均匀随机数U2生成角度Θ = 2π U2。于是：

Z1 = sqrt(-2 ln U1) * cos(2π U2)Z2 = sqrt(-2 ln U1) * sin(2π U2)

这就是 Box-Muller 变换。注意它一次性生成两个独立样本，不是生成一个。很多人图省事只取cos分量，等于把一半信息丢掉了，统计性质会打折扣。我在工程里见到过把这个公式抄错的情况，最典型的是把-2 ln U1写成-ln U1，导致方差少了一倍。公式里的sqrt(-2 ln U1)本质是把指数分布变成标准正态分布的半径，少了系数2，尾部就变瘦了。

在 Matlab 里实现这段逻辑特别直接，因为rand本身就给出均匀分布。需要注意U1不能取到 0，否则ln(0)是负无穷，后面会得到 NaN。Matlab 的rand在大多数版本里返回(0,1)开区间样本，但为了在旧版或混合引擎下更稳，我习惯在代码里主动做一次max(rand, eps)保护。这个细节放到避坑章细说。

2.2 基本形式与 Marsaglia 极坐标法：一个看中简单，一个看中速度

直接用cos和sin的 Box-Muller 叫基本形式，也叫三角函数形式。它的优点是代码和公式一一对应，几乎不可能写错。缺点是三角函数算得慢，尤其是在老版本 Matlab 里跑循环，几百万样本时能明显感受到卡顿。

另一个流派是 Marsaglia 极坐标法，也叫 Polar 方法。它绕开三角函数：先产生两个在[-1,1]上的均匀随机数V1, V2，计算S = V1² + V2²；如果S落在(0,1)内就接受，否则重来。然后用下面公式生成两个标准正态样本：

Z1 = V1 * sqrt(-2 ln(S) / S)Z2 = V2 * sqrt(-2 ln(S) / S)

这个做法的妙处在于S本身就是极坐标里的半径平方，而且从[-1,1]均匀分布中抽取落到单位圆内的概率是π/4，约等于 0.785，拒绝率不算高。换来的好处是整个生成过程只有四则运算和一次log，没有三角函数，在循环中更友好。

我用 Matlab 做对比时发现，向量化之后基本形式其实更快，因为cos和sin内置优化得很好。但极坐标法在产生单个随机数时更有优势，因为你可以预生成一批候选者。到底选哪个？我的习惯是：教学演示、代码可读性优先时用基本形式；性能敏感且要嵌入while循环时用极坐标法。

下面给极坐标法的典型循环实现：

function z = marsaglia_polar(n) z = zeros(1, n); i = 1; while i <= n v1 = 2 * rand - 1; v2 = 2 * rand - 1; s = v1^2 + v2^2; if s > 0 && s < 1 factor = sqrt(-2 * log(s) / s); z(i) = v1 * factor; if i + 1 <= n z(i+1) = v2 * factor; i = i + 1; end i = i + 1; end end end

这里的关键参数是v1和v2的取值范围。rand默认返回[0,1)，一定记得先平移缩放成[-1,1)。s > 0是为了防止log(0)，s < 1是极坐标法的接受条件。这个算法一次循环可能生成 0 个、1 个或 2 个样本，所以while的终止条件要仔细检查索引。我在早期写的时候就是因为漏了s > 0，最终在极小概率下得到 NaN。

2.3 为什么逆变换法在 Matlab 里不常用：高斯 CDF 没有解析逆

常见的随机数生成思路是逆变换采样：如果随机变量X的累积分布函数是F(x)，那么F(X)服从均匀分布，反过来F^{-1}(U)也就服从原分布。这条路线对指数分布、均匀分布都很好用，因为它们的F^{-1}是显式表达式。但高斯分布的F(x)是误差函数erf的变形，没有解析逆，只能依赖数值迭代求反。

Matlab 里确实提供了erfinv，理论上可以用sqrt(2) * erfinv(2*U-1)生成标准正态样本。实际测量下来，它的速度大约比 Box-Muller 慢 5 到 10 倍，而且数值误差在尾部附近会被放大。更关键的是，erfinv在U非常接近 0 或 1 时产生极大或极小的值，这会给仿真带来不稳定。Box-Muller 之所以是教科书默认方案，就是因为它用对数函数和三角函数绕开了高斯 CDF 求逆，同时还能成对产出样本。明白了这一点，你在选择生成策略时就不会被"逆变换万能"这句话带偏。

3. 用 Matlab 把 uniformgauss 写成可复用函数：向量化代码与两个必调参数

3.1 最小可用的 box_muller.m：参数 n、mu、sigma

很多uniformgauss压缩包里的核心代码就是一个函数：输入样本数量、均值和标准差，输出一段高斯分布随机数。我自己写的 min 版本长这样：

function x = box_muller(n, mu, sigma) if nargin < 2 mu = 0; end if nargin < 3 sigma = 1; end if sigma <= 0 error('sigma 必须大于 0'); end u1 = max(rand(1, n), eps); u2 = rand(1, n); z1 = sqrt(-2 * log(u1)) .* cos(2 * pi * u2); z2 = sqrt(-2 * log(u1)) .* sin(2 * pi * u2); x = [z1, z2] * sigma + mu; x = x(1:n); end

逻辑说明：n是请求的样本数，但我一次性生成了n对，得到2n个标准正态样本，再截断取前n个。这个做法浪费了一半计算量，不过代码最直观。mu和sigma是线性变换参数，标准正态样本乘sigma再加mu，就得到均值为mu、标准差为sigma的正态分布。

参数说明里有两个容易搞混的地方。第一，sigma是标准差，不是方差。如果你手里只有方差var，必须写sigma = sqrt(var)。第二，rand每次调用都会推进全局随机数流，所以u1和u2之间是独立的。只要不手动重置随机种子，这个独立性就是可靠的。这里的max(rand(1, n), eps)是为了防止u1等于 0，虽然概率极低，但一旦发生就是 NaN。

3.2 向量化批量生成：一次性生成 m 对高斯样本的关键

循环写法容易理解，但大量样本时性能不够好。Matlab 的强项是向量化运算。常见做法是把n个样本按两列排列，一次生成n/2对，然后用.*和sqrt同时处理整列数据：

function z = box_muller_vec(n) m = ceil(n / 2); u = rand(m, 2); u(:, 1) = max(u(:, 1), eps); r = sqrt(-2 * log(u(:, 1))); theta = 2 * pi * u(:, 2); z = zeros(m, 2); z(:, 1) = r .* cos(theta); z(:, 2) = r .* sin(theta); z = reshape(z, 1, []); z = z(1:n); end

逻辑说明：u是m行 2 列的矩阵，第一列当U1，第二列当U2。r和theta都是列向量，cos(theta)的结果也是列向量，所以用.*做逐元素乘法。最后用reshape把多行掰成一行，再截掉多余样本。

参数说明里我要特别提醒运算符。sqrt(-2 * log(u(:, 1)))里的*是标量乘矩阵，正确；但r .* cos(theta)的.*是逐元素乘法，如果写成r * cos(theta)，Matlab 会把它当成矩阵乘法，维度对不上时直接报 "Matrix dimensions must agree"。这是新手踩得最多的一个点。另一个参数是m = ceil(n / 2)：当n是奇数时，需要多算一对然后截断，这是保持输出长度和输入一致的最简单方式。

3.3 把 uniformgauss 封装成双接口函数：method 参数怎么选

我更推荐的做法是把两种算法都留在一个函数里，用method参数切换。这样调试时可以直接对比两路输出，验证有没有写错：

function x = uniformgauss(n, mu, sigma, method) if nargin < 4 method = 'basic'; end if nargin < 3 sigma = 1; end if nargin < 2 mu = 0; end switch lower(method) case 'basic' m = ceil(n / 2); u1 = max(rand(m, 1), eps); u2 = rand(m, 1); z1 = sqrt(-2 * log(u1)) .* cos(2 * pi * u2); z2 = sqrt(-2 * log(u1)) .* sin(2 * pi * u2); z = reshape([z1, z2]', 1, []); z = z(1:n); case 'polar' z = zeros(1, n); i = 1; while i <= n v1 = 2 * rand - 1; v2 = 2 * rand - 1; s = v1^2 + v2^2; if s > 0 && s < 1 fac = sqrt(-2 * log(s) / s); z(i) = v1 * fac; if i + 1 <= n z(i+1) = v2 * fac; i = i + 1; end i = i + 1; end end otherwise error('method 只能是 basic 或 polar'); end x = z * sigma + mu; end

逻辑说明：basic分支和上小节公式一致，polar分支用 Marsaglia 方法。所有分支最终都产生标准正态样本z，然后再统一做z * sigma + mu的线性变换。这样就把"生成标准正态"和"调整均值方差"两件事解耦了。

参数说明：method的默认值是'basic'，因为我平时测试时优先保证结果可预期。如果你要生成海量数据且不关心代码美观度，'polar'的while循环虽然看着笨，但在老版本 Matlab 里反而稳定。sigma和mu支持传入向量吗？这里是单变量生成器，不支持。如果你要多元高斯分布，需要走第 6 章的 Cholesky 路线。

4. 怎么确认你生成的是真高斯：直方图、QQ 图与四个统计量一起上

4.1 直方图叠加理论 PDF：一眼看出形状

先跑一段最容易的验证：生成 100 万样本，画直方图，然后把理论概率密度函数曲线叠上去。Matlab 里用histogram的'Normalization', 'pdf'选项，让直方图面积归一化为 1，才能和概率密度函数放在同一坐标系：

x = uniformgauss(1000000, 2, 3, 'basic'); histogram(x, 200, 'Normalization', 'pdf'); hold on; t = linspace(min(x), max(x), 500); pdf_theory = exp(-0.5 * ((t - 2) / 3).^2) / (3 * sqrt(2 * pi)); plot(t, pdf_theory, 'r-', 'LineWidth', 2); hold off;

逻辑说明：linspace从样本最小到最大取 500 个点，理论 PDF 公式里除以3 * sqrt(2 * pi)是因为标准差是 3。如果直方图的包络和红线下偏或上偏，形状不对，说明前半段生成逻辑有问题。200这个 bin 数量在百万样本下足够看清尾部，再大会引入大量毛刺。

这个验证只能看整体形状，看不出尾部细节。高斯分布的尾部决定了极端值的行为，而直方图在高 sigma 区域几乎是平的。所以还要上 QQ 图。

4.2 QQ 图比直方图更能暴露尾部问题

QQ 图把样本分位数和理论分位数放在一起比。如果样本来自标准正态，点会落在y=x直线上。自定义均值方差时，想先用(x - mu) / sigma标准化：

x = uniformgauss(500000, 5, 2, 'polar'); z = (x - 5) / 2; qqplot(z);

逻辑说明：qqplot默认对比标准正态。如果两端点往下掉或往上翘，说明尾部比真高斯更薄或更厚。Box-Muller 最常见的尾部问题来自U1被错误扩到 0：log(0)会产生特别大的负半径，样本里出现一个远超合理范围的极值，QQ 图的右上角会突然飞出去。

qqplot的横轴是理论分位数，纵轴是样本分位数。点整体偏离y=x但保持线性关系，说明均值方差标定错了；如果是 S 形或弧线，说明分布形状不对。我一般先跑这条命令，再决定要不要深挖。需要注意，样本量小于几千时 QQ 图尾部抖动很厉害，建议至少 10 万样本起步。

4.3 四个统计量与假设检验：不能只看均值方差

有人验证到均值接近 0、方差接近 1 就觉得大功告成，这不够。高斯分布的偏度（skewness）应为 0，峰度（kurtosis）应为 3。偏度反映对称性，峰度反映尾部厚重程度。Matlab 直接算：

x = uniformgauss(1000000, 0, 1, 'basic'); fprintf('mean = %.4f\n', mean(x)); fprintf('std = %.4f\n', std(x)); fprintf('skew = %.4f\n', skewness(x)); fprintf('kurt = %.4f\n', kurtosis(x));

参数说明：样本量越大，这些统计量越靠近理论值。100 万样本下，均值通常偏差在 0.01 以内，标准差偏差在 0.005 以内，偏度在 0.02 以内，峰度在 0.05 以内。如果偏差超出这个量级，多半是生成逻辑有问题。下表是我在工程里的容忍参考值，基于 100 万样本：

如果还想做正式假设检验，可以用jbtest（Jarque-Bera 检验）或kstest。[h, p] = jbtest(x)的h=0表示不能拒绝正态假设。注意大样本下这类检验特别敏感，百万样本只要有一点点尾部偏差就会拒绝正态，所以别把h=1一票否决，要结合 QQ 图看偏差是否在工程可接受范围。

4.4 样本量多大才算数：从 1e5 到 1e6 的经验

样本量太小，什么检验都看不出来。我见过有人拿 200 个点做 QQ 图，然后断言"不是正态"，这是把抽样误差当系统性偏差。经验是：功能正确性用 1e5 样本验证，统计参数用 1e6 样本确认，极端性能调优再上 1e7。1e6 样本在现在电脑上跑 Box-Muller 只需零点几秒，完全不是负担。

如果你在写课程作业或复现论文，1e5 也够用。但要做蒙特卡洛积分或小概率事件模拟，尾部形状是关键，那就需要更多样本并同时看 QQ 图的上下两端。另一个技巧是固定随机种子后跑多次，观察统计量波动范围。如果每次均值都偏向同一侧，说明代码有结构性偏差；如果正负交替，那纯属随机波动。

5. Box-Muller 的常见坑与排查：从 NaN 到维度翻车

5.1 所有元素变成 NaN：log(0) 的边界

现象：生成的数组里出现大量 NaN，尤其在样本量特别大时。

原因：sqrt(-2 * log(U1))中U1等于 0 或极小值。log(0)是负无穷，sqrt(Inf)是无穷，再乘cos时如果cos恰好是 0，就是Inf * 0，结果是 NaN。即使U1不是 0 而是接近 0，比如 1e-300，log也会给出约 -690，开方后约 37，虽不算 NaN，但已经属于异常极值。

解决：在采样后立刻做U1 = max(U1, eps)。eps是 Matlab 的浮点相对精度，约为 2.22e-16，取它作为下限既能阻止U1=0，又不会显著改变分布形状。如果你在极坐标法里，对应的是判断s > 0而不是s >= 0。

5.2 方差始终偏小：忘了先乘标准差再加均值

现象：生成样本的均值正确，但方差只有理论值的四分之一甚至更小。

原因：最常见的错误是把标准正态样本直接当成目标分布输出，或者只取了 Box-Muller 生成的两个样本中的一个。设计算Z2 = sqrt(-2 ln U1) * sin(2π U2)被省略后，剩下的Z1确实是标准正态，但如果你用Z1的极坐标半径去套目标分布，就可能算漏sqrt(sigma)那层。还有一种情况是用户传进来的是方差，却直接拿来当标准差用。

解决：统一走这条链：标准正态Z，然后X = mu + sigma * Z。检查代码里sigma有没有开方。如果调用方写的是uniformgauss(n, 0, 4, 'basic')，那么它期望的标准差就是 4，方差是 16。直接var(x)核对，偏差不超过 1% 才合格。

5.3 Marsaglia 极坐标法死循环或输出错误：均匀范围搞错了

现象：程序卡在while循环里无法退出，或生成的样本明显偏斜。

原因：Marsaglia 极坐标法的输入是[-1,1)上的均匀分布，不是[0,1)。如果直接用rand，那么v1^2 + v2^2几乎全部落在(0,1)内，虽然能出结果，但产生的不是单位圆均匀采样，公式里sqrt(-2 ln(S) / S)的几何意义被破坏，输出分布不是高斯。另一个问题是s = v1^2 + v2^2有可能等于 0，这时候log(s)是负无穷，并且0/0构不成有效样本。

解决：v1 = 2 * rand - 1; v2 = 2 * rand - 1;，然后严格用if s > 0 && s < 1做接受条件。不要改成s <= 1，因为边界上极坐标变换的雅可比不成立。在调试时可加一个计数器统计拒绝率，正常情况下拒绝率接近1 - pi/4 ≈ 0.21，如果偏差太大，检查随机数范围。

5.4 报错 “Matrix dimensions must agree”：运算符用错

现象：向量化代码报维度不匹配，或结果矩阵形状变成n×n。

原因：把.*写成*。sqrt(-2 * log(U1))中log(U1)是 1×n 行向量，直接乘cos(2π U2)会触发矩阵乘法规则，维度对不上。另一种情况是用reshape时顺序搞错，导致样本被重排错位。

解决：逐元素运算一律用.*和./。另外，reshape(z1, 1, [])和reshape([z1; z2], 1, [])的结果顺序不一样，前者把z1排完再接z2，后者交替排列。如果只需要前n个样本，两种排法都行，但统计独立性取决于U1和U2是否独立，与排列顺序无关。出错时先把中间矩阵size打出来，看各维度是否符合预期。

5.5 两组样本几乎不独立：种子和复用问题

现象：连续调用两次uniformgauss，得到的两段样本相关系数接近 1，或者完全一样。

原因：在函数内部手动调了rng('default')或rng(seed)，每次调用都从同一个初始状态开始。另一个原因是把同一个U1向量复用到多个U2上，破坏了独立性。

解决：调试时可以用rng(42)固定种子，但正式生成数据前必须去掉，或每次调用rng('shuffle')。rand在 Matlab 中是全局流，只要不在每次调用前重置，天然独立。如果你把rng写死在函数里，那么函数之间的独立性就没了。检查方法是生成两段各 10 万样本，计算互相关系数，理想情况应接近 0，偏差在 0.01 以内才算正常。

6. Box-Muller 进阶：从标准高斯到多元高斯与一个复现验证技巧

6.1 用 Box-Muller 生成多元高斯：Cholesky 矩阵乘法

单变量场景只需要均值加标准差，但做卡尔曼滤波、高斯过程、投资组合模拟时，样本之间要带协方差。此时先生成标准正态矩阵，再做 Cholesky 分解即可：

mu = [1; 2]; Sigma = [2, 0.5; 0.5, 1]; L = chol(Sigma, 'lower'); n = 100000; Z = uniformgauss(n * 2, 0, 1, 'basic'); Z = reshape(Z, n, 2); X = mu' + Z * L';

逻辑说明：chol(Sigma, 'lower')得到下三角矩阵L，满足Sigma = L * L'。标准正态样本矩阵Z每一行是两个独立样本，乘L'后每一行的协方差结构就变成Sigma。mu'是 1×2 的行向量，加到每一行上完成平移。

参数说明：Sigma必须是半正定矩阵，否则chol报错。如果你发现Z * L'和L * Z'两种写法容易混，记住一点：样本按行排列时用Z * L'，样本按列排列时用L * Z。验证方式是对生成结果直接cov(X)，看是否接近Sigma。

6.2 固定种子回归测试：让随机性不背黑锅

我调试随机算法时有个习惯：先把随机种子固定住，把中间变量U1、U2也存下来，确认第一次输出符合手算预期。比如取U1 = 0.5, U2 = 0.25，手动算一遍sqrt(-2*ln(0.5))*cos(0.5π)，和程序输出比对。这能一票否决公式抄错、运算符写错等问题。在函数里加一行assert(all(isfinite(z)), 'Box-Muller 输入存在边界问题')，可以防止 NaN 一路传到后面的滤波链。

我的第一个 Box-Muller 程序是在粒子滤波里翻的车：状态转移噪声生成后，滤波跑了五分钟突然全部变成 NaN。排查下来就是log(0)惹的祸，从那以后每个随机生成器入口我都放一个isfinite断言，也从那时开始，所有随机算法先固定种子跑回归测试，再放开随机性跑批量实验。这套习惯帮我省下了大量和随机性纠缠的调试时间。希望帮到你。

本文还有配套的精品资源，点击获取