多元概率单位模型

多元概率单位模型（multivariate probit model）用相关的潜在正态变量描述同一对象的多个二元响应。例如同一家企业可能同时选择两项服务，每项都是“采用或不采用”，两项选择又可能受到未观测共同因素的影响。模型分别描述各项采用概率，并用潜在相关结构描述联合选择。

这里“多元”指多个响应变量。一个二元响应含多个自变量的多变量概率单位回归，仍可是一项普通二元模型；多个互斥类别的多项概率单位模型又有不同的潜在效用与识别约束。若企业可以两项都采用，不能把这四种组合误当成两项互斥选择。

潜在方程与尺度约束

对对象 $i$ 的第 $j$ 个响应，设 $$\begin{array}{cc}{Z}_{ij}& ={\eta }_{ij}+{\epsilon }_{ij},{\eta }_{ij}=x_{ij}^{𝖳}{\beta }_j,\\ Y_{ij}& =\mathrm{𝟏}\{Z_{ij}>0\},j=1,\dots ,q.\end{array}$$ 给定解释变量，不同对象的误差向量独立，每个对象内部满足 $${\epsilon }_i\sim N_q(0,R),R_{jj}=1.$$ 本文采用正定相关矩阵 $R$，各方程可以使用不同的解释变量。单位对角线是识别所需的一项尺度规范：把某个潜在方程及其误差同时乘任意正常数，并不改变 $Z_{ij}>0$，所以只有二元输出时，不能再自由估计每个潜在边际的任意尺度。Cappellari与Jenkins的作者讲义给出这一多方程设定。

每个边际仍为普通概率单位形式 $Pr(Y_{ij}=1\mid x)=\Phi ({\eta }_{ij})$，但联合概率一般不是边际概率相乘。若 $q>2$，单独保证每对相关系数处于 $(-1,1)$ 也不够；它们必须共同构成正定矩阵。

联合似然的符号变换

对观察向量 $y_i\in \{0,1{\}}^q$，令 $s_{ij}=2y_{ij}-1$，$D_i=\mathrm{diag}(s_{i1},\dots ,s_{iq})$。观察到该向量的事件可写成 $D_i{Z}_i>0$，即各坐标都为正。由于 $-D_i{\epsilon }_i$ 仍为零均值正态，其协方差为 $D_{i}R{D}_i$，所以 $$Pr(Y_i=y_i\mid x_i)={\Phi }_q(D_i{\eta }_i;D_{i}R{D}_i).$$ 这里 ${\Phi }_q(a;C)=Pr(W_1\le a_1,\dots ,W_q\le a_q)$，$W\sim N_q(0,C)$。连续正态分布的边界概率为零，因此不等号的严格与否不改变这个概率。

将这些联合概率对独立对象相乘，才得到样本似然。二元情形的有符号参数形式也可见Stata官方手册：相关系数应随两个响应的符号一同变换。

零阈值二元情形的扇区证明

取两个潜在均值都为零，相关系数为 $\rho \in (-1,1)$。写 $$Z_1=U,Z_2=\rho U+\sqrt{1-{\rho }^2}V,$$ 其中 $U,V$ 独立标准正态。令 $\varphi =\arccos \rho \in (0,\pi )$，则 $Z_2=\cos \varphi U+\sin \varphi V$。

二维标准正态密度只依赖半径，极坐标的角度因此在整圈上均匀。事件 $U>0$ 对应角度位于一个半圆，事件 $\cos \varphi U+\sin \varphi V>0$ 对应中心方向转过 $\varphi$ 的另一个半圆。两个正半圆的重叠角宽为 $\pi -\varphi$。因此 $$\begin{array}{cc}Pr(Y_1=1,Y_2=1)& =\frac{\pi -\varphi }{2\pi }\\ & =\frac{1}{4}+\frac{\arcsin \rho }{2\pi }.\end{array}$$ 正负中心对称给出相同的 $p_{00}$；两个边际都为 $1/2$，故 $$p_{10}=p_{01}=\frac{1}{4}-\frac{\arcsin \rho }{2\pi }.$$ 这项推导只适用于这里的零阈值二元情形；有协变量造成非零均值后，一般需要计算相应的二维正态分布函数。

取潜在相关系数0.5，两个正半平面的重叠扇区角宽为120度，占整圈三分之一。所画圆只用来显示角度；概率包括该扇区内所有半径，不能把圆内面积直接当作正态概率。

当 $\rho =1/2$ 时，$\arcsin \rho =\pi /6$，从而 $p_{11}=p_{00}=1/3$、$p_{10}=p_{01}=1/6$。两个二元响应的可观测相关系数是 $$\begin{array}{cc}\mathrm{Corr}(Y_1,Y_2)& =\frac{p_{11}-1/4}{1/4}\\ & =\frac{2}{\pi }\arcsin \rho =\frac{1}{3}.\end{array}$$ 它并不等于潜在相关系数 $1/2$。在非零阈值下，两者关系还取决于各自边际概率。

潜在均值为零且相关系数为0.5时，同为零或同为一比一零或零一更常见；两项边际采用率仍各为二分之一。

三十六个对象的受限拟合

假设潜在均值预先固定为零，只估计共同相关系数；三十六个独立对象的四格计数为 $n_{11}=n_{00}=12$、$n_{10}=n_{01}=6$。记 $t=p_{11}=p_{00}$，则 $p_{10}=p_{01}=1/2-t$，$0<t<1/2$，对数似然为 $$\ell (t)=24\log t+12\log (1/2-t).$$ 一阶导数为 $24/t-12/(1/2-t)$，令其为零，得到 $t=1/3$。二阶导数严格为负，且端点使似然趋于零，所以这就是唯一极大值。反解概率公式： $$\hat{\rho }=\sin \{2\pi (t-1/4)\}=\sin (\pi /6)=\frac{1}{2}.$$ 这是一个只有一项相关参数的可复算例子；不能据此省略一般模型中回归系数、相关矩阵约束与多维积分的联合拟合。

变式。 若 $\rho =0$，四格概率均为 $1/4$，在这个联合正态模型内两项选择条件独立。若 $\rho =-1/2$，同号两格各为 $1/6$，异号两格各为 $1/3$，可观测相关系数为 $-1/3$。

多个二元响应的相关可以反映共同未观测因素，也可能来自遗漏变量或模型设定；单靠相关系数不能推出某项选择导致另一项选择。高维正态概率的数值积分或模拟还带有计算误差，应检查近似精度、随机种子和优化稳定性。尺度规范、可识别的设计矩阵、正定相关结构与独立观测单位，是解释估计结果前必须明确的条件。