Stata:随机数生成原理与实现

发布时间:2022-01-28 阅读 4314

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:王琳茗 (中山大学)
邮箱wanglm2326@outlook.com


目录


1. 随机数生成器介绍

很多游戏都会设置类似 “开箱” 的环节:比如《炉石传说》里,玩家打开一个卡包,可以获得特定卡牌库中的 5 张卡牌。但究竟是沉睡者伊瑟拉还是鱼人宝宝 (卡牌名称),则完全取决于玩家的运气,或者说计算机生成的随机数。

作为人类,我们或许可以随口说出一串数字。但是计算机并不能 “随口”,它必须有一个科学稳定的随机数来源,才能得到随机数。我们将这个来源称为随机数生成器 (Random number generators,RNG)。常见的随机数生成器有三种:

  • 真随机数生成器 (True Random Number Generator,TRNG),通过物理方法模拟自然界中的随机过程,真随机数无法预测且无周期性;
  • 伪随机数生成器 (Pseudo Random Number Generator,PRNG),通过数学方法生成和真随机数具有相似统计特征的伪随机数。如果能够通过统计检验,就可以被当成真随机数使用;
  • 随机数表法,即用真随机数生成器事先生成大量随机数,存到数据库中,使用时再从库中调用。在 20 世纪早期,这种方法被大量使用,比如 Rand 公司在 1955 年出版了一本《A Million Random Digits with 100,000 Normal Deviates》(百万乱数表)。

随机数表法需要占用大量存储空间,现在已经很少被采用。这里,我们主要介绍前两种方法:真随机数生成器与伪随机数生成器。

1.1 真随机数生成器

第一个真随机数生成器诞生于 1955 年,由我们上文提到的 Rand 公司创造。在 1999 年,Intel 发布 Intel810 芯片组时,就配备了硬件随机数发生器,原理是电阻和振荡器生成的热噪声。目前,大部分芯片厂商都集成了硬件随机数发生器,使用非常方便,而一系列为科研和信息安全设计的真随机数发生器也层出不穷。

1.1.1 基于电路的 TRNG

振荡器采样:如上文中提到的 Intel810RNG 芯片,利用热噪声 (由导体中电子的热震动引起的) 放大后,影响一个由电压控制的振荡器,再通过另一个高频振荡器来收集数据,得到随机数。在 Intel 815E 芯片组的个人电脑上安装 Intel Security Driver (ISD) 后,就可以通过编程读取寄存器获取 RNG 中的随机数。

直接放大电路噪声:直接以热噪声等电路噪声为随机源,通过运算放大,统计一定时间内达到阈值的信号数以此来得到随机数。

电路亚稳态:2010 年,德国研究团队开发出一种真随机数发生器,它使用的计算机内存双态触发器作为随机的一个额外层,触发器可随机的在 1 或 0 状态中切换。在切换之前,触发器处于行为无法预测的 “亚稳态”。在亚稳态结束时,内存中的内容为完全随机。研究人员对一个触发器单元阵列的实验显示,这种方法产生的随机数比传统方法 “随机” 约 20 倍。

混沌电路:混沌电路的输出结果对初始条件很敏感,不可预测,且在 IC 芯片中易集成,可产生效果不错的真随机数。

1.1.2 基于物理的 TRNG

如果你还记得薛定谔那只可怜的猫,你应该知道:在我们打开盒子前,这只猫同时处于两种状态,猫既是死的也是活的。这听上去像是一个悖论,然而微观世界的物理规则确确实实是这样的。微观世界粒子的空间分布和动量是完全不确定的 (即量子力学的不确定性原理),从本质上讲就是真正随机的,因此很适合用来做 TRN。

