随机数生成和随机变量生成

dcjszhr

已于 2024-04-27 14:55:13 修改

阅读量970

点赞数 30

分类专栏：系统模拟文章标签：开发语言

于 2024-02-12 14:52:45 首次发布

本文链接：https://blog.csdn.net/dcjszhr/article/details/136070122

版权

系统模拟专栏收录该内容

10 篇文章 0 订阅

订阅专栏

目标

研究随机数生成（RNG）技术
生成（0,1）区间上均匀分布的随机数，U(0,1)
蒙特卡洛模拟
测试RNGs
生成具有其他分布的随机数（随机变量）

RNG：动机

随机数在模拟研究（随机模型）中经常被使用
随机数在许多现实生活情况中被使用，例如国家彩票、游戏和赌博
任何系统的模拟，其中包含固有的随机组件，需要一种生成随机数（变量）的方法
生成随机变量的先决条件是能够生成均匀分布的随机数

随机数生成技术

机械方法
- 投掷骰子
- 抽取编号球
- 发牌
数学方法（伪随机数）
- 平方中方法
- 线性同余生成器（LCGs）
- 乘法生成器
- 非线性同余生成器随机数的属性

两个重要的统计属性：

均匀性 - 如果区间（0, 1）被分成n个等分的子区间，每个区间内预期的观测数为N/n，其中N是观测总数
独立性（随机性）- 在特定区间观察到一个值的概率与之前的值无关

其他期望属性

不相关
可再生性
能够生成多个独立的流
长（最好是完整）周期
快速且低内存消耗

中间平方方法

$X_0$

$X_i = f$ ，其中f是 $(X_{i-1})^2$ 的中间4位数字，i = 1, 2, ..., n

$U_i = 0.X_i$

示例：X0 = 7182
它有效吗？
如果X0=0怎么办？

线性同余生成器（LCG）

X0：种子
对于i=1, 2, ..., n
Xi = (a * Xi-1 + c) mod m
Ui = Xi / m
其中m>0, a<m, c<m, X0 <m;
如果c = 0，则为乘法生成器。

真实的随机数生成器（示例）

参数
$m = 2^{31}-1 = 2,147,483,647$
$a = 314,159,269$
$c = 453,806,245$
随机数生成器的周期：p, 如果p = m，则周期完整
随机数生成器的种子

线性同余生成器: 例子

• m = 2^31 -1=2,147,483,647

• a = 2^ 7 =16,807

• c = 0

----------------------------

• X i = a. X i-1 mod m; U i = X i /m

• X 0 = 123,457 (RNG 种子)

• X 1 = 2,074,941,799; U 1 = X 1 /m=0.9662

• X 2 = 559,872,160; U 2 = X 2 /m=0.2607

线性同余生成器（续）

定理7.1（Law & Kelton）：线性同余生成器周期完整当且仅当
唯一能同时整除m和c的正整数是1
如果4整除m，则4也要整除a-1
如果q是一个质数且整除m，则q也要整除a-1
示例：m = 16, a = 5, c = 3

生成k个随机数流

将LCG的周期除以k
例如，如果P = 1,000,000，我们可以将其分成10个流
流1：X0 - X99,999
流2：X100,000 - X199,999
流3：X200,000 - X299,999
流10：X900,000 - X999,999

回顾：mm1.c程序 - rand.h

/* 以下三个声明用于使用随机数生成器rand, randst和randgt。任何使用这些函数的程序都应包含此文件（rand.h）*/

float randno(int stream);
void randst(long zset, int stream);
long randgt(int stream)

