Skip to content
Lex Blog

Box-Muller 变换

Box-Muller 变换是一种利用服从均匀分布的随机变量构造服从高斯分布的随机变量的方法。

内容

Box-Muller 变换是一种从均匀分布的随机变量生成正态分布随机变量的方法。

具体来说，选取两个在区间 (0,1) 上均匀分布的随机变量 U1U_1 与 U2U_2。按如下方式构造随机变量 XX 与 YY：

X=cos⁡(2πU1)−2ln⁡U2X = \cos(2\pi U_1) \sqrt{-2\ln U_2}Y=sin⁡(2πU1)−2ln⁡U2Y = \sin(2\pi U_1) \sqrt{-2\ln U_2}

此时，XX 与 YY 服从均值为 0、方差为 1 的正态分布。

推导

设 XX 与 YY 是相互独立、均值为 0、方差为 1 的正态分布随机变量。记 p(X)p(X) 与 p(Y)p(Y) 分别为其概率密度函数。于是有：

p(X)=12πe−X22,p(Y)=12πe−Y22p(X) = \frac{1}{\sqrt{2\pi}} e^{-\frac{X^2}{2}}, \\ p(Y) = \frac{1}{\sqrt{2\pi}} e^{-\frac{Y^2}{2}}

由于 XX 与 YY 相互独立，它们的联合概率密度函数为：

p(X,Y)=12πe−X2+Y22p(X,Y) = \frac{1}{2\pi} e^{-\frac{X^2 + Y^2}{2}}

现在进行坐标变换，定义：

X=Rcos⁡(θ),Y=Rsin⁡(θ)X = R \cos(\theta), \\ Y = R \sin(\theta)

于是，在极坐标下，联合分布变为：

∫0∞∫02π12πe−R22R dθ dR=1\int_{0}^{\infty} \int_{0}^{2\pi} \frac{1}{2\pi} e^{-\frac{R^2}{2}} R \, d\theta \, dR = 1

由此，可推导出 RR 与 θ\theta 的分布函数，分别记为 PRP_R 与 PθP_\theta：

PR(R≤r)=∫0r∫02π12πe−R22R dθ dR=1−e−r22P_R(R \leq r) = \\ \int_{0}^{r} \int_{0}^{2\pi} \frac{1}{2\pi} e^{-\frac{R^2}{2}} R \, d\theta \, dR = 1 - e^{-\frac{r^2}{2}}Pθ(θ≤ϕ)=∫0ϕ∫0∞12πe−R22R dR dθ=ϕ2πP_\theta(\theta \leq \phi) = \\ \int_{0}^{\phi} \int_{0}^{\infty} \frac{1}{2\pi} e^{-\frac{R^2}{2}} R \, dR \, d\theta = \frac{\phi}{2\pi}

显然，θ\theta 在 [0,2π][0, 2\pi] 上服从均匀分布。设

FR=1−e−r22F_R = 1 - e^{-\frac{r^2}{2}}

其反函数为：

R=FR−1(z)=−2ln⁡(1−z)R = F_R^{-1}(z) \\ = \sqrt{-2 \ln(1 - z)}

当 zz 在 [0,1][0,1] 上服从均匀分布时，RR 的分布函数为 FR(r)F_R(r)。因此，我们可以选取两个在 [0,1][0,1] 上均匀分布的随机变量 U1U_1 与 U2U_2，使得：

θ=2πU1,1−z=U2,R=−2ln⁡U2\theta = 2\pi U_1, \\ 1 - z = U_2, \\ R = \sqrt{-2 \ln U_2}

将上述各式代入前面的表达式中，可得：

X=Rcos⁡(θ),Y=Rsin⁡(θ)X = R \cos(\theta), \\ Y = R \sin(\theta)

由此恢复了 XX 与 YY 的原始表达式，它们服从均值为 0、方差为 1 的正态分布。

验证

为了生成服从标准正态分布的随机变量 X∼N(0,1)X \sim N(0,1)，我们注意到误差函数 erf−1\mathrm{erf}^{-1} 的反函数不是初等函数。用级数展开来近似这个反函数会引入较大的计算误差。

另一方面，初等函数具有成熟的高精度算法。因此，我们希望用初等函数从均匀分布随机变量构造正态分布随机变量。

设 U,V∈(0,1)U, V \in (0,1) 为均匀分布的随机变量。考虑如下变换：

X=−2ln⁡Ucos⁡(2πV),Y=−2ln⁡Usin⁡(2πV)X = \sqrt{-2\ln U} \cos(2\pi V), \\ Y = \sqrt{-2\ln U} \sin(2\pi V)

接下来，我们验证 XX 与 YY 相互独立且服从标准正态分布。

定义：

R=X2+Y2,Θ=arctan⁡(YX)R = \sqrt{X^2 + Y^2}, \\ \Theta = \arctan\left(\frac{Y}{X}\right)

