背景

背景：在我的一篇文章中，我指出，如果一个选拔过程（如高等教育）要求成功的申请者在 2 个不同的变量上高于平均值，那么成功的申请者必然会比成功的申请者必须高于平均数时更少仅 1 个变量的平均值。（本例中的 2 个不同的变量是尽责性和智商的大五人格因素。）当 2 个变量的相关性为 0 时，这种减少最为显着，但即使相关性很大，仍然会导致许多申请人被过滤掉。到底有多少被过滤掉了？

简单的问题

好吧，如果过滤器适用于均值以上 2 个标准差，并且变量与 1 相关，那么 2.3% 的总体将通过；如果变量与 0 不相关，则 2.3% * 2.3% 或 5.29e-4% 的人口将通过。

相关

但是中间值呢？例如，心理学文献报道了责任心和智商之间的相关性为-0.21。

我在双变量正态分布上查阅了维基百科，但我不太了解。我发现的最接近的是相关正态随机变量的总和，但在这种情况下，我想要的是更接近min函数。

模拟

我能够进行一个 R 模拟来看看它是如何工作的，它似乎符合我的直觉：

install.packages("fMultivar")
library ("fMultivar")

x <- rnorm2d(10000000, rho=0.5)

xgreater <- length(subset(x, x[,1] > mean(x[,1])+2*sd(x[,1])))
xandygreater <- length(subset(x, x[,1] > mean(x[,1])+2*sd(x[,1]) & x[,2] > mean(x[,2])+2*sd(x[,2])))

c(xgreater, xandygreater); c(xgreater / length(x), xandygreater / length(x), xgreater / xandygreater) * 100

# example results for different values of 'rho='
0.1
[1] 454,664  17,570
[1] 2.273e+00 8.785e-02 2.588e+03

0.2
[1] 458,284  82,552
[1]   2.2914   0.4128 555.1458

0.5
[1] 454,484  80,872
[1]   2.2724   0.4044 561.9794

0.9
[1] 455,242 267,912
[1]   2.276   1.340 169.922

0.95
[1] 455,162 321,024
[1]   2.276   1.605 141.784

0.99
[1] 455,260 394,448
[1]   2.276   1.972 115.417

精确的pdf计算？

我真的希望有更多精确的分析解决方案，所以更多的搜索最终出现了一篇论文，“两个高斯随机变量的最大值/最小值的精确分布”，它给出了 2 个相关正态变量的最小值的定义。这似乎是我想要的；pg1 的顶部，第二列：

...其中和分别是标准正态分布的 pdf 和累积分布函数 (cdf)。已知的pdf为，其中 $\phi(.)$ $\Phi(.)$ $Y = \min(X_1, X_2)$ $f(y) = f_1(y) + f_2(y)$

(3) (4) $f_1(y) = \frac{1}{\sigma_1} \phi (\frac{y-\mu_1}{\sigma_1}) \times \Phi (\frac{p(y - \mu_1)}{\sigma_1 \sqrt{1 - p^2}} - \frac{y - \mu_2}{\sigma_2 \sqrt{1 - p^2}})$ $f_2(y) = \frac{1}{\sigma_2} \phi (\frac{y-\mu_2}{\sigma_2}) \times \Phi (\frac{p(y - \mu_2)}{\sigma_2 \sqrt{1 - p^2}} - \frac{y - \mu_1}{\sigma_1 \sqrt{1 - p^2}})$

他们在 pg6 上给出了一个 R 实现（第一列）；它似乎有一个pnorm错字，但我修复了它。一旦它开始工作，我尝试生成一个稍微 (0.1) 相关的双变量分布，看起来不错：

fmin<-function (y,mu1,mu2,sigma1,sigma2,rho)
     {t1<-dnorm(y,mean=mu1,sd=sigma1)
     tt<-rho*(y-mu1)/(sigma1*sqrt(1-rho*rho))
     tt<-tt-(y-mu2)/(sigma2*sqrt(1-rho*rho))
     t1<-t1*pnorm(tt)
     t2<-dnorm(y,mean=mu2,sd=sigma2)
     tt<-rho*(y-mu2)/(sigma2*sqrt(1-rho*rho))
     tt<-tt-(y-mu1)/(sigma1*sqrt(1-rho*rho))
     t2<-t2*pnorm(tt)
     return(t1+t2)}
