Stata：多元正态t截断分布的命令
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
理论 + 实证：从「读懂模型」到「折腾模型」
🎦 理论模型构建专题
📅 2022 年 4 月 23-24 日 (周六-周日)
🔑 郭凯明副教授 (中山大学)
🍓 课程主页：https://gitee.com/lianxh/emodel
作者：陈卓然 (中山大学岭南学院)
邮箱：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-
目录
1. 问题背景
2. 原理
2.1 多元正态分布
2.2 多元 t 分布
2.3 截断 MVN 和 MVT 分布
3. 命令介绍
4. 典例分析
4.1 密度比较
4.2 多元正态分布的拟合优度检验
5. 参考文献
6. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 问题背景
由于绝大多数被收集的数据都被假定为正态分布，多元正态分布在统计学中扮演着重要角色。中心极限定理告诉我们对于满足均值为 ，方差为 ，且独立同分布的随机变量，其任何一组随机样本 在 时， 依分布收敛于 。
与此同时， 分布及其在小样本未知方差的正态分布数据中应用也非常重要。然而在许多统计分析中使用的是多元变量，而非单变量。因此将一维的正态分布和 分布拓展到多维情形 (多元正态分布 (MVN) 和多元 分布 (MVT) ) 变得日益重要。
譬如，Jonathan (2016) 就提到在结构动态离散选择模型中，可以将最大值的期望分解为多元正态 CDF 的一个线性组合，从而可以提高计算速度。在多元 分布中也存在着类似的分解方式。再比如说 Grayling 使用多元正态积分优化了多期临床实验设计等。
在本推文中，我们将给出一些命令用以计算在存在变量截断和不存在变量截断条件下 MVN 和 MVT 分布的非退化情形。从实际操作层面来看，这些命令可以应用于如下四个方面：
生成服从上述四种分布的随机数，以进行蒙特卡洛模拟分析； 计算上述四种分布下的联合密度函数； 计算四种分布下的分布函数； 计算四种分布下的分位数。
2. 原理
2.1 多元正态分布
考虑 维具有非退化的 MVN 分布的随机变量 ，我们使用 来表示。其中 是一个位置参数， 是一个正定的协方差矩阵。
以下的计算都是围绕着 展开。其中 ，以及 。对于 ，$-\infty \leq a_{i}<b_{i} \leq="" \infty$。那么对于="" $\boldsymbol{x}="\left(x_{1}," \ldots,="" x_{k}\right)^{\top}\in="" \mathbb{r}^{k}$，<="" p="">
其中，。
我们将在对变量进行分割后，使用一种拟蒙特卡洛模拟的算法对一种随机化的格进行操作。整个过程将使用 pmvnormal 命令完成。这种算法需要对漂移项 (shifts) 和用于蒙特卡洛积分的样本数量进行具体设定。同时它也需要蒙特卡洛置信因子的值，然而我们提供的默认值就可以进行非常好的估计，因而不需要用户进行过多的调整。尽管漂移项或者样本数量的增加理论上会减少蒙特卡洛积分误差，但是会增大代码的运行时间。
在上述等式中， 是 的概率密度函数。我们稍后会介绍一个命令 mvnormalden，它可以直接对于任意指定的 计算上述密度。
值得注意的是，对于 Stata 15 的用户，可以使用命令 mvnormalcv() 去估计 。除此之外，命令 lnmvnormalden() 也可以用来进行密度函数的计算。
在多元分析中, 多元正态分布中的等坐标分位数 (equicoordinate quantiles) 也是研究的重点之一。变量 的第 100 位等坐标分位数要求我们求解如下等式：
为实现此目的，我们使用 invmvnormal 命令，它将使用 mvnormalcv() 或者 pmvnormal 命令计算上式右边的部分。
最后，为得到分布 的伪随机向量，我们将用到如下事实：
其中 。 可以通过 rnormal() 来实现。我们将使用 rmvnormal 对这一伪随机分布进行估计。
2.2 多元 t 分布
考虑一个 维具有 NCMVT 分布的随机变量 ，我们用 来代表。其中 是一个非中心参数向量， 是一个正定矩阵， 是自由度参数。
当然有很多种方式可以去定义有关的 NCMVT 积分，这里我们考虑位移 NCMVT 积分，详见「Multivariate T-Distributions and Their Applications」。这种 NCMVT 积分是伴随线性回归中对回归系数的贝叶斯后验分布中产生的。当 时，上述 NCMVT 就变为中心化的 MVT 分布。上述积分的重点在于：
我们可以通过 mvt.mvtden 来估计 。为得到等坐标分位数，我们需要解如下方程：
我们稍后会介绍命令 invmvt 可以进行上述计算，最后为得到一个位移 分布的伪随机向量，我们使用如下事实：如果 ，并且 ，那么 ，即具有上述我们所需要的分布。进一步地使用和上述多元正态分布相同的算法得到 的伪随机向量。为了获取伪随机数 ，我们可以使用 rchi2() 命令得到。实现上述效果我们使用 rmvt 命令。
2.3 截断 MVN 和 MVT 分布
我们进一步考虑截断情形下的多元正态分布和非中心化的多元 分布。我们分别使用 和 来表示。其中 ，，并且对于 ，$l_{i}<u_{i}$ 是="" $\boldsymbol{x}$="" 的上下截断点。上述变量的分布函数为：、<="" p="">
其中 和 是其各自的密度。我们可以使用 tmvnormal、tmvt、tmvnormalden、tmvtden 命令对上述分布和密度进行计算。
进一步使用 invmvnormal 和 invtmvt 命令确定 的值，使得下式成立：
或者，
最后在命令 rtmvnormal 和 rtmvt 中使用排斥抽样法 (rejection sampling) 产生具有 和 分布的随机变量。简单来说，我们是从 进行随机抽样。如果抽取的样本 对于 满足 ，那么我们就保留这次抽样结果，否则我们就拒绝这次结果，重复上述过程。
3. 命令介绍
命令安装和数据下载如下：
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 分布的位置参数 ，在 MVN 下就是分布的期望值；sigma(string)：指明上述的 矩阵。这个矩阵必须是对称且正定的，在 MVN 下该矩阵就是变量的协方差矩阵，然而在 NCMVN 下这一结论不成立；lowertruncation(numlist)：指明截断下限向量，用.代表 ；uppertruncation(numlist)：指明截断上限向量，用.代表 。
4. 典例分析
4.1 密度比较
为解释 mvnormalden、mvtden、tmvnormalden、tmvden 等命令的用法，我们绘制并比较在给定下列参数情形下的不同密度。
首先在区域 循环一系列值，在每一个网格点上我们计算一个 MVN、MVT、截断的 MVN、以及截断的 MVT。
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 分布和 MVN 分布更相近，这也解释了低自由度 MVT 分布的尾部概率要更大。与非截断相比，截断会增加在区域 内每一点的密度。
4.2 多元正态分布的拟合优度检验
我们考虑二元正态分布，其均值向量和协方差矩阵为：
我们同时也考虑两个截断 MVT 分布的混合，每一个MVT 分布都具有上述的协方差矩阵 ，但是其中一个具有如下的均值和截断条件：
另一个具有如下均值和截断条件：
我们额外考虑一些可能的 值。然后对于不同的 值，我们从 MVT 分布中抽取 个样本，这一抽取重复 1000 次，之后计算出 Doornik-Hansen 检验中的 值平均数。这一过程同样对于截断的 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
在上述过程中我们对协方差矩阵 进行了 Cholesky 分解。对于绝大多数分布而言，分解方法的选择对于生成随机变量的速度和质量影响不大，但是当维度 很大时，分解方法的选择就需要非常谨慎了。
从上述结果可以看出，MVT 分布中 的增加会增加 值的平均值，这一点对于所有的 都成立。这是因为随着 的增加，MVT 分布逐渐接近相应的 MVN 分布。而对于截断 NCMVT 变量的混合，平均的 值总是高于相应的 MVT 分布。在任何情形下 Doornik-Hansen 的检验都不会在 5% 的水平下显著，这意味着平均而言我们无法正确拒绝多元正态性的原假设。
5. 参考文献
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- Cappellari L, Jenkins S P. Calculation of multivariate normal probabilities by simulation, with applications to maximum simulated likelihood estimation[J]. The Stata Journal, 2006, 6(2): 156-189. -PDF- Eggleston J. An efficient decomposition of the expectation of the maximum for the multivariate normal and related distributions[J]. Journal of Econometrics, 2016, 195(1): 120-133. -PDF- Kotz S, Balakrishnan N, Johnson N L. Continuous multivariate distributions, Volume 1: Models and applications[M]. John Wiley & Sons, 2004. Kotz S, Nadarajah S. Multivariate t-distributions and their applications[M]. Cambridge University Press, 2004. Mazza C, Benaim M. Stochastic dynamics for systems biology[M]. CRC Press, 2014. Patel J K, Read C B. Handbook of the normal distribution[M]. CRC Press, 1996.
连享会 · 文本分析 | 爬虫 | 机器学习
6. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 分布, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 品头论足-distcomp：组间分布差异检验 reldist-相对分布：分布差异分析和检验 专题：Stata程序 Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数 专题：回归分析 Stata：局部回归分布估计量-lpdensity Stata：用负二项分布预测蚊子存活率
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【**百度一下：**连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。