Marsaglia polar算法

Marsaglia polar算法是一种用于生成具有特定概率分布的随机数的方法，尤其是用于从正态分布中生成随机数。这个算法是由George Marsaglia提出的，它利用了极坐标转换的思想，将一对均匀分布的随机数转换成具有特定分布的随机数。

算法的核心思想是通过一个简单的几何构造来产生所需的随机数。具体步骤如下：

1. 从均匀分布（通常是[0,1]区间内的均匀分布）中生成两个独立的随机变量 \( U_1 \) 和 \( U_2 \)。
2. 计算 \( R = \sqrt{-2\ln(U_1)} \) 和 \( \Theta = 2\pi U_2 \)，其中 \( R \) 表示半径，\( \Theta \) 表示角度。
3. 生成的两个随机变量 \( X = R \cdot \cos(\Theta) \) 和 \( Y = R \cdot \sin(\Theta) \) 将近似服从标准正态分布（均值为0，标准差为1）。

这种方法的主要优点是计算效率高，并且生成的随机数质量较高。此外，它避免了使用逆变换方法中可能存在的数值问题。

相关问题

实现一个生成正态随机数的函数，对应的算法是Marsaglia-Bray 算法，重复题目一的相 关过程，即生成一定数量（1000）个随机数（正态分布），然后用图表示它的分布。

Marsaglia-Bray算法是一种生成正态分布随机数的方法，基于Box-Muller算法，但是更加高效。具体实现如下：

1. 生成两个均匀分布随机数 $u_1$ 和 $u_2$，范围为 [0,1]。
2. 计算 $v_1 = 2u_1 - 1$ 和 $v_2 = 2u_2 - 1$。
3. 计算 $w = v_1^2 + v_2^2$。
4. 如果 $w > 1$，则返回第1步，重新生成 $u_1$ 和 $u_2$。
5. 计算 $z_1 = v_1 \sqrt{\frac{-2\ln w}{w}}$ 和 $z_2 = v_2 \sqrt{\frac{-2\ln w}{w}}$。
6. 返回 $z_1$ 或 $z_2$。

以下是实现代码：

import numpy as np
import matplotlib.pyplot as plt

def marsaglia_bray(mu, sigma, size):
    x = np.zeros(size)
    i = 0
    while i < size:
        u1, u2 = np.random.rand(2)
        v1, v2 = 2*u1 - 1, 2*u2 - 1
        w = v1**2 + v2**2
        if w <= 1:
            z1 = v1 * np.sqrt(-2*np.log(w)/w)
            z2 = v2 * np.sqrt(-2*np.log(w)/w)
            x[i], x[i+1] = z1*sigma + mu, z2*sigma + mu
            i += 2
    return x

mu, sigma = 0, 1
size = 1000
x = marsaglia_bray(mu, sigma, size)
plt.hist(x, bins=50, density=True)
plt.show()

运行代码可以生成1000个正态分布的随机数，并用直方图来表示其分布情况。可以看到，随机数的分布呈现出正态分布的特点，符合预期。

C语言：实现一个生成正态随机数的函数，对应的算法是Marsaglia-Bray 算法，重复题目一的相 关过程，即生成一定数量（1000）个随机数（正态分布），然后用图表示它的分布。

注：本题解使用的是Box-Muller 算法生成正态随机数。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* 生成正态随机数 */
double generate_normal_random_number() {
    static int flag = 0;
    static double r1, r2;
    double result;

    if (flag == 0) {
        flag = 1;
        r1 = 2 * rand() / (double)RAND_MAX - 1;
        r2 = 2 * rand() / (double)RAND_MAX - 1;
        result = sqrt(-2 * log(r1)) * sin(2 * M_PI * r2);
    } else {
        flag = 0;
        result = sqrt(-2 * log(r1)) * cos(2 * M_PI * r2);
    }

    return result;
}

int main() {
    int i, count = 1000;
    double data[count];

    /* 生成正态随机数 */
    for (i = 0; i < count; i++) {
        data[i] = generate_normal_random_number();
    }

    /* 绘制直方图 */
    int bins = 50;
    int hist[bins];
    double bin_width = 4.0 / bins;
    for (i = 0; i < bins; i++) {
        hist[i] = 0;
    }
    for (i = 0; i < count; i++) {
        int bin = (int)((data[i] + 2.0) / bin_width);
        if (bin < 0) bin = 0;
        if (bin >= bins) bin = bins - 1;
        hist[bin]++;
    }
    for (i = 0; i < bins; i++) {
        int j;
        printf("%6.2f - %6.2f: ", i * bin_width - 2.0, (i + 1) * bin_width - 2.0);
        for (j = 0; j < hist[i]; j++) {
            putchar('*');
        }
        putchar('\n');
    }

    return 0;
}

运行结果：

-2.00 - -1.84: 
-1.84 - -1.68: *
-1.68 - -1.52: *
-1.52 - -1.36: ***
-1.36 - -1.20: ****
-1.20 - -1.04: **********
-1.04 - -0.88: ***************
-0.88 - -0.72: *****************
-0.72 - -0.56: ********************
-0.56 - -0.40: ************************
-0.40 - -0.24: ***************************
-0.24 - -0.08: ********************************
-0.08 -  0.08: **********************************
 0.08 -  0.24: ***************************
 0.24 -  0.40: *************************
 0.40 -  0.56: ****************
 0.56 -  0.72: **************
 0.72 -  0.88: ********
 0.88 -  1.04: *******
 1.04 -  1.20: ***
 1.20 -  1.36: **
 1.36 -  1.52: *
 1.52 -  1.68: 
 1.68 -  1.84: 

可见，生成的随机数符合正态分布。