fmin(c(1:200),100,100,15,15,0.1)
  [1] 1.849e-11 2.864e-11 4.418e-11 6.784e-11 1.037e-10 1.578e-10 2.392e-10 3.608e-10 5.418e-10
 [10] 8.101e-10 1.206e-09 1.787e-09 2.636e-09 3.872e-09 5.663e-09 8.243e-09 1.195e-08 1.724e-08
 [19] 2.476e-08 3.542e-08 5.043e-08 7.148e-08 1.009e-07 1.417e-07 1.982e-07 2.760e-07 3.827e-07
 [28] 5.282e-07 7.257e-07 9.928e-07 1.352e-06 1.833e-06 2.475e-06 3.326e-06 4.449e-06 5.926e-06
 [37] 7.859e-06 1.037e-05 1.364e-05 1.784e-05 2.324e-05 3.014e-05 3.891e-05 5.002e-05 6.401e-05
 [46] 8.154e-05 1.034e-04 1.306e-04 1.641e-04 2.054e-04 2.558e-04 3.173e-04 3.917e-04 4.814e-04
 [55] 5.889e-04 7.173e-04 8.696e-04 1.049e-03 1.261e-03 1.507e-03 1.794e-03 2.125e-03 2.506e-03
 [64] 2.941e-03 3.435e-03 3.993e-03 4.620e-03 5.318e-03 6.093e-03 6.945e-03 7.878e-03 8.890e-03
 [73] 9.982e-03 1.115e-02 1.239e-02 1.370e-02 1.506e-02 1.647e-02 1.791e-02 1.938e-02 2.084e-02
 [82] 2.230e-02 2.371e-02 2.508e-02 2.636e-02 2.755e-02 2.863e-02 2.956e-02 3.034e-02 3.095e-02
 [91] 3.138e-02 3.162e-02 3.165e-02 3.149e-02 3.112e-02 3.056e-02 2.981e-02 2.889e-02 2.781e-02
[100] 2.660e-02 2.526e-02 2.383e-02 2.233e-02 2.077e-02 1.920e-02 1.762e-02 1.605e-02 1.452e-02
[109] 1.305e-02 1.164e-02 1.031e-02 9.063e-03 7.912e-03 6.857e-03 5.899e-03 5.039e-03 4.272e-03
[118] 3.595e-03 3.004e-03 2.491e-03 2.050e-03 1.675e-03 1.358e-03 1.093e-03 8.732e-04 6.923e-04
[127] 5.447e-04 4.254e-04 3.297e-04 2.535e-04 1.935e-04 1.466e-04 1.102e-04 8.220e-05 6.085e-05
[136] 4.470e-05 3.258e-05 2.357e-05 1.692e-05 1.205e-05 8.517e-06 5.973e-06 4.157e-06 2.870e-06
[145] 1.966e-06 1.337e-06 9.018e-07 6.036e-07 4.008e-07 2.641e-07 1.727e-07 1.120e-07 7.210e-08
[154] 4.604e-08 2.917e-08 1.834e-08 1.144e-08 7.078e-09 4.346e-09 2.647e-09 1.600e-09 9.594e-10
[163] 5.708e-10 3.369e-10 1.973e-10 1.146e-10 6.607e-11 3.778e-11 2.143e-11 1.207e-11 6.737e-12
[172] 3.733e-12 2.052e-12 1.119e-12 6.052e-13 3.248e-13 1.730e-13 9.137e-14 4.788e-14 2.489e-14
[181] 1.284e-14 6.571e-15 3.336e-15 1.680e-15 8.393e-16 4.160e-16 2.046e-16 9.979e-17 4.830e-17
[190] 2.319e-17 1.104e-17 5.218e-18 2.446e-18 1.137e-18 5.248e-19 2.402e-19 1.090e-19 4.911e-20
[199] 2.194e-20 9.727e-21

绘制pdf

现在，我将 PDF 理解为“描述该随机变量取给定值的相对可能性的函数。随机变量落入特定区域的概率由该变量的密度在地区”。所以我想我应该总结 pdf > 130 中的每个点（因为 130 是 2 个标准偏差，当我指定 SD=15 时按构造），这就是我随机偏差的概率min(130,130)。某人在两个变量上都超过 130 的总概率是多少？我认为那将是：

sum(fmin(c(1:200),100,100,15,15,0.1)[130:200])
[1] 0.001004

如果我将r增加到 0.9，结果是 0.01455，这是令人满意的大。