注意到 Θ=2πV\Theta = 2\pi V，R=−2ln⁡UR = \sqrt{-2\ln U}。下面证明在极坐标变换下，RR 与 Θ\Theta 相互独立，并且它们产生标准正态随机变量。

考虑其中一个变量固定时另一个变量的概率密度函数。Θ\Theta 的概率密度函数为：

PΘ(Θ)=12πP_\Theta(\Theta) = \frac{1}{2\pi}

这很直接。对于 RR，其累积分布函数为：

P(R≤R0)=P(U≤exp⁡(−R022))P(R \leq R_0) = \\ P\left(U \leq \exp\left(\frac{-R_0^2}{2}\right)\right)

因此，RR 的概率密度函数为：

PR(R)=Re−R22P_R(R) = R e^{-\frac{R^2}{2}}

因此，RR 与 Θ\Theta 的联合概率密度为：

P(R,Θ)=12πRe−R22P(R, \Theta) = \frac{1}{2\pi} R e^{-\frac{R^2}{2}}

由于极坐标变换给出：

R dR dΘ=dX dYR \, dR \, d\Theta = dX \, dY

我们得到 XX 与 YY 的联合概率密度：

P(X,Y)=12πe−X2+Y22=(12πe−X22)(12πe−Y22)P(X,Y) = \frac{1}{2\pi} e^{-\frac{X^2 + Y^2}{2}} \\ = \left(\frac{1}{\sqrt{2\pi}} e^{-\frac{X^2}{2}}\right) \left(\frac{1}{\sqrt{2\pi}} e^{-\frac{Y^2}{2}}\right)

这正是两个相互独立、服从标准正态分布的随机变量 XX 与 YY 的概率密度函数。

代码实现

function randInterval(): number {
  return Math.random();
}

function boxMuller(mu: number, sigma: number): [number, number] {
  const u = randInterval();
  const v = randInterval();
  const x = Math.cos(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  const y = Math.sin(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  return [x * sigma + mu, y * sigma + mu];
}

// 使用示例
const [value1, value2] = boxMuller(0, 1);
console.log(value1, value2);
// 输出两个服从正态分布的随机数

为了验证 boxMuller 函数生成的数服从正态分布，可以生成大量的随机数并绘制直方图与标准正态分布进行比较。你可以使用 plotly 等库将直方图可视化并评估输出。

下面是一个测试程序，它使用 boxMuller 函数生成数据集，进行基本的统计检验，并绘制直方图。假设你在 Node.js 环境中运行，可以使用 plotly-nodejs，或直接在浏览器中使用 plotly.js 生成图形。

首先，安装所需的依赖：

npm i plotly.js-dist

然后，执行以下代码：

import * as fs from "fs";
import * as plotly from "plotly.js-dist";

// 使用 Box-Muller 生成服从正态分布的随机数
function randInterval(): number {
  return Math.random();
}

function boxMuller(mu: number, sigma: number): [number, number] {
  const u = randInterval();
  const v = randInterval();
  const x = Math.cos(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  const y = Math.sin(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  return [x * sigma + mu, y * sigma + mu];
}

// 生成 n 个服从正态分布的随机数
function generateNormalDistribution(
  mu: number,
  sigma: number,
  n: number
): number[] {
  const values: number[] = [];
  for (let i = 0; i < n / 2; i++) {
    const [value1, value2] = boxMuller(mu, sigma);
    values.push(value1, value2);
  }
  return values;
}

// 生成数据
const mu = 0;
const sigma = 1;
const sampleSize = 10000; // 样本大小
const values = generateNormalDistribution(mu, sigma, sampleSize);

// 计算生成的随机数的均值与方差
const mean = values.reduce((acc, val) => acc + val, 0) / values.length;
const variance =
  values.reduce((acc, val) => acc + (val - mean) ** 2, 0) / values.length;

console.log(`均值：${mean}`);
console.log(`方差：${variance}`);

// 绘制直方图
const trace = {
  x: values,
  type: "histogram",
  xbins: {
    size: 0.1, // 直方图的组距
  },
  marker: {
    color: "blue",
  },
  opacity: 0.7,
  name: "Box-Muller 分布",
};

const data = [trace];

const layout = {
  title: "生成的正态分布",
  xaxis: { title: "数值" },
  yaxis: { title: "频数" },
  bargap: 0.05,
};

const graphOptions = {
  filename: "box-muller-distribution",
  fileopt: "overwrite",
};
fs.writeFileSync("box-muller.html", plotly.plot(data, layout, graphOptions));

console.log("Box-Muller 分布数据已保存为 box-muller.html。");

该程序使用 boxMuller 函数生成 10,000 个服从正态分布的随机数，并计算所生成数据集的均值与方差。随机数通过直方图可视化，便于与标准正态分布直接比较。直方图保存为名为 box-muller.html 的文件，可在浏览器中打开查看。