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

发布时间:2022-03-12 阅读 376

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 系列推文:

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,且独立同分布的随机变量,其任何一组随机样本 X1,...,Xn 在 n 时,n(i=1nXi/nμ) 依分布收敛于 N(0,σ2)

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

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

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

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

2. 原理

2.1 多元正态分布

考虑 k 维具有非退化的 MVN 分布的随机变量 X=(X1,,Xk),我们使用 XNk(δ,Σ) 来表示。其中 δRk 是一个位置参数,ΣMk×k 是一个正定的协方差矩阵。

以下的计算都是围绕着 Φk(a,b,δ,Σ) 展开。其中 a=(a1,,ak),以及 b=(b1,,bk)。对于 i=1,,k,$-\infty \leq a_{i}

其中,|Σ|=det(Σ)

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

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

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

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

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

最后,为得到分布 Nk(δ,Σ) 的伪随机向量,我们将用到如下事实:

其中 z=(z1,,zk)ziN(0,1) 可以通过 rnormal() 来实现。我们将使用 rmvnormal 对这一伪随机分布进行估计。

2.2 多元 t 分布

考虑一个 k 维具有 NCMVT 分布的随机变量 X=(X1,,Xk),我们用 XTk(δ,Σ,ν) 来代表。其中 δRk 是一个非中心参数向量,ΣMk×k 是一个正定矩阵,νR+ 是自由度参数。

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

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

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

2.3 截断 MVN 和 MVT 分布

我们进一步考虑截断情形下的多元正态分布和非中心化的多元 t 分布。我们分别使用 XTNk(δ,Σ,l,u) 和 XTTk(δ,Σ,ν,l,u) 来表示。其中 l=(l1,,lk)Rku=(u1,,uk)Rk,并且对于 i=1,...,k,$l_{i}

其中 ϕk(x,δ,Σ,l,u) 和 ψk(x,δ,Σ,ν,l,u) 是其各自的密度。我们可以使用 tmvnormaltmvttmvnormaldentmvtden 命令对上述分布和密度进行计算。

进一步使用 invmvnormalinvtmvt 命令确定 q 的值,使得下式成立:

或者,

最后在命令 rtmvnormalrtmvt 中使用排斥抽样法 (rejection sampling) 产生具有 TNk(δ,Σ,l,u) 和 TTk(δ,Σ,ν,l,u) 分布的随机变量。简单来说,我们是从 Nk(δ,Σ) 进行随机抽样。如果抽取的样本 Xi 对于 i=1,...,k 满足 liXiui,那么我们就保留这次抽样结果,否则我们就拒绝这次结果,重复上述过程。

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 密度比较

为解释 mvnormaldenmvtdentmvnormaldentmvden 等命令的用法,我们绘制并比较在给定下列参数情形下的不同密度。

首先在区域 [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 多元正态分布的拟合优度检验

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

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

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

我们额外考虑一些可能的 ν 值。然后对于不同的 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小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh