温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
作者: 崔颖(中央财经大学)
目录
Source: Survival Rates and Negative Binomial (Francis, 2012)
炎炎夏日,你是否也深受蚊虫叮咬困扰?
本篇推文将帮助你用 Stata 数值模拟预测在给定不同的死亡率(可能被捕食或遭遇意外等)情况下,蚊子的平均存活期限及存活概率。
若给定失败概率和淘汰所需失败次数,则预期存活天数服从负二项分布,这里就是第1次“死亡事件”出现时所需的天数服从负二项分布。根据负二项分布定义及概率密度函数可知,预期存活天数的期望为
其中,
给定若干死亡率:
死亡率=0.01时,
死亡率=0.015时,
死亡率=0.02时,
死亡率=0.025时,
死亡率=0.03时,
死亡率=0.035时,
首先设置 12,000 个观测值。
clear
set obs 12000
将其分为6组,每组200个观测值。6组观测值的死亡率取值分别为0.01,0.015,0.02,0.025,0.03,0.035。
gen mort_rate=mod(_n,6)/200+0.01
tab mort_rate
生成一个哑变量表示蚊子尚且存活:
gen living=1
生成一个连续变量表示蚊子存活的天数:
gen survived=0
连享会计量方法专题……
假设观测了 500 天,看每天还剩多少蚊子存活。
forvalue i=1/500{
*随机模拟伯努利试验,生成临时变量`died`表示特定某一天蚊子是否死亡
qui gen died=rbinomial(1,mort_rate) if living==1
*如果蚊子死亡,则替换其存活状态变量
qui replace living=0 if living==1 & died == 1
*如果蚊子还活着,替换survived变量来表示当前是第几天
qui replace survived=`i' if living==1
*每次循环报告出当前循环是第几天,以及当天蚊子尚存活的比例
qui sum living
di "Round `i' :" r(mean)
*删除临时变量`died`
drop died
}
根据不同的死亡率分组绘制蚊子存活天数的直方图。
hist survived, by(mort_rate)
从图上可以看出,死亡率越高的组别,蚊子存活数量衰减越快。
此外,我们还可以根据不同的死亡率分组看 suivived 变量的描述性统计信息。
bysort mort_rate: sum survived
-----------------------------------------------------------------------
-> mort_rate = .01
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
survived | 2,000 100.424 96.9187 0 500
-----------------------------------------------------------------------
-> mort_rate = .015
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
survived | 2,000 66.1455 67.45336 0 470
-----------------------------------------------------------------------
-> mort_rate = .02
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
survived | 2,000 49.517 50.21277 0 408
-----------------------------------------------------------------------
-> mort_rate = .025
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
survived | 2,000 40.9 40.77116 0 302
-----------------------------------------------------------------------
-> mort_rate = .03
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
survived | 2,000 32.2755 33.19746 0 276
-----------------------------------------------------------------------
-> mort_rate = .035
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
survived | 2,000 28.033 28.43426 0 215
从上表可得,6组观测值存活天数的数值模拟均值均与本文开篇通过负二项分布概率密度函数计算所得的期望近似,验证了结果的可靠性。
本文通过一个有趣的例子介绍了用 Stata 做数值模拟的方法,步骤可以大致概括为:确定分布函数和相关参数,生成观测值,循环模拟随机试验,绘图呈现结果。
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看,所有课程可以随时购买观看。
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
⭕ Stata数据清洗 | 游万海 | 直播, 2 小时,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD