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 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ 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]}$ 的序列：

z^{[k]} = (a z^{[k - 1]} + c) mod η

模数符号 $mod$ 表示， $z^{[k]}$ 是 $a z^{[k - 1]} + c$ 被 $η$ 除后的余数。接着，我们将 $z^{[k]}$ 除以 $η$ ，就得到了一个伪随机数 $u^{[k]}$ 的序列：

u^{[k]} = z^{[k]} / η

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

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

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

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

a (η - 1) + c < 2^{53}

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

1.2.2 梅森旋转法

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

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

初始化 $n$ 个状态：根据给定的种子点 $x_{0}$ ，通过移位、异或、乘法、加法等操作生成后续的 $n - 1$ 个状态 $x_{1}$ 到 $x_{n - 1}$ ，bit 位数为 $w$ ；
生成伪随机数：根据当前状态，通过移位、与、异或操作生成随机数；
更新 $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$ 已知，则该随机变量的随机数可以通过 $F^{- 1} (u)$ 获得。其中 $u$ 为 $U (0, 1)$ 的随机数。即连续型随机变量的随机数可以通过对 $u \sim U (0, 1)$ 求分布函数的逆运算获得。

对于任意给定的离散型随机变量，若概率分布为 $P {X = x_{j}} = P_{j}$ ，则该随机变量的随机数同样可以通过对 $u \sim U (0, 1)$ 求分布函数的逆运算获得。具体方式为：令 $U \sim U (0, 1)$ ，

X = {\begin{array}{cl} x_{1}, & if U < P_{1} \\ x_{2}, & if P_{1} < U < P + P_{2} \\ ⋮ \\ x_{j}, & if \sum_{1}^{j - 1} P_{j} < U < \sum_{1}^{j} P_{j} \\ ⋮ \end{array}

可以看出 $X$ 的概率分布为：

P {X = x_{j}} = P {\sum_{1}^{j - 1} P_{j} < U < \sum_{1}^{j} P_{j}} = P

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

专题：Stata绘图
- Stata绘图：随机推断中的系数可视化
专题：Stata程序
- Stata：随机森林算法简介与实现
- Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数
专题：IV-GMM
- Stata-双样本孟德尔随机化法 (MR)-mrrobust
专题：倍分法DID
- 倍分法：DID是否需要随机分组？
专题：内生性-因果推断
- Stata：ritest-随机推断（Randomization Inference）
专题：机器学习
- 机器学习：随机森林算法的Stata实现

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]