当然,这并不是说基于经典宏观物理学的 TRNG 就不存在。例如,「RANDOM.ORG」网站从 1998 年开始,在线提供基于大气噪音生成的真随机数。以下是一些基于量子力学的真随机数生成器实例:

  • 克罗地亚计算机科学家发明的 TRNG,全名是 Quantum Random Bit Generator Service (QRBGS),它依赖于半导体光子发散量子物理过程中内在的随机性,通过光电效应检测光子得到随机数;
  • 2010 年,比利时物理学家 S. Pironio 和同事利用纠缠粒子的随机性和非局域性属性创造了真随机数;
  • 2011 年,加拿大渥太华的物理学家 Ben Sussman 利用激光脉冲和钻石创造了真随机数。Sussman 的实验室使用持续几万亿分之一秒的激光脉冲照射钻石,激光进入和出来的方向发生了变化。Sussman 称改变与量子真空涨落的相互作用有关,在量子法则中这种作用是不可知的,他认为这可以用于创造真正的随机数;
  • 2012 年,澳大利亚国立大学的科学家从真空中的亚原子噪音获取随机数,创造了世界上最快的随机数发生器。量子力学中,亚原子对会持续自发的产生和湮灭,通过监听真空内亚原子粒子量子涨落产生的噪音,可以得到真正的随机数。

1.1.3 基于其他因素的 TRNG

除了以上方法,真随机数的来源还可以是击键的间隔时间、鼠标移动速度、特定中断的时间间隔和块 IO 请求的响应时间等。

1.2 伪随机数生成器

基于物理方法的 TRNG,生成的随机数无周期、不可预测、分布均匀。然而,这种随机数生成器技术要求高、生成效率低,难以满足计算机高速计算的需要。因此为了提高数据产生效率,TRNG 的结果往往被用作 PRNG 种子值,并以此生成伪随机的输出序列。

1.2.1 线性同余法

线性同余生成器 (linear congruential generators,LCGs) 是一种很常见的 PRNG,在很多编程语言中都有应用,如 Java 的 java.util.Random 类就是基于这种方法。该方法由 Lehmer (1951) 提出:给定非负整数 ηa、和 c,对于任意整数 z[0],我们都能通过递归公式得到一个整数 z[k] 的序列:

模数符号 mod 表示,z[k] 是 az[k1]+c 被 η 除后的余数。接着,我们将 z[k] 除以 η,就得到了一个伪随机数 u[k] 的序列:

我们可以从中看到伪随机数生成器的两个特点:

伪随机数由种子值决定:给定合适的 ηa、和 c 值后,得到的伪随机数序列完全取决于 z[0],这个初始值称为种子值 (seed value)。如果种子值不变,产生的随机数序列也不会变。这是伪随机数的优点之一:可以重复相同的模拟过程。

值得一提的是,目前 90% 的游戏都是使用内置时钟的时间作为种子值的。当你开箱时,程序会记录你按键时的时间,比如 12 点 34 分 17.745395 秒,使用小数点后的三四位 53 作为种子值。因此,当你总是开出 “蓝天白云” (稀有度很低的卡牌) 时,不妨换个时间试试。

伪随机数是周期性的:因为整数 z[k] 是非负的,并且以 η 为界,伪随机数序列一定会在 η 的周期内循环。为了保证在 32 位计算机上实现无舍入的浮点数计算,在过去,LCG 的参数通常会被设置为:

从而限制了最大周期 η。周期性是所有伪随机数生成器的特性。为了应对这一问题,我们需要使循环周期尽可能长,超过我们所可能用到的最长序列长度。

1.2.2 梅森旋转法

梅森旋转算法 (Mersenne twister) 由松本真和西村拓士 (1997) 提出,其周期可达到 2199371,比可观测宇宙内粒子总数的估计值 (1087) 还要高出上千个数量级。尽管如此长的周期并不必然意味着高质量的伪随机数,但短周期一定不是好的伪随机数。

该方法目前在各种编程语言和库中已普遍存在或作为默认的伪随机数生成器。在 Stata 中,runiform() 函数是所有随机数生成函数的基础,其默认使用的就是 64 位梅森旋转算法。梅森旋转算法主要分为 3 步:

  1. 初始化 n 个状态:根据给定的种子点 x0,通过移位、异或、乘法、加法等操作生成后续的 n1 个状态 x1 到 xn1,bit 位数为 w
  2. 生成伪随机数:根据当前状态,通过移位、与、异或操作生成随机数;
  3. 更新 n 个状态:每生成 n 个随机数后,在生成下一个随机数前,更新状态。
梅森旋转算法示意图
梅森旋转算法示意图

由于梅森旋转法的过程要比线性同余法复杂得多,此处就不详细展开,有兴趣可以阅读 Matsumoto 和 Nishimura (1998)。

2. Stata 中常用的随机数生成方式

如上文所提到,Stata 使用 64 位梅森旋转算法作为默认的随机数发生器。Stata 可以据此生成符合给定分布的伪随机数。

2.1 生成服从均匀分布的随机数

在 Stata 中,用内置函数 runiform() 生成均匀分布随机数。在生成随机数前,可以通过 set seed # 命令设定随机数的种子值。具体代码如下:

. * 生成 5 个 U(0, 1) 的随机数序列
. // 每个序列包含 5 个随机数
. clear
. set obs 5
. forvalues i = 1/5 {
  2.    gen unif`i' = runiform()
  3. }
 
. // 每个随机数都设置相同的种子
. clear
. set obs 5
. forvalues i = 1/5 {
  2.    set seed 233
  3.    gen unif`i' = runiform()
  4. }

2.2 生成服从其他分布的随机数

其他分布的随机数,一般都可在均匀分布随机数的基础上,通过适当变换生成。对于任意给定的连续型随机变量,若分布函数 F 已知,则该随机变量的随机数可以通过 F1(u) 获得。其中 u 为 U(0,1) 的随机数。即连续型随机变量的随机数可以通过对 uU(0,1) 求分布函数的逆运算获得。

对于任意给定的离散型随机变量,若概率分布为 P{X=xj}=Pj,则该随机变量的随机数同样可以通过对 uU(0,1) 求分布函数的逆运算获得。具体方式为:令 UU(0,1)

可以看出 X 的概率分布为:

2.3 Stata 内置的随机数生成函数

Stata 很贴心的内置了很多随机数函数 (几乎所有典型分布都有覆盖),使得我们大多数情况下并不需要 “对 U(0,1) 求分布函数的逆运算” 手动获得服从其他分布的随机数。

注:除下文列出的一些典型内置函数外,更详细的内容可以通过 help random 命令查看。

2.3.1 常见的离散型随机变量

  • 二项分布:rbinomial(n, p)
  • 超几何分布:rhypergeometric(N, K, n)
  • 泊松分布:rpoisson (m)

具体代码如下:

. set obs 10000
. * 二项分布 B(5, 0.5) 的随机数序列
. gen rb = rbinomial(5, 0.5) 
. * 超几何分布 H(10, 3, 5) 的随机数序列
. gen rh = rhypergeometric(10, 3, 5)
. * 泊松分布 Poisson(10) 的随机数序列
. gen rp = rpoisson(10)

2.3.2 常见的连续型随机变量

  • 正态分布:rnormal(m, s)
  • 指数分布:rexponential(b)
  • 卡方分布:rchi2(df)
  • t 分布:rt(df)

具体代码如下:

. clear
. set obs 10000
. * 标准正态分布的随机数序列
. gen rn = rnormal()
. * 指数分布 Exp(10) 的随机数序列
. gen re = rexponential(10)
. * 自由度为 10 的卡方分布的随机数序列
. gen rc = rchi2(10)
. * 自由度为 5 的 t 分布的随机数序列
. gen rt = rt(5)

3. 参考文献

  • 解密随机数生成器 (1)——真随机数生成器 -Link-
  • 林适雨. 随机数不随机[M]. 科学之谜, 2019-06. -Link-
  • 随机数:真随机数和伪随机数 -Link-
  • 伪随机数生成算法-梅森旋转算法介绍 -Link-
  • Lehmer D H. Mathematical methods in large-scale computing units[J]. Annu. Comput. Lab. Harvard Univ., 1951, 26: 141-146. -Link-
  • Holton G A . Value-at-Risk: Theory and Practice[J]. financial risk, 2003. -Link-
  • Matsumoto M, Nishimura T. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator[J]. ACM Transactions on Modeling and Computer Simulation (TOMACS), 1998, 8(1): 3-30. -PDF-

4. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 随机, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,400+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh