温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
作者:陈勇吏 (上海交通大学安泰经济与管理学院)
目录
本文讨论随机数的产生原理,以及随机数的代码实现方式。在介绍随机数之前,我们先简单介绍一下蒙特卡洛模拟和随机数的概念。
假设
但现实生活中,我们感兴趣的
近似计算方法:
上述近似计算方法称为蒙特卡洛模拟方法。即通过从相同的概率分布中多次抽取独立同分布的样本,使用 样本平均值 近似 总体平均值。将每一次抽取到的独立同分布的样本记为
随机变量的分布形式多种多样,包括连续性随机变量、离散型随机变量。其中连续性随机变量包括均匀分布、指数分布等可以写出分布函数形式的随机变量,以及正态分布、卡方分布等只存在概率密度函数形式的随机变量。针对不同的概率分布形式,生成随机数的方式有所不同。但生成任意分布的随机数的基础,是生成服从
产生
然而,数字计算机通常不使用这种方式来产生随机数,而是使用伪随机数来代替随机数。伪随机数的产生原理是取模。即给定一个初始值
Stata 有很多用来产生随机数的内置函数,其中生成均匀分布随机数的函数为runiform()
。在生成随机数前,可以通过set seed #
命令设定随机数的初始值(种子)。具体代码如下:
*===== 代码 =====
clear
set obs 10
forvalues i = 1/3 {
gen unif`i' = runiform()
}
*===== 生成结果 =====
unif1 unif2 unif3
.1387824 .56928813 .32682486
.92432766 .7060023 .30766002
.3676381 .81376687 .21982671
.11971952 .06189653 .69186391
.91902139 .27757641 .04663675
.60349285 .41005086 .60378652
.89752506 .22520046 .52267927
.48974351 .34788238 .23261151
.21933775 .17111017 .92006581
.88946164 .72442917 .39863705
/* 三次模拟,产生三个不同的随机数序列 */
*===== 代码 =====
clear
set obs 10
forvalues i = 1/3 {
set seed 10
gen unif`i' = runiform()
}
*===== 生成结果 =====
unif1 unif2 unif3
.60128311 .60128311 .60128311
.91370436 .91370436 .91370436
.26673195 .26673195 .26673195
.60736577 .60736577 .60736577
.03607367 .03607367 .03607367
.75318004 .75318004 .75318004
.45917266 .45917266 .45917266
.90569935 .90569935 .90569935
.31343658 .31343658 .31343658
.3755541 .3755541 .3755541
/* 设定相同的随机数,三次模拟都产生相同的随机数序列。*/
对于任意给定的连续型随机变量,若分布函数
证得,通过
生成连续型随机变量的随机数函数包括:均匀分布 runiform(a,b)
、正态分布 rnormal(m,s)
、指数分布 rexponential(b)
、卡方分布 rchi2(df)
、t 分布 rt(df)
等。详细内容可以通过 help random
命令查看。
生成随机数的具体代码如下所示:
*========== 代码:生成随机数序列 ==========
* 1.生成随机数序列
clear
set obs 1000
*1.1 U(1,5)均匀分布的随机数序列
gen runif = runiform()
*1.2 指数分布 Exp(5) 的随机数序列
gen exp = rexponential(5)
*1.3 标准正态分布的随机数序列
gen rnormal = rnormal()
*1.4 N(3,2)正态分布的随机数序列
gen rnormal_3_2 = rnormal(3, 2)
*1.5 自由度为 3 的 t 分布的随机数序列
gen rt = rt(3)
*1.6 自由度为 43 的卡方分布的随机数序列
gen rchi2 = rchi2(43)
* 2.随机数序列的分布图
foreach dist in runif exp rnormal rnormal_3_2 rt rchi2 {
histogram `dist', name(`dist', replace) nodraw //kdensity
}
graph combine runif exp rnormal rnormal_3_2 rt rchi2
由于
*========== 代码 ==========
clear
set obs 10
* 1.指数分布 Exp(5) 的随机数序列
set seed 20
gen exp_f = -5 * log(runiform()) //分布函数求逆的方式
set seed 20
gen exp = rexponential(5) //使用内置函数的方式
* 2.均匀分布 U(1,5) 的随机数序列
set seed 20
gen unif_f = 1 + (5-1)*runiform() //分布函数求逆的方式
set seed 20
gen unif = runiform(1,5) //使用内置函数的方式
*========== 生成结果 ==========
exp_f exp unif_f unif
1.6727648 1.6727648 3.8626318 3.8626318
.47869168 .47869168 4.634807 4.634807
1.1330894 1.1330894 4.1889014 4.1889014
7.785189 7.785189 1.8430378 1.8430378
1.0036905 1.0036905 4.2725067 4.2725067
.02883674 .02883674 4.976997 4.976997
4.806301 4.806301 2.5296427 2.5296427
8.2741926 8.2741926 1.7644917 1.7644917
7.1020147 7.1020147 1.9664666 1.9664666
10.250074 10.250074 1.514932 1.514932
/* 可以看出,两种方式得到的随机数序列相同。 */
当感兴趣的随机变量在 Stata 中没有相应的内置函数来生成随机数时,可以通过分布函数求逆的方式计算。
对于任意给定的离散型随机变量,若概率分布为
可以看出
生成离散型随机变量的随机数函数包括:二项分布 rbinomial(n,p)
、超几何分布 rhypergeometric(N,K,n)
、泊松分布 rpoisson(m)
、负二项分布 rnbinomial(n,p)
等。详细内容可以通过 help random
命令查看。
生成随机数的具体代码如下所示:
*========== 代码:生成随机数序列 ==========
* 1.生成随机数序列
clear
set obs 1000
*1.1 二项分布 B(10,0.5) 的随机数序列
gen rbinomial = rbinomial(10,0.5)
*1.2 超几何分布 H(20,4,10) 的随机数序列
/*
其中:20 代表总体数量
4 代表总体中指定种类的数量
10 代表抽取的样本数量
*/
gen rhypergeometric = rhypergeometric(20,4,10)
*1.3 泊松分布 Poisson(20) 的随机数序列
gen rpoisson = rpoisson(20)
*1.4 负二项分布 NB(123,0.2) 的随机数序列
gen rnbinomial = rnbinomial(123,0.2)
* 2.随机数序列的分布图
foreach dist in rbinomial rhypergeometric rpoisson rnbinomial {
histogram `dist', name(`dist', replace) nodraw //kdensity
}
graph combine rbinomial rhypergeometric rpoisson rnbinomial
Stata 提供了 drawnorm
命令,可以从多元正态分布中生成随机数序列。具体语法为
drawnorm newvarlist [, options]
,常用的 option
选项包括:
corr(matrix)
:设定多元正态分布的协方差矩阵mean(vector)
:设定多元正态分布的均值向量,默认为 mean(0)
。
具体规则可以查看 help drawnorm
。Pearson 相关系数和 Spearman 相关系数是衡量变量相关性最常用的两个指标。
Stata 举例如下:
clear
set obs 1000
* 1.生成相关系数矩阵
/*
相关系数矩阵是对称矩阵,矩阵“维数”与待生成的“随机序列个数”一致。
其中:
(第 i 个)主对角元素:代表(第 i 个)随机变量的【方差】,
(第 i 行 j 列 的)非主对角元素:代表(第 i,j)两个随机变量间的【协方差】
*/
matrix C = (1, 0.5 \ 0.5, 1)
matrix list C
*======= Stata 结果 ========
symmetric C[2,2]
c1 c2
r1 1
r2 .5 1
*===========================
* 2.生成二元正态分布样本(生成两个随机数序列:均值为 0,协方差矩阵为 C)
drawnorm x y, corr(C) //means(0,0)选项设定均值向量为(0,0),可省略不写
* 3.查看两个随机序列的相关性
*3.1 Pearson 相关系数
corr x y
*========== Stata 结果 ===========
(obs=1,000)
| x y
-------------+------------------
x | 1.0000
y | 0.5218 1.0000
*=================================
*3.2 Spearman 相关系数
*3.2.1 使用 Pearson 系数 计算 Spearman 相关系数
sort x
gen rank_x = _n
sort y
gen rank_y = _n
corr rank_x rank_y
*========== Stata 结果 ===========
(obs=1,000)
| rank_x rank_y
-------------+------------------
rank_x | 1.0000
rank_y | 0.5011 1.0000
*=================================
*3.2.2 使用 内置函数 计算 Spearman 相关系数
spearman x y
*========== Stata 结果 ===========
Number of obs = 1000
Spearman's rho = 0.5011
Test of Ho: x and y are independent
Prob > |t| = 0.0000
*=================================*
/* 两种方式得到的 Spearman 相关系数 相等,都等于 0.5011。*/
从 Pearson 和 Spearman 两个相关系数的定义及计算公式可以看出:对变量做正的线性变换,变量间的 Pearson 相关系数不变;对变量做非负的函数变换, 变量间的 Spearman 相关系数保持不变。
产生具有相关性的正态分布随机序列,可以使用 Stata 的 drawnorm
命令。产生具有相关性的其他分布的随机数序列没有现成的 Stata 命令,可以先将 drawnorm
得到的正态随机序列转换为均匀分布,然后通过分布函数求逆的方式得到想要的随机数序列。原理如下:
若
* 1.生成均匀分布随机序列
clear
set obs 10000
*1.1 生成两个相关的正态随机序列:x,y
matrix C = (1, 0.5 \ 0.5, 1)
drawnorm x y, corr(C)
*1.2 根据正态分布函数,生成与 x,y 对应的均匀分布序列:x_unif,y_unif
gen x_unif = normal(x)
gen y_unif = normal(y)
*1.3 作图展示生成的随机数序列的分布情况
foreach dist in x y x_unif y_unif {
hist `dist', name(`dist', replace) nodraw
}
graph combine x y x_unif y_unif
/* 可以看出,生成的 x_unif y_unif 两个随机序列在 (0,1) 上均匀分布 */
x_unif
和 y_unif
之间的 Spearman 相关系数,与正态序列 x
和 y
之间的 Spearman 相关系数是一样的。* 1.1 正态序列的 Spearman 值
spearman x y
/*============= Stata 结果 =============
Number of obs = 10000
Spearman's rho = 0.4917
Test of Ho: x and y are independent
Prob > |t| = 0.0000
*======================================*/
* 1.2 生成的均匀分布序列的 Spearman 值
spearman x_unif y_unif
/*============= Stata 结果 =============
Number of obs = 10000
Spearman's rho = 0.4917
Test of Ho: x_unif and y_unif are independent
Prob > |t| = 0.0000
*======================================*/
但 Pearson 相关系数有所不同:
* 2.1 正态序列的 Pearson 值
corr x y
/*============= Stata 结果 =============
(obs=10,000)
| x y
-------------+------------------
x | 1.0000
y | 0.5046 1.0000
*======================================*/
* 2.2 生成的均匀分布序列的 Pearson 值
corr x_unif y_unif
/*============= Stata 结果 =============
(obs=10,000)
| x_unif y_unif
-------------+------------------
x_unif | 1.0000
y_unif | 0.4912 1.0000
*======================================*/
/* 可以看出,两组变量的 Spearson 系数值完全一样,
两组变量的 Pearson 系数值有些微不同 */
invgamma()
和 invexponential()
,带入正态分布产生的均匀分布 伽马 Gamma
随机序列和 指数 Exponential
随机序列:* 1.生成 Gamma 随机序列和 Exponential 随机序列
gen x_gamma = invgammap(4, x_unif)
gen y_exponential = invexponential(1, y_unif)
* 2.做序列的分布图
hist x_gamma, name(x_gamma, replace)
hist y_exponential , name(y_exponential , replace)
x_gamma
和指数随机序列 y_exponential
之间的 Spearman 相关系数,与正态序列 x
和 y
之间的 Spearman 相关系数是一样的。spearman x y
/*============= Stata 结果 =============
Number of obs = 10000
Spearman's rho = 0.4758
Test of Ho: x and y are independent
Prob > |t| = 0.0000
*======================================*/
spearman x_gamma y_exponential
/*============= Stata 结果 =============
Number of obs = 10000
Spearman's rho = 0.4758
Test of Ho: x_gamma and y_exponential are independent
Prob > |t| = 0.0000
*======================================*/
但 Pearson 相关系数有所不同:
spearman x y
/*============= Stata 结果 =============
(obs=10,000)
| x y
-------------+------------------
x | 1.0000
y | 0.4894 1.0000
*======================================*/
spearman x_gamma y_exponential
/*============= Stata 结果 =============
(obs=10,000)
| x_gamma y_expo~l
-------------+------------------
x_gamma | 1.0000
y_exponent~l | 0.4515 1.0000
*======================================*/
如此,drawnorm
命令不仅可以生成 正态
随机序列,也可以进一步通过分布函数求逆的方式得到相关性相同的 其他分布
的随机序列。使用 drawnorm
命令生成具有相关性的 非正态分布
随机序列的 Stata 完整代码如下所示:
clear
set obs 10000
* 1.定义相关系数矩阵
matrix C = ( 1, 0.7, -0.3, 0.2 \ ///
0.7, 1, 0.2, -0.1 \ ///
-0.3, 0.2, 1, 0.3 \ ///
0.2, -0.1, 0.3, 1)
* 2.生成 4 个相关系数矩阵为 C 的正态随机向量
drawnorm x1 x2 x3 x4, corr(C)
* 3.根据 逆分布函数 方法生成四个其他分布的随机序列
/*
四个分布包括:
自由度为 5 的卡方分布:chi2(5),
均值为 5 的泊松分布:poisson(5),
自由度为 5,5 的 F 分布:F(5,5)
在(-5,5)之间的均匀分布:uniform(-5,5)
*/
gen x_chi2 = invchi2(5, normal(x1))
gen x_poisson = invpoisson(5, 1-normal(x2)) // normal(x2) 和 1-normal(x2) 都是 U(0,1)均匀分布的随机序列
gen x_F = invF(5, 5, normal(x3))
gen x_uniform = -5 + normal(x4)*10
* 4.查看四个生成序列的分布图
foreach dist in x_chi2 x_poisson x_F x_uniform {
hist `dist', name(`dist', replace) nodraw
}
graph combine x_chi2 x_poisson x_F x_uniform, imargin(large)
* 5.查看相关系数
mat list C
spearman x?
spearman x_*
corr x?
corr x_*
连享会计量方法专题……
在计量模型中,随机数可以很好的模拟扰动项的随机性,以及我们对数据特征的高质量要求。这有助于我们按照预计的分布特征去抽样,进而验证计量模型的有效性,比较计量模型的优劣。我们将引入一个简单的例子,介绍随机数如何应用于 似不相关回归(SUR) 和 OLS 的模型比较。
我们考虑应用 似不相关回归(SUR) 方法来决定是否开发小行星铂金矿物资源的例子,具体内容参见 Pricing an Asteroid。
在本例中,我们尝试估计三个存在相关性的模型。三个模型具有相同的解释变量: 市场现有铂金数量 (X1),投资者们的平均年龄 (X2),宇宙飞船的新闻报道数量(X3)。被解释变量分别为:铂金价格 (Y1),宇宙飞船的投资金额 (Y2),宇宙飞船的新闻报道数量 (Y3)。
模型设定如下:
. clear
. set obs 10000
. gen x1 = 10 + rnormal()
. gen x2 = 50 + rnormal()
. gen x3 = rpoisson(10)
. matrix C = (12, -5, 5 \ -5, 12, -3 \ 5, -3, 6 )
. drawnorm u1 u2 u3, cov(C)
若模型中的系数矩阵为
则生成被解释变量的 Stata 命令如下:
. gen y1 = 20 + 3*x1 + 0*x2 - 0.1*x3 + u1
. gen y2 = 5 + 3*x1 + 3.5*x2 + 5.2*x3 + u2
. gen y3 = 10 + 3*x1 + 5.3*x2 - 10*x3 + u3
. reg y1 x1 x2 x3
. reg y2 x1 x2 x3
. reg y3 x1 x2 x3
*-不加约束:
. sureg (y1 = x1 x2 x3) (y2 = x1 x2 x3) (y3 = x1 x2 x3)
*-加约束:三个模型的 x1 系数相等
. constraint 1 [y1]x1 = [y2]x1
. constraint 2 [y1]x1 = [y3]x1
. sureg (y1 = x1 x2 x3) (y2 = x1 x2 x3) (y3 = x1 x2 x3), const(1 2)
OLS 估计结果:
SUR 估计结果(不加约束):