Box-Muller法生成正态分布

Pγ

已于 2024-09-12 16:21:09 修改

阅读量233

点赞数 8

文章标签：算法

于 2024-09-05 17:14:19 首次发布

本文链接：https://blog.csdn.net/weixin_45121008/article/details/141932278

版权

理论推导过程

二维标准正态分布的变量相互独立，可以写为：

$f(x_1,x_2)=\frac{1}{2\pi}e^{-\frac{x_1^2+x_2^2}{2}}=f(x_1)f(x_2)$

令:

$r^2=x_1^2+x_2^2$

此时可以把笛卡尔坐标系的(x1,x2)转化为极坐标系的(r,θ)：

$\left\{\begin{matrix}x_1=rcos\theta \\ x_2=rsin\theta \end{matrix}\right.$

$\left\{\begin{matrix}r=\sqrt{x_1^2+x_2^2} \\ \theta=arctan\frac{x_2}{x_1} \end{matrix}\right.$

改写分布为：

$f(r,\theta)=\frac{1}{2 \pi}e^{-\frac{r^2}{2}}$

需要注意的是在极坐标下求概率分布（求面积）是：

$\int_{}^{}f(r,\theta)rdrd\theta$

此时角度θ服从[0,2π)的均匀分布，r2服从k=2的卡方分布，也就是λ=1的指数分布。作者提出如果一个分布U服从[0,1)的均匀分布，那么-2lnU服从卡方分布。证明如GPT所示：

同时需要注意到z的定义域为(0,+∞)，也是符合卡方分布或者指数分布的。

那么可以利用 $r^2=x_1^2+x_2^2$ 服从k=2的卡方分布，反推得到一个U1服从均匀分布，刚好这个U1的形式是：

$U_1=e^{-\frac{x_1^2+x_2^2}{2}}$

同时前面提到θ服从[0,2π)的均匀分布，那么 $\frac{1}{2\pi}arctan\frac{x_2}{x_1}$ 服从[0,1)的均匀分布。现在问题就从正态分布转换为了两个均匀分布：

$\left\{\begin{matrix}U_1=e^{-\frac{x_1^2+x_2^2}{2}} \\ U_2=\frac{1}{2\pi}arctan\frac{x_2}{x_1} \end{matrix}\right.$

反变换为：

$\left\{\begin{matrix}x_1=\sqrt{-2ln(U_1)}cos(2\pi U_2) \\ x_2=\sqrt{-2ln(U_1)}sin(2\pi U_2) \end{matrix}\right.$

对计算机来说用table来实现均匀分布要好实现多了。而标准二维分布又可以拆成两个一维的标准正态分布相乘，那么一维的分布就可以得到了。

Matlab代码

在Matlab中rand函数是生成(0, 1)之间均匀分布的随机数，因此生成u1和u2，再相乘，可以得到一个一维高斯分布的平方。分布是怎么样的可以利用直方图函数hist观察。

% 参数
N = 1000;  % 样本数
u1 = rand(1, N);  % 生成均匀分布u1
u2 = rand(1, N);  % 生成均匀分布u2

% Box-Muller变换
z1 = sqrt(-2 * log(u1)) .* cos(2 * pi * u2);  % 第一维高斯分布
z2 = sqrt(-2 * log(u1)) .* sin(2 * pi * u2);  % 第二维高斯分布

% 提取其中一个维度的数据作为一维正态分布
z = z1;  % 选择第一维

% 绘制一维高斯分布的直方图
figure;
histogram(z, 30, 'Normalization', 'pdf');  % 绘制正态分布的直方图
hold on;

% 绘制理论正态分布曲线
x = linspace(min(z), max(z), 100);
y = normpdf(x, 0, 1);  % 标准正态分布的概率密度函数
plot(x, y, 'r', 'LineWidth', 2);  % 绘制理论正态分布曲线
title('One-dimensional Gaussian Distribution');
xlabel('Value');
ylabel('Probability Density');
legend('Sampled Data', 'Theoretical Gaussian');
hold off;