加载中…
个人资料
方程equationln
方程equationln
  • 博客等级:
  • 博客积分:0
  • 博客访问:10,320
  • 关注人气:9
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

【钟形曲线那些事(下)】正态随机数的生成

(2014-12-08 08:26:56)
标签:

数学

分类: 有时讲理
距《钟形曲线那些事(中)》,也一年过去了。这种程度的拖延,我也不知该说点啥……为自己还能回来再写一篇而鼓掌吧!
【钟形曲线那些事(下)】正态随机数的生成


在Mathematica中输入这行代码

Random[NormalDistribution[0, 1]]

会返回一个服从标准正态分布的随机数。
随机数可以由随机乱数表生成。但乱数表只能生成简单随机数;生成正态随机数需要依赖别的办法。直接能想到的方法是,模仿乱数表,根据频数做个正态分布的随机数表;表中,负数和正数各一半、落在[-1,1]区间的数有68.2%左右……

当然,一般的方法没那么麻烦。更常见的是下面的算法:
  1. 生成12个服从均匀分布的随机数(简单随机数):x1、x2、……、x12;
  2. 计算u=sum(x1,x2,…,x12)-6;
  3. 计算y=E+D×u;
最后得到的y是均值为E、标准差为D的正态随机数。
看上去有点诡异,但可以用一句话将它说破:算法依据是中心极限定理。

本系列第一篇里,谈过中心极限定理。当时说的是:
多次重复相同的实验,每次间互相独立,并且每次的均值与方差已知,那么,如果实验次数较多,则作为全体看待的实验结果,趋近服从正态分布。

本文的算法中,单次实验不再是服从二点分布的「掷硬币」,而是服从连续的均匀分布的「区间[0,1]上的随机数生成」。若将每一次的实验结果记为x1、x2、x3……,则每一次独立地看,均值都是1/2,方差都是1/12。
多次重复实验,实验结果作一次累加:设结果的累加为y=x1+x2+x3+……,按运算法则,随机变量的和的均值和方差,都累加了相应次数。本算法的重复次数是12,为的是,累加后得到的y——方差是12×(1/12)=1,整齐好看,减少运算量。
下一步就很有概率论的味道了。
若记y=sum(x1,x2,…,x12),按中心极限定理,在计算的意义上,它可以近似地视为均值6、标准差1(方差1²)的正态分布。但是我们要的不是服从N(6,1)的随机数啊,我们要的是服从N(E,D²)的随机数——这就要继续下一步了。处理方法也不难:从概率函数的角度入手,作一些变换即可。
使用上的逻辑是,如果一个实验结果x小于指定数值a的概率为Pr(x≤a),那么记y=kx+c小于指定数值b=ka+c的概率也是Pr(y≤b)=Pr(kx+c≤ka+c)=Pr(x≤a)。
回到算法,假设指定数值a是E+d、则求的是Pr(y≤E+D×u)。这时做点变换,即可使得y*=(y-E)/D是服从均值0、方差1的正态分布量,转为求Pr(y*≤u),数值上等于Pr(y≤E+D×u);因为这是个可逆的变换,可以由计算y*得到y的概率——算法在第二步和第三步,通过这个手段,在N(6,1)、N(0,1)和N(E,D²)之间作转换。
这就是这个算法的原理。
——————————————————————————————
【钟形曲线那些事(下)】正态随机数的生成

——————来个xkcd版的钟形曲线作为文章中间分割线——————


来小小地验算一下这个算法。

使用Mathematica生成正态分布的随机数,画出频率直方图。
某次获得的1000个N(0,1)随机数直方图
【钟形曲线那些事(下)】正态随机数的生成
Histogram[
  RandomVariate[NormalDistribution[0, 1], 1000], Automatic, "PDF", Epilog ->        First@Plot[PDF[NormalDistribution[0, 1], x], {x, -4, 4}, PlotStyle -> Blue]]

换用本文提到的算法,还是生成1000个N(0,1)随机数
【钟形曲线那些事(下)】正态随机数的生成
Histogram[
  Table[Sum[RandomReal[{0, 1}], {i, 1, 12}] - 6, {x, 1, 1000}], Automatic, "PDF", Epilog ->    First@Plot[PDF[NormalDistribution[0, 1], x], {x, -4, 4}, PlotStyle -> Red]]

目测,感觉上结果的符合程度还是可以接受的。如有量化的需要,可用统计学软件测量两个算法结果的「距离」。但我觉得关系不大。大部分统计学软件的正态随机数生成器原理都如此,只是我不太清楚Mathematica的是不是这个算法。(用//FullForm没出来生成器的具体内部结构)



这个系列,三篇文章,写作时间相隔太久了,以致于主题几乎不相关……
下次我避免这种程度的拖延吧。钟形曲线的话题,就这样告终了。

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有