线性同余法生成伪随机数

本文作者：刘子艳，河南大学经济学院

本文编辑：温和铭

技术总编：李婷婷

Stata&Python云端课程来啦！

       为了平衡团队运营成本，维系公众号的运营，也与国内动辄数千元的Stata课程缩短差距，我们的网课不得不上调价格，我们决定于11月1日起调价，Python课程的价格调整为249.9元，Stata进阶课程调整为249.9元，Stata基础课程调整为299.9元。大家可以告知一下身边想要购买的小伙伴，欲购从速哦，对报名有任何疑问欢迎在公众号后台和腾讯课堂留言~我们在这篇推文的最后提供了每门课程的课程二维码，大家有需要的话可以直接扫描二维码查看课程详情并进行购买哦~

导读

随机数生成器，顾名思义，作用就是随机产生数字，不能根据已经产生的数字预测下次所产生的数。真正的随机数生成器产生的随机数具有随机性、不可预测性、不可重现性。在现实的应用中，除了关键性的应用（比如在密码学中），人们一般使用真正的随机数，而其它应用往往使用伪随机数就足够了。今天就来为大家介绍一种生成伪随机数的常用方法——线性同余法。

一、初步了解线性同余法

线性同余法(linear congruential generator，LCG)是目前应用广泛的伪随机数生成算法，其基本思想是通过对前一个数进行线性运算并取模从而得到下一个数。具体的设计涉及四个参数：乘数a、增量c、模数m和初始值。

当a=0时为和同余法，当c=0时为乘同余法，c≠0时为混合同余法。乘数、增量和模数的选取可以多种多样，只要保证产生的随机数有较好的均匀性和随机性即可，一般采用 混合同余法。

线性同余法的最大周期是m，但一般情况下会小于m。要使周期达到最大，应该满足以下条件：

(1) c和m互质；

(2) m的所有质因子的积能整除a-1；

(3) 若m是4的倍数，则a-1也是；

(4) a，c，都比m小；

(5) a，c是正整数。

掌握了线性同余法的原理以后，我们可以把该方法作为一种程序设计技能的训练，通过选取线性同余的参数来编制自己的随机数生成器。当然，我们的关键目的在于通过生成的随机数了解如下的特性：

(1) 伪随机数的非随机性；

(2) 参数选取对随机数生成的影响；

(3) 伪随机数的自循环性等特征。

下面我们通过几个小“栗子”说明线性同余法生成的伪随机数的性质。

二、“栗子”相伴，消化不难

（一）线性同余法生成伪随机数1.0

我们用线性同余公式，试图生成100个伪随机数，选取的参数分别是：a=13, m=17, c=0。具体程序如下：

