Stata：多元正态t截断分布的命令

发布时间：2022-03-12 阅读 1478

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：陈卓然 (中山大学岭南学院)
邮箱：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}$ ，且独立同分布的随机变量，其任何一组随机样本 $X_{1}, . . ., X_{n}$ 在 $n \to \infty$ 时， $\sqrt{n} (\sum_{i = 1}^{n} X_{i} / n - μ)$ 依分布收敛于 $N (0, σ^{2})$ 。

与此同时， $t$ 分布及其在小样本未知方差的正态分布数据中应用也非常重要。然而在许多统计分析中使用的是多元变量，而非单变量。因此将一维的正态分布和 $t$ 分布拓展到多维情形 (多元正态分布 (MVN) 和多元 $t$ 分布 (MVT) ) 变得日益重要。

譬如，Jonathan (2016) 就提到在结构动态离散选择模型中，可以将最大值的期望分解为多元正态 CDF 的一个线性组合，从而可以提高计算速度。在多元 $t$ 分布中也存在着类似的分解方式。再比如说 Grayling 使用多元正态积分优化了多期临床实验设计等。

在本推文中，我们将给出一些命令用以计算在存在变量截断和不存在变量截断条件下 MVN 和 MVT 分布的非退化情形。从实际操作层面来看，这些命令可以应用于如下四个方面：

生成服从上述四种分布的随机数，以进行蒙特卡洛模拟分析；
计算上述四种分布下的联合密度函数；
计算四种分布下的分布函数；
计算四种分布下的分位数。

2. 原理

2.1 多元正态分布

考虑 $k$ 维具有非退化的 MVN 分布的随机变量 $X = {(X_{1}, \dots, X_{k})}^{⊤}$ ，我们使用 $X \sim N_{k} (δ, Σ)$ 来表示。其中 $δ \in R^{k}$ 是一个位置参数， $Σ \in M^{k \times k}$ 是一个正定的协方差矩阵。

以下的计算都是围绕着 $Φ_{k} (a, b, δ, Σ)$ 展开。其中 $a = {(a_{1}, \dots, a_{k})}^{⊤}$ ，以及 $b = {(b_{1}, \dots, b_{k})}^{⊤}$ 。对于 $i = 1, \dots, k$ ，$-\infty \leq a_{i}

\begin{aligned} Φ_{k} (a, b, δ, Σ) & = \int_{a_{1}}^{b_{1}} \dots \int_{a_{k}}^{b_{k}} ϕ_{k} (x, δ, Σ) d x \\ = \int_{a_{1}}^{b_{1}} \dots \int_{a_{k}}^{b_{k}} \frac{1}{\sqrt{(2 π)^{k} | Σ |}} \exp {- \frac{1}{2} (x - δ)^{⊤} Σ^{- 1} (x - δ)} d x \end{aligned}

其中， $| Σ | = \det (Σ)$ 。

我们将在对变量进行分割后，使用一种拟蒙特卡洛模拟的算法对一种随机化的格进行操作。整个过程将使用 pmvnormal 命令完成。这种算法需要对漂移项 (shifts) 和用于蒙特卡洛积分的样本数量进行具体设定。同时它也需要蒙特卡洛置信因子的值，然而我们提供的默认值就可以进行非常好的估计，因而不需要用户进行过多的调整。尽管漂移项或者样本数量的增加理论上会减少蒙特卡洛积分误差，但是会增大代码的运行时间。

在上述等式中， $ϕ_{k} (x, δ, Σ)$ 是 $X$ 的概率密度函数。我们稍后会介绍一个命令 mvnormalden，它可以直接对于任意指定的 $x, δ, Σ$ 计算上述密度。

值得注意的是，对于 Stata 15 的用户，可以使用命令 mvnormalcv() 去估计 $Φ_{k} (a, b, δ, Σ)$ 。除此之外，命令 lnmvnormalden() 也可以用来进行密度函数的计算。

在多元分析中, 多元正态分布中的等坐标分位数 (equicoordinate quantiles) 也是研究的重点之一。变量 $X$ 的第 100 位等坐标分位数要求我们求解如下等式：

p = Φ_{k} {(- \infty, \dots, - \infty)^{⊤}, (q, \dots, q)^{⊤}, δ, Σ}

为实现此目的，我们使用 invmvnormal 命令，它将使用 mvnormalcv() 或者 pmvnormal 命令计算上式右边的部分。

最后，为得到分布 $N_{k} (δ, Σ)$ 的伪随机向量，我们将用到如下事实：

{R R}^{⊤} = Σ, x = δ + R z

其中 $z = {(z_{1}, \dots, z_{k})}^{⊤}$ 。 $z_{i} \sim N (0, 1)$ 可以通过 rnormal() 来实现。我们将使用 rmvnormal 对这一伪随机分布进行估计。

