利用Box-Muller变换生成正态分布的随机数（附代码）

       前言

       之前的文章中笔者有提到利用独立同分布的中心极限定理生成正态分布的随机数，对随机数的分布特性要求越高，利用独立同分布的中心极限定理生成正态分布的随机数的算法时间复杂度也越高，如果需要大量的正态分布随机数“利用独立同分布的中心极限定理生成正态分布的随机数”的做法是不可取的。本文将为大家带来时间复杂度较低的算法——利用Box-Muller变换生成正态分布的随机数。If you want the proof process, please see the end of the article.

      Box-Muller变换

      设  是在[0,1)上遵从均匀分布的随机数（生成[0,1)上遵从均匀分布的随机数通常使用的是梅森旋转算法），令

 

则 必然遵守二元正态分布，即利用两个独立的遵从均匀分布的随机数得到两个独立的正态分布的随机数

定理证明过程

      设

     令

       

由概率密度函数的变换公式可得

（雅可比式的定义可以参见数学分析教材）

因为 在[0,1)上遵从均匀分布,故

通过计算雅可比式的值，可得

联立（1）（2）可得

故 ​​​​​​​​​​​​​​

此式恰为二元标准正态分布的概率密度函数

证毕

C++实现过程

//作者cclplus 如有疑问可请联系作者邮箱maxwell970710@gmail.com，作者会在后续内容中进行解释
#include <iostream>
#include <cmath>
#include <cstdlib>

using namespace std;
const double pi = 3.1415926897932384;
int main() {
	ios::sync_with_stdio(false);
	double x1, x2,y1,y2;
	int seed;//手动输入种子
	cout << "手动输入种子(任意给出一个整数)" << endl;
	cin >> seed;
	srand(seed);
	x1 = rand()% RAND_MAX /(double) RAND_MAX;
	x2= rand() % RAND_MAX / (double)RAND_MAX;
	cout << "产生的均匀分布的随机数种子为" << endl;
	cout << x1 <<" "<< x2 << endl;
	y1 = sqrt(-2 * log(x1))*cos(2 * pi*x2);
	y2 = sqrt(-2 * log(x2))*sin(2 * pi*x1);
	cout << "输出的正态分布的随机数为" << endl;
	cout << y1 << " " << y2 << endl;
	return 0;
}

利用matlab设计验证性实验

生成一百万个正态分布的函数，并对结果进行验证

%作者cclplus 有疑问请联系maxwell970710@gmail.com
clc
clear all
close all
n=1000000;
x1=rand(n,1);
x2=rand(n,1);
y1=sqrt(-2.*log(1.*x1)).*cos(2*pi.*x2);%生成正态分布函数
figure;
normplot(y1);%样本数据在图中用“+”显示；如果数据来自正态分布，则图形显示为直线，而其它分布可能在图中产生弯曲。

 

测试结果

如果想要生成其它类型的正态分布随机数，可以参考我之前写的文章，利用独立同分布的中心极限定理生成正态分布的随机数

Proof process 

It is assumed that the uniformly distributed random numbers  are obeyed on [0, 1) (the generation of [0, 1) random numbers that follow a uniform distribution is usually the Mason rotation algorithm)

assume

then，   is inevitable to follow the binary normal distribution, that is, to obtain two independent normal distribution random numbers by using two independent random numbers that follow the uniform distribution.

assume

The transformation formula of the probability density function can be used to obtain the following formula

  Uniform distribution on [0,1)

By calculating the value of the Jacobian formula, you can get

Calculate the above formula, you can get

The formula exactly the probability density function of the binary standard normal distribution.

Q.E.D.

 

拓展阅读   生成0-1分布的均匀变量

```*

#include <iostream>
using namespace std;
double rand_1(double * r) {
    int m;
    double s, u, v, ret;
    s = 65536.0;
    u = 2053.0;
    v = 13849.0;
    m = (int)(*r / s);
    *r = *r - m * s;
    *r = u * *r + v;
    m = (int)(*r / s);
    *r = *r - m * s;
    ret =*r / s;
    return ret;
}

int main()
{
    double r = 5;//种子（只有改变种子，所得随机数组才能发生变化）
    for (int i = 0; i < 10; i++) {
        cout << rand_1(&r) << endl;
    }
}
 

```