健全性检查 - 当相关性达到 1.0 时，应该不会减少。因此，我们对以相同方式定义的单个正态分布执行相同的问题：

sum(dnorm(c(1:200), 100, 15)[130:200])
[1] 0.02459

# the function blows NaN chunks on 1.0, so we'll try a lot of 9s:
sum(fmin(c(1:200),100,100,15,15,0.9999999999)[130:200])
[1] 0.02459

我想这也有效。

问题

所以我的问题是：

我的 R 模拟正确吗？
我的版本和使用fmin对吗？
是否有一些更直接的，甚至可能是纸笔计算的方法来计算我最初问题的答案？

3个回答

如果我正确理解了您的问题，您需要二元正态分布的 CDF。也就是说，对于标准化的情况：

Φ (b b, ρ) = \frac{1}{2 π \sqrt{1 - ρ^{2}}} \int_{- \infty}^{b_{1}} \int_{- \infty}^{b_{2}} \exp [- (x^{2} - 2 ρ x y + y^{2}) / (2 (1 - ρ^{2})] d y d x

$\Phi(\mathrm{\pmb{b}},\rho) = \frac{1}{2\pi\sqrt{1-\rho^{2}}}\int_{-\infty}^{b_{1}}{\int_{-\infty}^{b_{2}}}\exp\left[{-(x^{2}-2\rho xy+y^{2}})/(2(1-\rho^{2})\right]\mathrm{d}y\mathrm{d}x$

这没有封闭形式的解决方案，必须进行数值积分。使用现代软件，这是非常微不足道的。

这是 R 中具有完全相关正态分布（即）的示例，其中和协方差矩阵。我们计算两个变量高于 SD 的概率： $\rho = 1$ $\pmb{\mu}=(100, 40)^\intercal$ $\Sigma = \begin{bmatrix} 225 & 75 \\ 75 & 25 \end{bmatrix}$ $2$

library(mvtnorm)
library(MASS)

corr.mat <- matrix(c(1, 1, 1, 1), 2, 2, byrow = TRUE) # correlations
sd.mat <- matrix(c(15, 0, 0, 5), 2, 2, byrow = TRUE) # standard deviations

cov.mat <- sd.mat %*% corr.mat %*% sd.mat # covariance matrix

mu <- c(100, 40) # means

pmvnorm(lower = mu + 2*diag(sd.mat), upper = Inf, mean = mu, sigma = cov.mat)

[1] 0.02275013
attr(,"error")
[1] 2e-16
attr(,"msg")
[1] "Normal Completion"

正如你所说的：当它们完全相关时，概率约为 2.3%。

-0.21 的相关性怎么样？

corr.mat <- matrix(c(1, -0.21, -0.21, 1), 2, 2, byrow = TRUE) # correlations
sd.mat <- matrix(c(15, 0, 0, 5), 2, 2, byrow = TRUE) # standard deviations

cov.mat <- sd.mat %*% corr.mat %*% sd.mat # covariance matrix

mu <- c(100, 40) # means

pmvnorm(lower = mu + 2*diag(sd.mat), upper = Inf, mean = mu, sigma = cov.mat)

[1] 0.0001228342
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

概率要低得多，即 0.0001228342。

我们可以通过模拟来验证我们的计算。对于上面的例子：

set.seed(21)
dat <- rmvnorm(1e7, mean = c(100, 40), sigma = cov.mat)

sum(dat[, 1] > (100+2*15) & dat[, 2] > (40+2*5))/dim(dat)[1]

[1] 0.0001261

这与数值积分的结果非常接近。

这些计算可以很容易地扩展到多元正态分布。

好问题。你的直觉和方法是正确的。但是，对于您所追求的联合概率，没有简单的分析公式。可以用 R 编写一个简短的程序来计算概率，但它必须使用数值积分。

这是#3的有用解释：多元正态分布

请参阅该部分中的“二元正态分布”，您可以看到具有相关系数的二元正态分布的 pdf。您可以将-infinity与 + 2*之间的pdf集成。或者也许这是足够的信息来找到CDF的方程，可以在 + 2*处进行评估，从而为您提供所需的东西。 $\mu$ $\sigma$ $\mu$ $\sigma$

其它你可能感兴趣的问题

上一篇具有非恒定采样间隔的时间序列预测下一篇非线性回归的指数（在 R 中）？

相关双变量正态分布：找到比平均值高 2 个标准差的数据百分比？

背景

简单的问题

相关

模拟

精确的pdf计算？

问题