2.2 多元 t 分布

考虑一个 $k$ 维具有 NCMVT 分布的随机变量 $X = {(X_{1}, \dots, X_{k})}^{⊤}$ ，我们用 $X \sim T_{k} (δ, Σ, ν)$ 来代表。其中 $δ \in R^{k}$ 是一个非中心参数向量， $Σ \in M^{k \times k}$ 是一个正定矩阵， $ν \in R^{+}$ 是自由度参数。

当然有很多种方式可以去定义有关的 NCMVT 积分，这里我们考虑位移 NCMVT 积分，详见「Multivariate T-Distributions and Their Applications」。这种 NCMVT 积分是伴随线性回归中对回归系数的贝叶斯后验分布中产生的。当 $δ = 0$ 时，上述 NCMVT 就变为中心化的 MVT 分布。上述积分的重点在于：

\begin{aligned} Ψ_{k} (a, b, δ, Σ, ν) & = \int_{a_{1}}^{b_{1}} \dots \int_{a_{k}}^{b_{k}} ψ_{k} (x, δ, Σ, ν) d x \\ = \int_{a_{1}}^{b_{1}} \dots \int_{a_{k}}^{b_{k}} \frac{Γ {(ν + k) / 2}}{Γ (ν / 2) \sqrt{(ν π)^{k} | Σ |}} {1 + \frac{1}{ν} (x - δ)^{⊤} Σ^{- 1} (x - δ)}^{- (ν + k) / 2} d x \end{aligned}

我们可以通过 mvt.mvtden 来估计 $ψ_{k} (x, δ, Σ, ν)$ 。为得到等坐标分位数，我们需要解如下方程：

p = Ψ_{k} {(- \infty, \dots, - \infty)^{⊤}, (q, \dots, q)^{⊤}, δ, Σ, ν}

我们稍后会介绍命令 invmvt 可以进行上述计算，最后为得到一个位移 $T_{k} (δ, Σ, ν)$ 分布的伪随机向量，我们使用如下事实：如果 $y \sim N_{k} (0, Σ)$ ，并且 $v \sim χ_{ν}^{2}$ ，那么 $x = δ + y \sqrt{\frac{ν}{v}}$ ，即具有上述我们所需要的分布。进一步地使用和上述多元正态分布相同的算法得到 $y$ 的伪随机向量。为了获取伪随机数 $ν$ ，我们可以使用 rchi2() 命令得到。实现上述效果我们使用 rmvt 命令。

2.3 截断 MVN 和 MVT 分布

我们进一步考虑截断情形下的多元正态分布和非中心化的多元 $t$ 分布。我们分别使用 $X \sim T N_{k} (δ, Σ, l, u)$ 和 $X \sim T T_{k} (δ, Σ, ν, l, u)$ 来表示。其中 $l = {(l_{1}, \dots, l_{k})}^{⊤} \in R^{k}$ ， $u = {(u_{1}, \dots, u_{k})}^{⊤} \in R^{k}$ ，并且对于 $i = 1, . . ., k$ ，$l_{i}

\begin{aligned} Φ_{k}^{'} (a, b, δ, Σ, l, u) & = \int_{max (a_{1}, l_{1})}^{min (b_{1}, u_{1})} \dots \int_{max (a_{k}, l_{k})}^{min (b_{k}, u_{k})} ϕ_{k}^{'} (x, δ, Σ, l, u) d x \\ = \int_{max (a_{1}, l_{1})}^{min (b_{1}, u_{1})} \dots \int_{max (a_{k}, l_{k})}^{min (b_{k}, u_{k})} \frac{ϕ_{k} (x, δ, Σ)}{Φ_{k} (l, u, δ, Σ)} d x \\ Ψ_{k}^{'} (a, b, δ, Σ, ν, l, u) & = \int_{max (a_{1}, l_{1})}^{min (b_{1}, u_{1})} \dots \int_{max (a_{k}, l_{k})}^{min (b_{k}, u_{k})} ψ_{k}^{'} (x, δ, Σ, ν, l, u) d x \\ = \int_{max (a_{1}, l_{1})}^{min (b_{1}, u_{1})} \dots \int_{max (a_{k}, l_{k})}^{min (b_{k}, u_{k})} \frac{ψ_{k} (x, δ, Σ, ν)}{Ψ_{k} (l, u, δ, Σ, ν)} d x \end{aligned}

其中 $ϕ_{k}^{'} (x, δ, Σ, l, u)$ 和 $ψ_{k}^{'} (x, δ, Σ, ν, l, u)$ 是其各自的密度。我们可以使用 tmvnormal、tmvt、tmvnormalden、tmvtden 命令对上述分布和密度进行计算。

进一步使用 invmvnormal 和 invtmvt 命令确定 $q$ 的值，使得下式成立：

