高等数理统计实验报告

《高等数理统计》实验报告一、实验目的不同分布的随机数生成原理及其算法实现二、实验内容1.针对均匀分布的随机数生成原理及算法实现2.以指数分布为例的反函数方法的随机数生成3.针对正态分布的随机数生成原理及算法实现4.基于Box–Muller算法的随机数生成及算法实现5.基于舍选抽样法算法的随机数生成及算法实现三、实验过程1.均匀分布随机数生成方法一：迭代取中法原理：这里在迭代取中法中介绍平方取中法 , 其迭代式如下 :X n+1=X n210s(mod102s)R n+1=X n+1/102s其中，X n+1是迭代算子，而R n+1则是每次需要产生的随机数。

第一个式子表示的是将X n平方后右移 s 位，并截右端的 2s 位。

而第二个式子则是将截尾后的数字再压缩 2s 倍，显然 :0≤R n+1≤1。

注：迭代取中法有一个不良特性就是若初始随机数参数 s 以及初始值 X0选择不恰当，最后结果容易退化成0.附结果验证：（1）当初始值为123456，s=2时，部分结果如下，所有结果均不退化为0（2）当初始值为12345，s=2时，部分结果如下：注：当迭代次数达到48次时，随机数则退化为0（3）当初始值为12345，s=1时，部分结果如下：注：当迭代次数达到4次时，随机数则退化为0方法二：乘同余法X n+1= λ∗X n(mod M)R n+1=X n/M其中，X n+1是迭代算子，而R n+1则是每次需要产生的随机数，λ为参数；注：这里的参数选取是需要一定理论基础的，否则所产生的随机数的周期将较小，相关性会较大。

附图：当迭代次数为300次时，p值为0.54：当迭代次数为1000次时，p值几乎为1方法三：混合同余法（又称线性同余法LCG ）混合同余法是加同余法和乘同余法的混合形式,其迭代式如下 X n+1=(λ∗X n +μ)(mod M) R n+1=X n /M其中，X n+1是迭代算子，而 R n+1 则是每次需要产生的随机数，λ 为参数； 一般而言，高 LCG 的 m 是2的指数次幂(一般2^32或者2^64)，因为这样取模操作截断最右的32或64位就可以了在 M =2q 的条件下，参数 λ,μ,X 0按如下方式选取： λ=2c +1,c 取 q/2 附近的数μ=(12+√3)/MX 0为任意非负整数注：参数 λ,μ, X 0 选取直接影响了伪随机数产生的质量该随机数的生成方法在某一周期内成线性增长的趋势，但是在大多数场合，这种极富“规律”型的随机数是不适用的使用场合：LCG 不能用于随机数要求高的场合，例如不能用于Monte Carlo 模拟，不能用于加密应用。

不过有些场合LCG 有很好的应用，例如内存很紧张的嵌入式中，电子游戏控制台用的小整数，使用高位可以胜附结果验证：当迭代次数为1000时，p 值几乎为12.反函数法：采用概率积分变换原理,对于随机变量X的分布函数. F(X)可以求其反函数，得:X i=G(R i )其中, R i为一个[0,1]区间内的均匀分布的随机变量.当 F(X)较简单时，求解较易，当F(X)较复杂时，可能需要用到较为复杂的变换技巧。

这里以指数分布为例，则指数分布的 λ 取1.附结果及 p 值：3.利用中心极限定理生成正态分布随机数独立同分布、且数学期望和方差有限的随机变量序列的标准化和以标准正态分步为极限。

若ξ1,ξ2,⋯,ξn是一串独立同分布的随机变量序列，且有E(ξ)=m,D(ξk)=σ2标准化随机变量和ξk=∑(k) nk=1σ√n若0<σ2<∞，则有lim n→∞P(ξn<x)=1√2πe−t2/2dtx−∞又ξk~U̅(0,1)所以E(ξk)=12,σ2=112所以n=12附结果：4.Box –Muller 算法 原理：Box –Muller 算法实际上是依据瑞利分布来求标准正态分布的反函数。

我们知道标准正太分布的反函数是求不了的，但标准正态分布经过极坐标变换后却是可以求得反函数的 ，有2种形式：（1）基本形式：用 (0,1) 均匀分布随机数，需要计算三角函数 sin 和 cos 值 （2）极坐标形式：用(−1,1)均匀分布随机数，且不需要计算三角函数 我们计算积分I =∫e −x2/2dx +∞−∞,可以先取它的平方并用极坐标表示，I 2=∫∫e−(x 2+y 2)/2+∞−∞dx +∞−∞dy =∫∫re r2/2drdθ∞02π可以发现极角 θ 服从(0,2π)均匀分布，径向距离 r 服从 re r 2/2分布函数（即 r 2 服从 χ2分布），如果将 r 的累积分布函数写出来，即F (x )= ∫re r2/2dr x0=1−e −x2/2我们就可以通过 F 的逆函数 F −1(u ) ,将均匀分布 u~U (0,1) 映射到目标分布，即 x =√−2ln(1−u),也等价于x =√−2ln (u)（注意到如果 u~U (0,1)，那么 1−u 也是），则 产生相互独立的服从[0,1]均匀分布的随机数U 1，U 2，令α=(−2lnU 1)1/2cos 2πU 2 β=(−2lnU 1)1/2sin 2πU 2则 α，β是相互独立的 N (0,1) 随机数 其中，由(−2lnU 1)1/2生成服从Rayleigh 分布的随机变量 由2πU 2生成服从（0,2π）的均匀分布的随机变量 附结果验证：5.舍选抽样法(Rejection Method)用反函数法生成随机数时，如果求不出 F−1(x)的解析形式或者F(x)就没有解析形式，则可以用 F−1(x)的近似公式代替。

