Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:陈卓然 (中山大学岭南学院)
邮箱:chenzhr25@mail2.sysu.edu.cn
编者按:本文主要摘译自下文,特此致谢!
Source:Grayling M J, Mander A P. Calculations involving the multivariate normal and multivariate t distributions with and without truncation[J]. The Stata Journal, 2018, 18(4): 826-843. -PDF-
目录
由于绝大多数被收集的数据都被假定为正态分布,多元正态分布在统计学中扮演着重要角色。中心极限定理告诉我们对于满足均值为
与此同时,
譬如,Jonathan (2016) 就提到在结构动态离散选择模型中,可以将最大值的期望分解为多元正态 CDF 的一个线性组合,从而可以提高计算速度。在多元
在本推文中,我们将给出一些命令用以计算在存在变量截断和不存在变量截断条件下 MVN 和 MVT 分布的非退化情形。从实际操作层面来看,这些命令可以应用于如下四个方面:
考虑
以下的计算都是围绕着
其中,
我们将在对变量进行分割后,使用一种拟蒙特卡洛模拟的算法对一种随机化的格进行操作。整个过程将使用 pmvnormal
命令完成。这种算法需要对漂移项 (shifts) 和用于蒙特卡洛积分的样本数量进行具体设定。同时它也需要蒙特卡洛置信因子的值,然而我们提供的默认值就可以进行非常好的估计,因而不需要用户进行过多的调整。尽管漂移项或者样本数量的增加理论上会减少蒙特卡洛积分误差,但是会增大代码的运行时间。
在上述等式中,mvnormalden
,它可以直接对于任意指定的
值得注意的是,对于 Stata 15 的用户,可以使用命令 mvnormalcv()
去估计 lnmvnormalden()
也可以用来进行密度函数的计算。
在多元分析中, 多元正态分布中的等坐标分位数 (equicoordinate quantiles) 也是研究的重点之一。变量
为实现此目的,我们使用 invmvnormal
命令,它将使用 mvnormalcv()
或者 pmvnormal
命令计算上式右边的部分。
最后,为得到分布
其中 rnormal()
来实现。我们将使用 rmvnormal
对这一伪随机分布进行估计。
考虑一个
当然有很多种方式可以去定义有关的 NCMVT 积分,这里我们考虑位移 NCMVT 积分,详见「Multivariate T-Distributions and Their Applications」。这种 NCMVT 积分是伴随线性回归中对回归系数的贝叶斯后验分布中产生的。当
我们可以通过 mvt.mvtden
来估计
我们稍后会介绍命令 invmvt
可以进行上述计算,最后为得到一个位移 rchi2()
命令得到。实现上述效果我们使用 rmvt
命令。
我们进一步考虑截断情形下的多元正态分布和非中心化的多元
其中 tmvnormal
、tmvt
、tmvnormalden
、tmvtden
命令对上述分布和密度进行计算。
进一步使用 invmvnormal
和 invtmvt
命令确定
或者,
最后在命令 rtmvnormal
和 rtmvt
中使用排斥抽样法 (rejection sampling) 产生具有
命令安装和数据下载如下:
net sj 18-4 st0542 // 查看命令列表
net install st0542.pkg, replace // 命令安装
net get st0542.pkg, replace // 下载作者提供的范例 dofiles
命令清单和主要功能如下:
mvnormal
:Multivariate normal Distribution
mvnormalden
:Density of the Multivariate Normal Distribution;pmvnormal
:Distribution function of the Multivariate Normal Distribution;invmvnormal
:Quantiles of the Multivariate Normal Distribution;rmvnormal
:Draw n random vectors from the Multivariate Normal Distribution。mvt
:Multivariate t Distribution
mvtden
:Density of the Multivariate t Distribution;mvt
:Distribution function of the Multivariate t Distribution;invmvt
:Quantiles of the Multivariate t Distribution;rmvt
:Draws n random vectors from the Multivariate t Distribution。tmvnormal
:Truncated Multivariate Normal Distribution
tmvnormalden
:Density of the Truncated Multivariate Normal Distribution;tmvnormal
:Distribution function of the Truncated Multivariate Normal Distribution;invtmvnormal
:Quantiles of the Truncated Multivariate Normal Distribution;rtmvnormal
:Draws n random vectors from the Truncated Multivariate Normal Distribution。tmvt
:Truncated Multivariate t Distribution
tmvtden
:Density of the Truncated (noncentral) Multivariate t Dstribution;tmvt
:Distribution function of the Truncated (noncentral) Multivariate t Distribution;invtmvt
:Quantiles of the Truncated (noncentral) Multivariate t Distribution;rtmvt
:Draws n random vectors from the Truncated Multivariate t Distribution。主要选项的具体含义如下:
x(numlist)
:指定一个分位数向量,在此处求解密度;p(real)
:指定一个概率,real
必须严格介于 0 和 1 之间;delta(numlist)
:指明一个 NCMVT 分布或者一个截断 NCMVT 分布的非中心性参数 lower(numlist)
:指明一个下界向量,用 .
代表 upper(numlist)
:指明一个上界向量,用 .
代表 mean(numlist)
:指明一个 MVN 或者一个截断 MVN 分布的位置参数 sigma(string)
:指明上述的 lowertruncation(numlist)
:指明截断下限向量,用 .
代表 uppertruncation(numlist)
:指明截断上限向量,用 .
代表
为解释 mvnormalden
、mvtden
、tmvnormalden
、tmvden
等命令的用法,我们绘制并比较在给定下列参数情形下的不同密度。
首先在区域
matrix Sigma = J(2, 2, 0.5) + 0.5*I(2)
matrix mvnormalden = J(101, 101, 0)
matrix mvtden = J(101, 101, 0)
matrix tmvnormalden = J(101, 101, 0)
matrix tmvtden = J(101, 101, 0)
local i = 1
foreach x1 of numlist -3(0.06)3{
local j = 1
foreach x2 of numlist -3(0.06)3{
quietly mvnormalden, x(`x1', `x2') mean(0, 0) sigma(Sigma)
matrix mvnormalden[`i', `j'] = r(density)
quietly mvtden, x(`x1', `x2') delta(0, 0) sigma(Sigma) df(1)
matrix mvtden[`i', `j'] = r(density)
quietly tmvnormalden, x(`x1', `x2') mean(0, 0) sigma(Sigma) ///
lowertruncation(-1.5, -1.5) uppertruncation(1.5, 1.5)
matrix tmvnormalden[`i', `j'] = r(density)
quietly tmvtden, x(`x1', `x2') delta(0, 0) sigma(Sigma) ///
lowertruncation(-1.5, -1.5) uppertruncation(1.5, 1.5) df(1)
matrix tmvtden[`i', `j'] = r(density)
local `j++'
}
local `i++'
}
然后将上述参数传递到 twoway contour
。
quietly svmat mvnormalden
quietly svmat mvtden
quietly svmat tmvnormalden
quietly svmat tmvtden
generate row = _n
quietly reshape long mvnormalden mvtden tmvnormalden tmvtden, i(row) j(col)
local opt "aspect(2) ccuts(0 0.03 0.06 0.09 0.12 0.15 0.18 0.21) xtitle('x{subscript:1}') ytitle('x{subscript:2}')"
twoway contour mvnormalden row col, nodraw saving(g1,replace) `opt'
twoway contour mvtden row col, nodraw saving(g2,replace) `opt'
twoway contour tmvnormalden row col, nodraw saving(g3,replace) `opt'
twoway contour tmvtden row col, nodraw saving(g4,replace) `opt'
graph combine g1.gph g2.gph g3.gph g4.gph, common imargin(0 0 0 0) rows(2) cols(2)
由上图可知,
我们考虑二元正态分布,其均值向量和协方差矩阵为:
我们同时也考虑两个截断 MVT 分布的混合,每一个MVT 分布都具有上述的协方差矩阵
另一个具有如下均值和截断条件:
我们额外考虑一些可能的
set seed 2
matrix Sigma = J(2, 2, 0.5) + 0.5*I(2)
matrix gof_pvaluesum_mvn = J(4, 5, 0)
matrix gof_pvaluesum_mix = J(4, 5, 0)
local i = 1
foreach n of numlist 10 25 50 100 {
local j = 1
foreach df of numlist 2 5 10 50 100 {
foreach rep of numlist 1/1000 {
rmvt, delta(0,0) sigma(Sigma) df(`df') n(`n') method(chol)
matrix randsamp = r(rmvt)
quietly svmat randsamp
quietly mvtest normality randsamp1 randsamp2
matrix gof_pvaluesum_mvn[`i',`j'] = gof_pvaluesum_mvn[`i',`j'] + r(p_dh)
drop randsamp1 randsamp2
rtmvt, delta(-2,-2) sigma(Sigma) df(`df') n(`n') method(chol) ///
lowertruncation(-2.5, -2.5) uppertruncation(-1.5, -1.5)
matrix randsamp = r(rtmvt)
rtmvt, delta(2,2) sigma(Sigma) df(`df') n(`n') method(chol) ///
lowertruncation(1.5, 1.5) uppertruncation(2.5, 2.5)
matrix randsamp = (randsamp \ r(rtmvt))
quietly svmat randsamp
quietly mvtest normality randsamp1 randsamp2
matrix gof_pvaluesum_mix[`i',`j'] = gof_pvaluesum_mix[`i',`j']+r(p_dh)
drop randsamp1 randsamp2
}
local `j++'
}
local `i++'
}
matrix gof_pvaluesum_mvn = gof_pvaluesum_mvn/1000
matrix gof_pvaluesum_mix = gof_pvaluesum_mix/1000
matrix rownames gof_pvaluesum_mix ="n = 10" "n = 25" "n = 50" "n = 100"
matrix colnames gof_pvaluesum_mix = "nu = 2" "nu = 5" "nu = 10" "nu = 50" "nu = 100"
mat list gof_pvaluesum_mix
matrix rownames gof_pvaluesum_mvn ="n = 10" "n = 25" "n = 50" "n = 100"
matrix colnames gof_pvaluesum_mvn = "nu = 2" "nu = 5" "nu = 10" "nu = 50" "nu = 100"
mat list gof_pvaluesum_mvn
Mixture of MVT distributions
nu = 2 nu = 5 nu = 10 nu = 50 nu = 100
n = 10 .29478428 .40692869 .45395242 .48513255 .48824155
n = 25 .06537599 .27592039 .38071174 .50969908 .51034664
n = 50 .00839741 .15764171 .31944565 .48992866 .51603237
n = 100 .00014414 .05122232 .22498638 .4728084 .50657935
Mixture of truncated NCMVT distributions
nu = 2 nu = 5 nu = 10 nu = 50 nu = 100
n = 10 .66120812 .65019745 .63979867 .64024316 .64828921
n = 25 .62598445 .62452952 .6175105 .58672504 .58852881
n = 50 .42229776 .45829575 .48254353 .46979639 .44504671
n = 100 .14557859 .19336741 .22840076 .24817368 .23413819
在上述过程中我们对协方差矩阵
从上述结果可以看出,MVT 分布中
Note:产生如下推文列表的 Stata 命令为:
lianxh 分布, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh