文章目录
- 1.伪随机数介绍
- 1.1.伪随机产生的意义
- 1.2.伪随机产生的过程
- 2.产生U(0,1)的乘除同余法
- 2.1.原始的乘同余法
- 2.2.改进的乘同余法
- 3.产生正态分布的伪随机数
- 4.基于逆变法产生伪随机数
1.伪随机数介绍
1.1.伪随机产生的意义
1.随机数的产生是进行随机优化的第一步也是最重要的一步,随机优化算法运行过程中需要大量随机数。
2.传统手工方法:抽签,掷骰子,抽牌,摇号等,无法满足产生大量随机数的需求。
3.伪随机数方法:利用计算机通过某些数学公式计算而产生,从数学意义上说不是随机的,但只要通过随机数的一系列统计检验,就可以作为随机数来使用。
1.2.伪随机产生的过程
∙ \bullet ∙Step1:确定一个数学模型或某种规则。
∙ \bullet ∙Step2:规定几个初始值。
∙ \bullet ∙Step3:按照上述模型产生第一个随机数。
∙ \bullet ∙Step4:用产生的上一个随机数作为新的初值,按照相同的步骤产生下一个随机数,重复之,得到一个伪随机数序列。
2.产生U(0,1)的乘除同余法
2.1.原始的乘同余法
均匀的随机数的产生是生成其他随机数的基础,U(0,1)型随机数的产生主要是通过乘同余法进行设计生成,乘同余法的计算公式如下所示:
S k + 1 = ( A ⋅ S k ) m o d ( M ) S_{k+1}=(A\cdot S_k)\mathrm{~mod~}(M) Sk+1=(A⋅Sk) mod (M)
在公式中A表示整数常数, mod表示取模运算, M表示较大的整数
通过数论理论能够证明:
当 M = 2 L ( L > 2 ) M=2^{L}(L>2) M=2L(L>2)时,若 A = 4 k + 1 或者 A = 8 k + 3 A=4k+1或者A=8k+3 A=4k+1或者A=8k+3,且 S 0 S_{0} S0为奇数时,可以获得的最长随机数序列长度为 2 L − 2 2^{L-2} 2L−2.
算法实现如下所示:
%% 伪随机数的生成1--乘同余法
%乘同余法
%设定参数
clc;
S0=1;
A=3;
L=4;
M=2^L;
s=S0;
% 循环计算
fprintf('参数A=%d,M=%d,s0=%d的乘除同余法计算结果如下:\n',A,M,S0)
%得到的随机数循环周期为2^(L=2)个
for i =1:2^(L-2)s=mod(A*s,M);fprintf('第%d个随机数为:%d\n', i,s);
end
%如果想产生U(0,1)则需要将求出的s/M即可。
2.2.改进的乘同余法
对于普通的乘同余法只能获得周期为 2 L − 2 2^{L-2} 2L−2的随机数序列,这是远远不够的,所以我们通过添加一个与M互为质数的C来使得乘同余法获得周期为 2 L 2^{L} 2L的随机数数列,这种方法就被称作混合同余法,计算公式如下所示:
S k + 1 = ( A ⋅ S k + C ) m o d ( M ) S_{k+1}=(A\cdot S_k+C)\mathrm{~mod~}(M) Sk+1=(A⋅Sk+C) mod (M)
算法实现如下所示:
%% 伪随机数的生成2--混合乘同余法
%混合乘同余法
%设定参数
clc;
S0=1;
A=3;
L=4;
M=2^L;
s=S0;
C=3;
% 循环计算
fprintf('参数A=%d,M=%d,C=%d,S0=%d的乘除同余法计算结果如下:\n',A,M,C,S0)
%得到的随机数循环周期为2^(L)个
for i =1:2^(L)s=mod(A*s+C,M);fprintf('第%d个随机数为:%d\n', i,s);
end
%如果想产生U(0,1)则需要将求出的s/M即可。
3.产生正态分布的伪随机数
产生正态分布的伪随机数的基本原理:若 Y 1 , Y 2 , Y 3 . . . . . . Y n Y_{1},Y_{2},Y_{3}......Y_{n} Y1,Y2,Y3......Yn 是独立同分布,均值和方差分别为 μ \mu μ和 σ \sigma σ ,且 n n n较大,则 X = Y 1 + Y 2 + Y 3 . . . . . . + Y n X=Y_{1}+Y_{2}+Y_{3}......+Y_{n} X=Y1+Y2+Y3......+Yn 近似于正态分布,且满足 μ x = μ 1 + μ 2 + μ 3 . . . . . + μ n = n μ \mu_{x}=\mu_{1}+\mu_{2}+\mu_{3}.....+\mu_{n}=n\mu μx=μ1+μ2+μ3.....+μn=nμ 及 σ x 2 = σ 1 2 + σ 2 2 + σ 3 2 . . . . . + σ n 2 = n σ 2 \sigma_{x}^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}+\sigma_{3}^{2}.....+\sigma_{n}^{2}=n\sigma_{}^{2} σx2=σ12+σ22+σ32.....+σn2=nσ2,即 x ∈ N ( n μ , n σ 2 ) x\in N( n \mu,n\sigma^{2}) x∈N(nμ,nσ2).
于是正态分布可以由多个U(0,1)来近似。
对于 Y ∈ U ( 0 , 1 ) Y\in U( 0,1) Y∈U(0,1) 来说,对于Y的均值有:
μ y = 1 2 \mu_y=\frac12 μy=21
对于Y的方差,计算如下所示:
σ y 2 = E ( Y 2 ) − ( E ( Y ) ) 2 = ∫ − ∞ + ∞ f ( y ) d y − ( 1 2 ) 2 = ∫ 0 1 y 2 d y − 1 4 = y 3 3 ∣ 1 0 − 1 4 = 1 12 \begin{aligned}\sigma_y^2&=E\color{r}{\left(Y^2\right)-\left(E(Y)\right)^2=\int_{-\infty}^{+\infty}f(y)dy-\left(\frac12\right)^2\\}=\int_0^1y^2dy-\frac14=\frac{y^3}3|\frac10-\frac14=\frac1{12}\end{aligned} σy2=E(Y2)−(E(Y))2=∫−∞+∞f(y)dy−(21)2=∫01y2dy−41=3y3∣01−41=121
令 z = x − μ x σ x z=\frac{x-\mu_x}{\sigma_x} z=σxx−μx,则 z ∈ N ( 0 , 1 ) z\in N(0,1) z∈N(0,1),则z的公式如下所所示:
z = ∑ y i − μ x σ x = ∑ y i − n μ y n σ y 2 = ∑ y i − n 2 n / 12 z=\frac{\sum y_i-\mu_x}{\sigma_x}=\frac{\sum y_i-n\mu_y}{\sqrt{n\sigma_y^2}}=\frac{\sum y_i-\frac n2}{\sqrt{n/_{12}}} z=σx∑yi−μx=nσy2∑yi−nμy=n/12∑yi−2n
一般取n=12,则z的计算公式为:
z = ∑ i = 1 12 y i − 6 ∈ N ( 0 , 1 ) z=\sum_{i=1}^{12}y_i-6\in N(0,1) z=i=1∑12yi−6∈N(0,1)
若想产生服从一般正态分布 N ( μ , σ 2 ) N(\mu,\sigma^{2}) N(μ,σ2) 的随机数x ,则只需产生 z ∈ N ( 0 , 1 ) z\in N(0,1) z∈N(0,1) ,再按公式 x = μ + σ z x=\mu+\sigma z x=μ+σz,即可获得 x ∈ N ( μ , σ 2 ) x\in N(\mu,\sigma^2) x∈N(μ,σ2).
算法实现如下所示(取 n = 1200 n=1200 n=1200计算结果好):
%% 伪随机数的生成3--正态分布方法
%正态分布方法
%设定参数
clc
%生成12个U(0,1)分布的随机数就直接调用包来解决
N=1200;
for i=1:1000Y=rand(N,1);z=(sum(Y)-N*0.5)/sqrt(N/12);fprintf('第%d个属于N(0,1)的随机数为:%.2f\n',i,z)
end
4.基于逆变法产生伪随机数
逆变法产生伪随机数的基本原理是设Y是(0,1)上均匀分布随机变量,F为任意一个连续分布函数,定义随机变量 X = F − 1 ( Y ) X=F^{-1}(Y) X=F−1(Y),则 X具有分布函数F。
证明如下:
F X ( a ) = P { X ≤ a } = P { F − 1 ( Y ) ≤ a } = P { Y ≤ F ( a ) } \begin{aligned}F_X(a)&=P\{X\leq a\}=P\{F^{-1}~(Y)\leq a\}=P\{Y\leq F(a)\}\end{aligned} FX(a)=P{X≤a}=P{F−1 (Y)≤a}=P{Y≤F(a)}
这里 Y ∈ U ( 0 , 1 ) Y\in U( 0,1) Y∈U(0,1) ,有
f ( y ) = 1 , F ( y ) = P { Y ≤ y } = ∫ − ∞ y f ( y ) d y = y f(y)=1,\quad F(y)=P\{Y\leq y\}=\int_{-\infty}^yf(y)dy=y f(y)=1,F(y)=P{Y≤y}=∫−∞yf(y)dy=y
最后能够推出:
F X ( a ) = F ( a ) F_X(a)=F(a) FX(a)=F(a)