p = Φ_{k}^{'} {(- \infty, \dots, - \infty)^{⊤}, (q, \dots, q)^{⊤}, δ, Σ, l, u}

或者，

p = Ψ_{k}^{'} {(- \infty, \dots, - \infty)^{⊤}, (q, \dots, q)^{⊤}, δ, Σ, ν, l, u}

最后在命令 rtmvnormal 和 rtmvt 中使用排斥抽样法 (rejection sampling) 产生具有 $T N_{k} (δ, Σ, l, u)$ 和 $T T_{k} (δ, Σ, ν, l, u)$ 分布的随机变量。简单来说，我们是从 $N_{k} (δ, Σ)$ 进行随机抽样。如果抽取的样本 $X_{i}$ 对于 $i = 1, . . ., k$ 满足 $l_{i} \leq X_{i} \leq u_{i}$ ，那么我们就保留这次抽样结果，否则我们就拒绝这次结果，重复上述过程。

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)：指明一个下界向量，用 . 代表 $- \infty$ ；
upper(numlist)：指明一个上界向量，用 . 代表 $+ \infty$ ；
mean(numlist)：指明一个 MVN 或者一个截断 MVN 分布的位置参数 $δ$ ，在 MVN 下就是分布的期望值；
sigma(string)：指明上述的 $Σ$ 矩阵。这个矩阵必须是对称且正定的，在 MVN 下该矩阵就是变量的协方差矩阵，然而在 NCMVN 下这一结论不成立；
lowertruncation(numlist)：指明截断下限向量，用 . 代表 $- \infty$ ；
uppertruncation(numlist)：指明截断上限向量，用 . 代表 $+ \infty$ 。

4. 典例分析

4.1 密度比较

为解释 mvnormalden、mvtden、tmvnormalden、tmvden 等命令的用法，我们绘制并比较在给定下列参数情形下的不同密度。

δ = (\begin{array}{l} 0 \\ 0 \end{array}), Σ = (\begin{array}{cc} 1 & 0.5 \\ 0.5 & 1 \end{array}), ν = 1, l = (\begin{array}{l} - 1.5 \\ - 1.5 \end{array}), u = (\begin{array}{l} 1.5 \\ 1.5 \end{array})

首先在区域 $[- 3, 3]^{2}$ 循环一系列值，在每一个网格点上我们计算一个 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)

由上图可知， $ν = 1$ 阶自由度的 MVT 分布和 MVN 分布更相近，这也解释了低自由度 MVT 分布的尾部概率要更大。与非截断相比，截断会增加在区域 $[- 1.5, 1.5]^{2}$ 内每一点的密度。

4.2 多元正态分布的拟合优度检验

我们考虑二元正态分布，其均值向量和协方差矩阵为：

δ = (\begin{array}{l} 0 \\ 0 \end{array}), Σ = (\begin{array}{cc} 1 & 0.5 \\ 0.5 & 1 \end{array})

我们同时也考虑两个截断 MVT 分布的混合，每一个MVT 分布都具有上述的协方差矩阵 $Σ$ ，但是其中一个具有如下的均值和截断条件：

δ = (\begin{array}{l} - 2 \\ - 2 \end{array}), l = (\begin{array}{l} - 2.5 \\ - 2.5 \end{array}), u = (\begin{array}{l} - 1.5 \\ - 1.5 \end{array})

另一个具有如下均值和截断条件：

δ = (\begin{array}{l} 2 \\ 2 \end{array}), l = (\begin{array}{l} 1.5 \\ 1.5 \end{array}), u = (\begin{array}{l} 2.5 \\ 2.5 \end{array})

我们额外考虑一些可能的 $ν$ 值。然后对于不同的 $n$ 值，我们从 MVT 分布中抽取 $n$ 个样本，这一抽取重复 1000 次，之后计算出 Doornik-Hansen 检验中的 $p$ 值平均数。这一过程同样对于截断的 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 分解。对于绝大多数分布而言，分解方法的选择对于生成随机变量的速度和质量影响不大，但是当维度 $k$ 很大时，分解方法的选择就需要非常谨慎了。

从上述结果可以看出，MVT 分布中 $ν$ 的增加会增加 $p$ 值的平均值，这一点对于所有的 $n$ 都成立。这是因为随着 $ν$ 的增加，MVT 分布逐渐接近相应的 MVN 分布。而对于截断 NCMVT 变量的混合，平均的 $p$ 值总是高于相应的 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

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Stata：多元正态t截断分布的命令

1. 问题背景

2. 原理

2.1 多元正态分布

2.2 多元 t 分布

2.3 截断 MVN 和 MVT 分布

3. 命令介绍

4. 典例分析

4.1 密度比较

4.2 多元正态分布的拟合优度检验

5. 参考文献

6. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