clear all 
cap mkdir E:\线性同余法
cd E:\线性同余法 //切换到你的工作路径
set obs 100 //设置观测值个数为100
local m= 17 //设置参数m=17
local a= 13 //设置参数a=13
gen int y = 1 if _n==1
gen t = _n
replace y=mod(`a'*y[_n-1], `m') if _n>1 //为y循环赋值
gen x = y/`m' //生成伪随机变量x
list 
twoway line x t //绘制x随着t变化的趋势图

为了更直观地说明：当参数太小时，线性同余法生成的伪随机数达到一定的数量以后会出现重复、循环的特点，程序在最后绘制了伪随机数x随着观测值序号t变化的趋势图。

从上面两张图可以很直观地看出，生成的伪随机数从第5个数字开始循环，原因在于我们选取的模数m太小。一个整数除以m的余数最多有m个不同的取值，太小的m会使得生成的随机数m步之内必然会重复。因此，合理地选取较大的参数m至关重要。

（二）线性同余法生成伪随机数2.0

我们选择更大的模数m和乘数a，重复上面的问题。具体程序如下：

clear 
set obs 100
local m= 2^13-13
local a= 7^10-17
gen int y = 1 if _n==1
gen t = _n
replace y=mod(`a'*y[_n-1], `m') if _n>1 //为y循环赋值
gen x = y/`m'
list x
twoway line x t
runtest x //对随机数进行随机性检验

伪随机数x随着观测值序号t变化的趋势图如下图：

从趋势图中和数据集中不难看出，在参数改变之后，伪随机数序列周期变长了，图中的随机数到第95个数字以后才出现循环。说明我们的尝试是有效的。尽管数据的随机性和第一个例子相比有较大改进，但是仍不能满足随机模拟的需要。随机性检验结果如下图所示，检验的P-value为0.11，没有达到需要的显著性水平。

（三）线性同余法生成伪随机数3.0

要达到模拟的需要，我们需要数量级达百万以上的m和a。接下来我们通过一段程序，选择更大的参数生成伪随机数列，观察数据的随机性。

clear
set mem 100m
set obs 1000
local m= 10^12 -11
local a= 427419669081
gen int y = 1 if _n==1
gen t = _n
replace y=mod(`a'*y[_n-1], `m') if _n>1
gen x = y/`m'
twoway scatter x t //绘制伪随机数散点图
runtest x

程序绘制的伪随机数散点图及随机性检验结果如下图所示：

从结果中可以看出，程序生成的1000个随机数并没有和前两个例子一样出现循环，散点图的分布很均匀，同时也通过了随机性检验。

比较这三个小“栗”子，我们很容易理解，参数的选取会直接决定生成的伪随机数的随机性程度，故而在我们运用线性同余法生成随机数时必须认真考虑参数的选择。

今天的“小课堂”到这里就结束啦，不知道这篇推文有没有帮助你更好的理解线性同余法呢？如果有的话，希望可以得到你滴点赞、评论、转发哟~

最后，我们为大家揭秘雪球网(https://xueqiu.com/)最新所展示的沪深证券和港股关注人数增长Top10。

腾讯课堂课程二维码

 对我们的推文累计打赏超过1000元，我们即可给您开具发票，发票类别为“咨询费”。用心做事，不负您的支持！

往期推文推荐
Stata&Python云端课程来啦！Seminar|总统的朋友：政治关联与企业价值
爬虫实战 | 爬取中国天气网Seminar|媒体关联董事对融资和外部治理的影响神奇的组内交叉合并 PDF分章节转TXT并实现可视化——以胡景北知青日记1971至1978年为例
order命令——快速改变变量顺序的利器 Ajax应用场景——以获取雪球网港股代码及公司名称为例
Meta Analysis in Stata17      
cngdf——名义GDP与实际GDP之间的摆渡船
合并，“纵”享新丝滑：frameappend & xframeappend
什么是全局配置项？|从零开始的Pyecharts（二）帮你拿下数据可视化|从零开始的Pyecharts 这十年，《金融研究》的编委和读者偏爱哪些研究话题和文章?
学好这一手，英语词典常在手 

关于我们 

   微信公众号“Stata and Python数据分析”分享实用的Stata、Python等软件的数据处理知识，欢迎转载、打赏。我们是由李春涛教授领导下的研究生及本科生组成的大数据处理和分析团队。

   武汉字符串数据科技有限公司一直为广大用户提供数据采集和分析的服务工作，如果您有这方面的需求，请发邮件到[e-mail]，或者直接联系我们的数据中台总工程司海涛先生，电话：[phone]，wechat: super4ht。海涛先生曾长期在香港大学从事研究工作，现为知名985大学的博士生，爬虫俱乐部网络爬虫技术和正则表达式的课程负责人。

此外，欢迎大家踊跃投稿，介绍一些关于Stata和Python的数据处理和分析技巧。

投稿邮箱：[e-mail]投稿要求：
1）必须原创，禁止抄袭；
2）必须准确，详细，有例子，有截图；
注意事项：
1）所有投稿都会经过本公众号运营团队成员的审核，审核通过才可录用，一经录用，会在该推文里为作者署名，并有赏金分成。

2）邮件请注明投稿，邮件名称为“投稿+推文名称”。
3）应广大读者要求，现开通有偿问答服务，如果大家遇到有关数据处理、分析等问题，可以在公众号中提出，只需支付少量赏金，我们会在后期的推文里给予解答。