/* the following three declarations are for use with random numbers */
/* define the constants */
#define MODLUS 2147483647
#define MULT1 24112
#define MULT2 26143
/* set the default seeds for all 100 streams */
static long zrng[ ] =
{ 0,
1973272912, 281629770, 20006270, 1280689831, 2096730329, 1933576050,
913566091, 246780520, 1363774876, 604901985, 1511192140, 1259851944,
824064364, 150493284, 242708531, 75253171, 1964472944, 1202299975,
233217322, 1911216000, 726370533, 403498145, 993232223, 1103205531,
762430696, 1922803170, 1385516923, 76271663, 413692397, 726466604,
336157058, 1432650381, 1120463904, 595778810, 877722890, 1046574445,
68911991, 2088367019, 748545416, 622401386, 2122378830, 640690903,
1774806513, 2132545692, 2079249579, 78130110, 852776735, 1187867272, 
1351423507, 1645973084, 1997049139, 922510944, 2045512870, 898585771,
243649545, 1004818771, 773686062, 403188473, 372279877, 1901633463,
498067494, 2087759558, 493157915, 597104727, 1530940798, 1814496276,
536444882, 1663153658, 855503735, 67784357, 1432404475, 619691088,
119025595, 880802310, 176192644, 1116780070, 277854671, 1366580350,
1142483975, 2026948561, 1053920743, 786262391, 1792203830, 1494667770, 
1923011392, 1433700034, 1244184613, 1147297105, 539712780, 1545929719,
190641742, 1645390429, 264907697, 620389253, 1502074852, 927711160,
364849192, 2049576050, 638580085, 547070247
};

蒙特卡洛模拟

对于难以数学解决的概率问题非常有用
基于频率理论的概率 : 如果一个随机实验重复n次，且x是事件X发生的次数，则X发生的相对频率是x/n。频率理论的概率声明，这个相对频率随着n趋向于∞而收敛于X的实际概率

要计算掷两个骰子得到点数和为7的概率，我们可以直接根据您提供的信息来进行计算。

总样本空间中的情况数为36（每个骰子有6个面，因此有6×6=366×6=36种可能的结果）。
得到点数和为7的情况有6种，即：(1,6), (2,5), (3,4), (4,3), (5,2), (6,1)。

因此，得到7的概率P(7)为：P(7)=6/36

让我们计算这个概率的确切值。

掷两个骰子得到点数和为7的概率是1/6或大约16.67%。这意味着如果我们模拟掷骰子6000次，我们预期大约会有1000次出现点数和为7。

示例：通过蒙特卡洛模拟计算π

方法：在一个正方形内形成一个圆。圆的半径为单位长度，所以正方形的边长为2×2。
E1：对于i=1到n执行
在(0, 1)内生成两个随机数xi和yi
用m表示落在圆内的点(xi, yi)的数量。
E2：对于非常大的n，m/n接近于π/4

测试随机数生成器

需要测试均匀性和独立性

独立性测试

独立性意味着生成的随机数之间没有相关性。独立性的测试通常涉及分析随机数序列，以寻找任何可预测的模式或重复出现的序列。

均匀性测试

均匀性意味着生成的随机数在其定义域内等概率地出现。测试这一特性的假设（H0）是随机数在(0,1)区间上均匀分布。
与已知结果比较：可以通过与简单的随机实验比较来检验随机数的均匀性，例如：
- 抛硬币：理想情况下，硬币的正反面出现概率应该相等，可以用来比较随机数生成器生成0和1的分布是否均匀。
- 掷骰子：掷骰子的每一个面出现的概率应该相同，用来检验随机数生成器在一个有限集合内生成数字的均匀性。
卡方检验（Chi-square test）：这是一种统计检验方法，用于比较观察频率和期望频率之间的差异，以判断随机数的分布是否与期望的均匀分布相符。
Kolmogorov-Smirnov（K-S）检验：这是一种非参数检验方法，用于比较两个样本分布或一个样本分布与参考分布之间的差异。K-S检验特别适用于检验随机数的连续分布的均匀性。
序列相关性分析：检查随机数序列中的数值是否存在相关性，比如通过计算序列的自相关函数。
游程测试（Runs test）：这是一种用于检测随机数序列中随机性的非参数检验，通过分析连续数字出现的游程长度和分布来评估序列的随机性。

均匀性测试：卡方检验

