Stata: 蒙特卡洛模拟(Monte Carlo Simulation)没那么神秘
作者:侯新烁 湘潭大学 [编译] (知乎 | 简书 | 码云)
Stata连享会 精彩推文1 || 精彩推文2
资料参考来源: The Stata Blog » Monte Carlo simulations using Stata
- Stata连享会 精彩推文1 || 精彩推文2
通过计算机模拟,从总体抽取大量随机样本的计算方法统称为“蒙特卡罗法(Monte Carlo Method,简记为MC)”,这个名字来源于摩纳哥的蒙特卡洛赌场(Monte Carlo Casino),这是最早使用这个方法的一位美国物理学家的叔叔常去的赌场。蒙特卡罗方法长用来确定统计量的小样本性质。,我们知道许多统计量的精确分布没有解析解,通常的处理方法就是应用大样本理论,用渐近分布来近似真实分布,但现实中的样本容量可能不够大,此时蒙特卡罗法的价值就更加凸显了。
在蒙特卡罗模拟中通常需要设计 loop
语句,可能使用的命令有如 forvalues
、 foreach
、while
; 也可能需要生成一些随机数据,比如 help random_number_functions
、 help simulate
等,读者可自行查阅进行扩展。
在本次推文中,我们将使用服从卡方分布的随机数据进行蒙特卡洛模拟,分别展示随机数据的均值分布,多个随机均值的生成与保存,以及样本标准误差的蒙特卡洛模拟过程。
1.概览
一个估计量的 蒙特卡罗模拟( Monte Carlo simulation ,MCS)是基于特定的 数据生成过程 ( data-generating process, DGP )和样本量,使用模拟的方法去近似估计量的抽样分布的过程。使用MCS可以考察一个估计方法对特定的数据生成过程的估计表现。在这篇推文中,我们将展示如何对 Stata 中的估计量进行 MCS 模拟以及如何解释结果。
大样本理论告诉我们,当真正的 DGP 是一个来自自由度为 1 的 chi^2 分布的随机样本时(标记为 χ 2 ( 1 ) \chi^2(1) χ2(1) ),样本平均值是均值的一个很好的估计量。但一位朋友声称这个估计量对于这个 DGP 不会很好,因为 χ 2 ( 1 ) \chi^2(1) χ2(1) 分布会产生异常值。在这篇文章中,使用 MCS 方法,看是否大样本理论可以很好地用于本 DGP 的 500 个观察样本情形。
2.第一个 MCS
我首先从展示如何从 χ 2 ( k ) \chi^2(k) χ2(k) 分布中抽取一个样本量为 500 的随机样本,以及如何估计平均值和平均值的标准误差开始。
示例1:模拟数据的平均值
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
设置种子 为 12345 ,设置随机数生成器的种子可允许样本结果是可重复的。该随机样本的平均样本平均估计值为 0.91 ,估计的标准误差为 0.055。
若有很多估计,每个估计都来自一个独立的随机样本,我就可以估计估计量的抽样分布均值和标准差。为了获得许多估计,我需要多次重复以下过程:
- 1.从 DGP 中抽取
- 2.计算估算值
- 3.存储估算值
我需要知道如何存储这些估计结果以继续这个过程。我还需要知道如何多次重复该过程以及如何访问 Stata 估计值,鉴于许多读者已经熟悉这些主题,我将这些细节信息附在了附录 I 和 II 中,我们将专注于如何存储如此多的抽样结果。
我想把这些估计结果保存在某处,它们将成为随后可分析数据集的一部分。我将使用命令 postfile
, post
和 postclose
将估计值存储在内存中,并在完成后将所有存储的估计值写入一个数据集。示例 2 对有三次抽样过程的情形进行了说明。
示例2:三次抽奖的估计平均值
set seed 12345
postfile buffer mhat using mcs, replaceforvalues i=1/3 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean ypost buffer (_b[y])}postclose bufferuse mcs, clear
list
命令 postfile buffer mhat using mcs, replace
在内存中创建了一个空间叫做 buffer
,在其中我们可以存储将最终写入数据集的结果。mhat
为估计变量的名称,数据集将保存为 mcs.dta
。关键词 using
将新变量名称与新数据集的名称分开。同时,该例中,定义选项 replace
以将 msc.dta
的任何先前版本替换为此处创建的版本。如图可见,我们可以得到 mhat
的三个结果。
具体而言,此例子使用如下循环,重复了计算均值的过程三次。(更多信息可参阅附录 I )。
forvalues i=1/3 {
}
此外,命令
quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y
分别执行了 删除以前的数据,从 χ 2 ( 1 ) \chi^2(1) χ2(1) 分布中抽取大小为 500 的样本,并估算平均值的操作。命令语句前的 quietly
表示抑制输出。
命令 post buffer (_b[y])
将当前抽取的估计平均值存储在缓冲区中,以便下一次观察 mhat
。命令 postclose buffer
将存储在缓冲区中的内容写入文件 mcs.dta
。命令 use mcs, clear; list
从内存中删除最后一个 χ 2 ( 1 ) \chi^2(1) χ2(1) 样本,读入 msc.dta
数据集,并列出数据集。
示例 3 是示例 2 的修改版本; 增加了抽样数量并对结果进行了描述统计。
示例3:2 000 个估计平均值的均值
set seed 12345
postfile buffer mhat using mcs, replaceforvalues i=1/2000 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean ypost buffer (_b[y])}
postclose buffer
use mcs, clear
summarize
估计值的平均值是估计量的抽样分布均值的估计值,它接近 1.0 的真值。2 000 个估计的样本标准差是估计量的采样分布的标准差的估计量,并且接近真实值 sqrt {sigma^2/N} = sqrt{2/500}
,约为 0.0632 ),其中 (sigma^2)是 χ 2 ( 1 ) \chi^2(1) χ2(1) 随机变量的方差。
3.包含标准误
基于均值统计量报告的标准误差是估计量的抽样分布的标准误差估计值。若大样本分布在做近似估计的抽样分布时表现良好,估计的标准误均值
应当接近这些平均估计值的样本标准差。
为了将估计量的 标准差 与估计的 标准误 的平均值进行比较,修改了示例3,以存储标准误差信息。
示例4:2 000个标准误差的平均值
set seed 12345
postfile buffer mhat sehat using mcs, replaceforvalues i=1/2000 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean ypost buffer (_b[y]) (_se[y])}
postclose buffer
use mcs, clear
summarize
其中,命令 postfile buffer mhat sehat using mcs, replace
在缓冲区为新变量 mhat
和 sehat
腾出空间,并且使用命令
post buffer (_b[y]) (_se[y])
存储的每个估计的均值 mhat
,和估计的标准误 sehat
。(如示例 3 所示,命令 postclose buffer
将存储在内存中的内容写入新数据文件。)
可见,这 2,000 个估计的样本标准偏差为 0.0625 ,接近 估计的标准误均值 0.0630。
你可能会认为我应该写“非常接近”,但 0.0625 距离 0.0630 究竟有多近?老实说,我们不知道这两个数字是否彼此足够接近,因为它们之间的距离并不能自动地告诉我们结果推论有多么可靠。
4.估计拒绝率
在频率统计中,如果 p
值低于某指定数值,我们拒绝零假设。若大样本分布很好的近似于有限样本分布,对真实的零假设检验,其拒绝率应接近指定的数值大小。
为了比较拒绝率在 5% 时的情形,本推文修改了例子 4 来计算和存储是否拒绝零假设的 Wald 检验 相关指标。
测试与真正的零假设。(机制讨论见附录 Ⅲ )
示例5:估计拒绝率
set seed 12345
postfile buffer mhat sehat reject using mcs, replaceforvalues i=1/2000 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean yquietly test _b[y]=1local r = (r(p)<.05)post buffer (_b[y]) (_se[y]) (`r')}postclose buffer
use mcs, clear
summarize
本例中,通过 local r = (r(p)<.05)
识别了每次估计中 p
是否显著。可见, 0.048 的拒绝率非常接近 0.05 的真实大小。
5.总结
在这篇文章中,我展示了如何在 Stata 中执行估计量的 MCS 。我讨论了使用 post
命令来存储大量估计结果的机制,并讨论了如何解释大量估计的均值和标准误均值内涵。同时,还介绍了使用估计的拒绝率来评价给定DGP和样本量条件下的大样本近似抽样分布的可用性。
示例结果证明,抽样均值表现为与大样笨原理预测的估计量均值一致。当然,这一结论并不意味着我朋友对异常值的担忧完全是错误的。其他对异常值更稳健的估计量可能具有更好的性质。以后的推文中将尝试为这种权衡提供一些证据。
附录
附录I:同一过程的多次重复
本附录简要介绍了局域宏以及如何使用它们进行多次重复命令操作; 相关更多详细信息,请参见 [P] macro and [P] forvalues for more details。(可以 help macro
help forvalues
)
我们可以在局域宏中存储和访问字符串信息。如将字符串 hello
存储在本地宏命名值中。命令为 local value "hello"
;为了访问存储的信息,需要修饰局域宏的名称。具体来说,在它前面加上单引号 (实际上是英文状态下的 Tab 按键上面按键符号),并用右单引号跟随它。如访问和展示局域宏 value
中的值,可利用 display
在结果窗中展示。也可以将数字存储为字符串。
左边 ` 右边 ' 合起来 `value'
local value "hello"
display "`value'"local value "2.134"
display "`value'"
要多次重复某些命令,需将它们放在一个 forvalues
循环中。例如,下面的代码重复三次 display
命令。
forvalues i=1/3 {display "i is now `i'"}
上面的示例说明 forvalues
定义了一个局域宏,它包含了值列表中的每个值。在前面的示例中,局域宏的名称是 i
,指定的值是 1/3 = {1,2,3 } 。
附录II:获取估计结果
在 Stata 估计命令之后,您可以通过键入 _b[y]
来访问名为 y
的参数的点估计值,也可以通过键入 _se[y]
来访问其标准误。下面的示例说明了此过程。
示例6:访问估计值
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
display _b[y]
display _se[y]
附录Ⅲ:通过检验计算获取p值
本附录解释了创建指示 Wald 检验是否会在一个特定的数值上拒绝零假设的指标的机理,首先是要生成一些数据,并对真实的零假设进行 Wald 检验。
示例7:Wald测试结果
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
test _b[y]=1
检验的结果被保存在 r()
中。接下来,使用 return list
来查看它们,输入 help return list
以获取详细信息。
示例8:通过测试存储的结果
return list
可以看到,p 值存储在 r(p)
中,接下来生成一个 p 值是否小于 0.05 的 0/1
指示器指标,存储在局域宏 r
中(局域宏的介绍可见附录II)。本文通过展示局域宏包含的信息是 0
来完成说明。
local r = (r(p)<.05)
display "`r'"
代码汇总
drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean yset seed 12345
postfile buffer mhat using mcs, replaceforvalues i=1/3 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean ypost buffer (_b[y])}postclose bufferuse mcs, clear
listset seed 12345
postfile buffer mhat using mcs, replaceforvalues i=1/2000 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean ypost buffer (_b[y])}
postclose buffer
use mcs, clear
summarizeset seed 12345
postfile buffer mhat sehat using mcs, replaceforvalues i=1/2000 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean ypost buffer (_b[y]) (_se[y])}
postclose buffer
use mcs, clear
summarizeset seed 12345
postfile buffer mhat sehat reject using mcs, replaceforvalues i=1/2000 {quietly drop _allquietly set obs 500quietly generate y = rchi2(1)quietly mean yquietly test _b[y]=1local r = (r(p)<.05)post buffer (_b[y]) (_se[y]) (`r')}postclose buffer
use mcs, clear
summarizequi local value "hello"
display "`value'"qui local value "2.134"
display "`value'"forvalues i=1/3 {display "i is now `i'"}drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
display _b[y]
display _se[y]drop _all
set obs 500
set seed 12345
generate y = rchi2(1)
mean y
test _b[y]=1return listlocal r = (r(p)<.05)
display "`r'"
关于我们
- 【Stata 连享会(公众号:StataChina)】由中山大学连玉君老师团队创办,旨在定期与大家分享 Stata 应用的各种经验和技巧。
- 公众号推文同步发布于 CSDN-Stata连享会 、简书-Stata连享会 和 知乎-连玉君Stata专栏。可以在上述网站中搜索关键词
Stata
或Stata连享会
后关注我们。 - 点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
- Stata连享会 精彩推文1 || 精彩推文2
联系我们
- 欢迎赐稿: 欢迎将您的文章或笔记投稿至
Stata连享会(公众号: StataChina)
,我们会保留您的署名;录用稿件达五篇
以上,即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。 - 意见和资料: 欢迎您的宝贵意见,您也可以来信索取推文中提及的程序和数据。
- 招募英才: 欢迎加入我们的团队,一起学习 Stata。合作编辑或撰写稿件五篇以上,即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
- 联系邮件: StataChina@163.com
往期精彩推文
- Stata连享会推文列表1
- Stata连享会推文列表2
- Stata连享会 精彩推文1 || 精彩推文2
Stata: 蒙特卡洛模拟(Monte Carlo Simulation)没那么神秘相关推荐
- correlation 蒙特卡洛_蒙特卡洛模拟法及其matlab案例(Monte Carlo simulation method and its matlab case).doc...
蒙特卡洛模拟法及其matlab案例(Monte Carlo simulation method and its matlab case) 蒙特卡洛模拟法及其matlab案例(Monte Carlo s ...
- 蒙特卡洛积分(Monte Carlo Integration)应用:利用蒙特卡洛积分生成 McBeth表
蒙特卡洛积分(Monte Carlo Integration)应用 蒙特卡洛积分 通常函数f(x)的积分: 可以解释为计算函数曲线下方的面积: 而我们的蒙特卡洛积分则是通过近似的方式来获取一个函数的积 ...
- matlab monte carlo,Monte Carlo Simulation
Monte Carlo simulation is a technique used to study how a model responds to randomly generated input ...
- 蒙特卡洛(Monte Carlo)法求定积分
蒙特卡洛(Monte Carlo)法是一类随机算法的统称.随着二十世纪电子计算机的出现,蒙特卡洛法已经在诸多领域展现出了超强的能力.在机器学习和自然语言处理技术中,常常被用到的MCMC也是由此发展而来 ...
- 蒙特卡洛模拟(Monte Carlo simulation)
1.蒙特卡罗模拟简介 蒙特卡罗模拟,也叫统计模拟,这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的,其基本思想很早以前就被人们所发现和利用.早在17世纪,人们就知道用事 ...
- 蒙特卡罗模拟概述(Monte Carlo Simulation)
Overview 蒙特卡洛(Monte Carlo)方法,或称计算机随机模拟方法,是一种基于"随机数"的计算方法.这一方法源于美国在第二次世界大战中研制原子弹的"曼哈顿计 ...
- 如何用R进行蒙特卡罗模拟(Monte Carlo Simulation with R)
本文所讲的蒙特卡罗模拟是建立在正态分布的基础上.假设我们给定一只股票的初始价格P0, 并且从历史日度数据估计出该股票的每日期望回报率为 mean.logret, 标准差为sd.logret. ( 在估 ...
- 计算机软件及应用stata,蒙特卡洛模拟及其Stata应用实现
蒙特卡洛模拟及其Stata应用实现 出版时间:2015年版 丛编项:海南大学经济管理系列丛书 内容简介 <蒙特卡洛模拟及其Stata应用实现>的第1章是Stata软件基础,主要介绍了Sta ...
- MATLAB蒙特卡洛(Monte Carlo)方法求椭圆面积
MATLAB蒙特卡洛方法求椭圆面积 代码 代码 在某个规定的范围内随机打点,找到满足条件的点,并数一下这些点的数量与总的随机点数量的比,就OK了.关键是设置条件. 代码 clear;clc; n=10 ...
最新文章
- 笔记本建立wifi热点!(主要针对XP)【直播】
- CAN总线基础(三)
- mysql中Invalid default value for 'stime'问题
- Struts ActionForm简单理解
- java 获取内存地址_Java 的多态在 JVM 里原来是这样的
- 2.10 m 个样本的梯度下降
- 关于STVP写保护等级2的问题
- 柯桥托业TOEIC考试和PETS哪个含金量高?
- 华为OD机试题:统计射击比赛成绩
- 噜噜噜啦啦啦啦啦啾啾啾~
- vue+element+上传图片插件
- Python字符串日常练习(基础向)
- oracle 19c ORA-00942: 表或视图不存在 ORA-02063: 紧接着 line
- 企业如何做软文营销推广?
- 从零开始开发一个全栈Web应用实录
- 路径中 / ./ ../ ../../ ../../../ 的含义
- 湖北武汉劳务员报考劳务员的工程劳务管理建筑七大员报考
- 一些做“飞卡智能车”时候的总结
- 【算法1-3】暴力枚举——First Step
- vc实现魔兽3改键程序
热门文章
- c++中为什么要用引用?
- C++修改WIN10的ip、子网掩码、网关、DNS
- Ubuntu18.04.4服务器 安装mysql5.7,navicat远程连接
- VirtualBox使用vboxmanage clonehd实现虚拟机克隆
- 微信小程序开发中,onLoad() 和onShow()的区别
- 一只小蜜蜂(hdoj 2044,动态规划递推)
- html li bd,bd.html
- java rowmapper 通用实现_06.Teambiz中三种特殊RowMapper的使用
- 计算机毕业设计Java校园流浪猫图鉴管理系统的设计与实现(源码+系统+mysql数据库+Lw文档)
- PX4官方动态 | ROS2最新版本Galactic现已支持PX4