可以考虑由Von Neumann提出的舍选抽样法。

原理：构造一个新的密度函数 q(x)，使它的形状接近 p(x)，并选择一个常数 k 使得kq(x)≥p(x)对于 x 定义域内的值都成立。

对应下图，首先从分布q(z)中生成随机数 z0，然后按均匀分布从区间[ 0, kp(z0)]生成一个随机数 u0。

如果 u0>p(z0),则拒绝该样本,否则保留.下图中灰色区域内的点都要舍弃。

由于随机点 u0只出现在曲线 kq(x)之下，且在 q(x)较大处出现次数较多，从而大大提高了采样效率。

显然 q(x)形状越接近 p(x)，则采样效率越高。

算法过程：1.从分布 q(x)中随机取一个数 x2.从均匀分布[0,1]中随机抽取一个数 y3.如果 y <p(x)/kq(x),则接受该样本(x,y),否则就拒绝该样本。

4.重复上述步骤附结果图：附python代码：(1)import numpy as npdef function(x_n,s,iteration):for i in range(iteration):x_n=((x_n**2)/np.power(10,s))%np.power(10,2*s)R_n=x_n/np.power(10,2*s)print(R_n)def main(x_n,iteration,s=2):function(x_n,s,iteration)if __name__ == '__main__':main(x_n=123456,iteration=100,s=2)(2)import numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import kstestset=[]def function(x_n,iteration,lambd,M,epsilon=1e-6):for i in range(iteration):x_n=(x_n*lambd)%MR_n=x_n/Mprint(R_n)set.append(R_n)def main():lambd=np.power(5,5)M=np.power(2,12)function(x_n=12345,iteration=300,lambd=lambd,M=M)if __name__ == '__main__':main()plt.hist(set, bins=10)plt.show()plt.plot(np.arange(1, 301, 1), set)plt.show()p=kstest(set,'uniform')print(p)(3)import numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import kstestdef function(x_n,iteration,lambd,M,mu,epsilon=1e-6):for i in range(iteration):x_n=(x_n*lambd+mu)%MR_n=x_n/Mprint(R_n)def main():lambd=np.power(2,16)+1M=np.power(2,32)mu=(0.5+np.sqrt(3)/6)/Mfunction(x_n=12345,iteration=1000,lambd=lambd,M=M,mu=mu) if __name__ == '__main__':main()plt.hist(set, bins=10)plt.show()plt.plot(np.arange(1, 1001, 1), set)plt.title('迭代次数为1000次')plt.show()p=kstest(set,'uniform')print(p)(4)import numpy as npnp.random.seed(666)#随机数种子n=100y=np.random.uniform(0,1,n)#(0,1)均匀分布随机数def exponential_random_numbers(Lambda):try:x=-1/Lambda*np.log(1-y)return xexcept:return np.zeros(len(y))np.set_printoptions(suppress=True)ern=exponential_random_numbers(Lambda=1)s=np.var(ern)print('方差var:',s)(5)import matplotlib.pyplot as pltimport seaborn as snssns.distplot(ern,bins=20,kde=True, fit=expon)plt.title('指数分布随机数直方图，n=1000,var=%.4f' %s) plt.ylabel('频率')plt.xlabel('x')plt.show()(6)import numpy as npfrom scipy.stats import normnorm_num=[]n=10000for i in range(n):s=0for n in range(12):x=np.random.random()s=s+xy=s-6norm_num.append(y)print(len(norm_num))import matplotlib.pyplot as pltimport seaborn as snssns.distplot(norm_num,bins=20,kde=False,fit=norm)plt.title('标准正态分布随机数直方图，n=10000')plt.ylabel('频率')plt.xlabel('x')plt.show()from scipy.stats import kstestp=kstest(norm_num,'norm')print(p)(6)import numpy as npimport mathfrom scipy.stats import normnorm_nums=[]for i in range(1000):x=np.random.random()y = np.random.random()z=(-2*math.log(x))**(1.0/2)*math.cos(2*math.pi*y) norm_nums.append(z)import matplotlib.pyplot as pltimport seaborn as snssns.distplot(norm_nums,bins=20,kde=False,fit=norm) plt.title('Box–Muller算法，n=10000')plt.ylabel('频率')plt.xlabel('x')plt.show()from scipy.stats import kstestp=kstest(norm_nums,'norm')print(p)(7)from scipy.stats import kstestfrom scipy import *import matplotlib.pyplot as pltp = lambda x: np.exp(-x) # our distributiong = lambda x: 1/(x+1)# our proposal pdfCDFg = lambda x: np.log(x +1)xmin = 0xmax = 8umin = CDFg(xmin)umax = CDFg(xmax)N = 10000accepted = 0samples = np.zeros(N)count = 0while (accepted < N):u = np.random.uniform(umin, umax)xproposal = np.exp(u) - 1y = np.random.uniform(0, 1)if y < p(xproposal)/g(xproposal):samples[accepted] = xproposalaccepted += 1count +=1hinfo = np.histogram(samples,30)plt.hist(samples,bins=30)xvals=np.linspace(xmin, xmax, 1000)plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)') plt.legend()plt.show()p=kstest(samples,'expon') print(p)。