给定Ui，i = 1, …, n
H0：Ui' 在[0,1]上均匀分布
将[0,1]分成k个子区间I1, …, Ik
让Oj：Ij中Ui的数量，j=1,2,…,k
$\chi ^2=\sum_{j=1}^n \frac{(O_j-E_j)^2}{E_j}$
如果 $\chi 2 > \chi ^2_\alpha ,k-1$ ，则拒绝H0。
其中临界值 $\chi ^2_\alpha$ ,k-1在Banks的表A6（或统计表）中查找

示例：卡方检验

0.34 0.90 0.25 0.89 0.17 0.44 0.12 0.21 0.46 0.67

0.83 0.76 0.79 0.64 0.70 0.11 0.44 0.74 0.22 0.74

0.06 0.99 0.77 0.67 0.56 0.41 0.52 0.73 0.98 0.12

这些数字由RNG生成。这些数字在[0, 1]中均匀分布吗？

根据表A6，χ^2_0.05,4 = 9.49，因此H0不被拒绝

均匀性测试：Kolmogorov-Smirnov（K-S）检验

H0：Ri在[0, 1]上均匀分布
对于给定的显著性水平α（例如，α=0.05），如果D > Dα，则拒绝H0
其中D+ = Max{i/N - Ri}；
D-=Max{Ri - (i-1)/N}；
D = max(D+, D-)
Dα的值可以从教科书的表A8中读取
H0在D > Dα时被拒绝
例如，0.08, 0.11, 0.55, 0.03
D+ = 0.64; D- =0.03; D=0.64;
根据表A8，Dα=0.624
结论：H0被拒绝

注：Ri 必须按升序排序

独立性测试：游程测试

测试序列中的“上升和下降的游程”或“高于和低于平均值的游程”，通过比较实际值与预期值
检查序列中数字的排列，以测试独立性假设

游程上升和下降

游程定义为一个事件的连续发生，前后由不同的事件界定。
游程的长度是发生在游程中的事件数量
例如，抛硬币10次 H ｜ T T ｜ H H ｜ T T T ｜ H ｜ T，游程数 = 6

例如： 0.87 0.15 0.23 0.45 0.69 0.32 0.30 0.19 0.24 0.18 0.65 0.82 0.93 0.22 0.81
数字根据后续数字是更大还是更小分别标记为“+”或“-”。序列为：- + + + - - - + - + + + - +
共有8个游程，四个上升游程和四个下降游程。

如果序列中的数字总数为N，则游程的最大数量为N-1，最小数量为1。理想情况下，游程的数量会在这两个极端之间。
设a为序列中游程的总数，可以证明，如果序列中的数字是随机选取的，则
$E(a) = \mu_a = (2N-1)/3$
$V(a) = \sigma_a^2 = (16N-29)/90$
对于N > 20，游程的分布可以近似为正态分布 $N(\mu_a, \sigma_a^2)$

使用标准正态分布
Z0 = (a - μa) / σa,
其中Z0是N(0,1)，a是观察到的游程数量
如果-zα/2 ≤ Z0 ≤ zα/2，则不拒绝独立性假设（Ho），其中α是显著性水平

例子

给定以下40个数字序列，当α = 0.05时，可以拒绝独立性假设吗？

0.41 0.68 0.89 0.94 0.74 0.91 0.55 0.62 0.36 0.27

0.19 0.72 0.75 0.08 0.54 0.02 0.01 0.36 0.16 0.28

0.18 0.01 0.95 0.69 0.18 0.47 0.23 0.32 0.82 0.53

0.31 0.42 0.73 0.04 0.83 0.45 0.13 0.57 0.63 0.29

产生其他分布的随机数：逆变换法

假设变量X的累积分布函数（cdf）为F(X)，
F(x) = Prob(X≤x)；1 ≥ F(x) ≥ 0
F(x)是连续且非递减的，即x2 > x1 ⇒ F(x2) ≥ F(x1)
目标：生成具有F(X)为cdf的Xi

在[0,1]中生成Ui
设Ui = F(Xi)
计算Xi = F^-1(Ui)
如果需要更多数字，则转到1

生成指数分布的随机数

示例：服务器的服务时间是指数分布的，其服务率=λ（或其平均服务时间=1/λ）

生成指数分布的随机变量

生成指数分布随机变量的函数

float expon(float mean) /* 指数变量生成函数 */
{
    float u;
    /* 生成一个U(0,1)随机变量 */
    u = rand(seed); /* 例如，seed=2000119 */
    /* 返回均值为"mean"的指数随机变量 */
    return -mean*log(u);
}

示例：X服从均值为5的指数分布
调用expon(5)将提供序列的下一个元素

逆变换法生成离散分布

$p(x_i)=p(X=x_i)\\ F(x_i)=p(X\leq x_i) =\sum_{k=1}^ip(x_k)$

从U(0, 1)生成X
找到最小的整数i，使得F(xi) ≥ X;
返回xi
示例：等概率生成硬币的正面和反面。用x1表示正面，x2表示反面；p(x1)= p(x2)=0.5
- 从U(0, 1)生成X
- 假设X=0.76 => i=2;
- 返回反面

离散分布

示例生成x1, x2, …, x6，其中 p(x1)= p(x2)= p(x3)= p(x4)= p(x5)= p(x6)= 1/6=0.166
问题：生成具有上述分布的X
解决方案
- 从U(0, 1)生成X
- 假设X=0.45
- 找到最小的整数i，使得F(xi) ≥ X
- X=0.45==> i = 3
- 返回x3

M/M/1模拟

M/M/1是一个单服务站队列，假设
- 到达时间是IID且服从指数分布
- 服务时间是IID且服从指数分布
时间测量指标
- 客户在队列中的平均等待时间
- 队列中客户的平均时间加权数量
- 服务器的平均利用率
如何生成到达事件？

M/M/1中生成下一个到达

void arrive(void) /* 到达事件函数 */
{
    float delay;
    /* 安排下一个到达 */
    time_next_event = current_time + expon(mean_interarrival);
    /* 检查服务器是否繁忙 */
    if (server_status == BUSY) 
    {
        /* 服务器繁忙，因此增加队列中的客户数量 */
        ++num_in_q;
        /* 检查是否存在溢出条件 */
        if (num_in_q > Q_LIMIT)
        { 
            /* 队列已溢出，因此停止模拟 */
            printf("\nOverflow of the array time_arrival at time %f\n", current_time);
            exit(2);
        }
        /* 队列中还有空间，因此在time_arrival的（新）末端存储到达客户的到达时间 */
        time_arrival[num_in_q] = current_time;
    }
    else

练习：生成具有三角分布的随机变量

错误简化模拟模型

考虑一个单服务系统，λ = 9且µ = 10

使用平均值计算平均队列长度会怎样？
使用M/M/1公式计算平均队列长度会怎样？

M/M/1理论指标

λ = 到达率（每个时间单位到达的平均客户数）
µ = 服务率（每个时间单位服务的平均客户数）
ρ = 服务器利用率 = λ /µ
w = 客户在队列中的平均等待时间
q = 队列中客户的（时间加权）平均数量

$\rho = \frac{\lambda}{\mu}\\ W=\frac{\rho}{\mu-\lambda}=\frac{\lambda}{\mu(\mu-\lambda)}\\ q=\frac{\lambda^2}{\mu(\mu-\lambda)}$

如何使用理论指标？

示例1：λ = 4; µ = 5 ρ =0.8; w=0.8; q=3.2;
示例2：λ = 8; µ = 10 ρ =0.8; w=0.4; q=3.2;
示例：人们倾向于简化模拟模型以简化计算。考虑一个单服务系统，λ = 9且µ = 10
- 使用平均值，队列将不存在 ==> q = 0
- 使用M/M/1公式得出： $q=\frac{\lambda^2}{\mu(\mu-\lambda)}=\frac{81}{10}=8.1$