Which orthant probabilities are the largest? (For a multivariate normal distribution)

Question

I have a $k$-dimensional multivariate normal distribution $X∼N(0,\Sigma)$ with covariance matrix $\Sigma$. $\Sigma$ has two distinct eigenvalues, say $\lambda_1 > \lambda_2$, with orthogonal eigenspaces $V_1$ and $V_2$. I am interested in the orthant probabilities; given an orthant defined by $\epsilon = (\epsilon_1,\ldots,\epsilon_k)$, with each $\epsilon_i \in \{1,-1\}$, there is an orthant probability $p_\epsilon=Pr(\forall i: \;\epsilon_i X_i \geq 0)$.

There seems to be a literature on finding closed forms for these in special cases, but I do not necessarily need a closed form, but rather I just want to know when $p_{\epsilon} \geq p_{\epsilon'}$. One can write $\epsilon$ as a sum $\epsilon = v_1 + v_2$ with $v_1 \in V_1$ and $v_2 \in V_2$, and $k = \lVert\epsilon\rVert^2 = \lVert v_1 \rVert^2+\lVert v_2\rVert^2$. My conjecture is that the greater $\lVert v_1\rVert$ is, the greater $p_\epsilon$ is (recall that $\lambda_1 > \lambda_2$). Geometrically, this means that the diagonal vector of the octant defined by $\epsilon$ is closer to the longer axes of the ellipsoids which are the equidensity contour of the distribution.

Is this true? And if so, how does one prove it? I have been stumped (though I know little about probability).

(In the particular case I am interested in, I have a complete graph $K_\ell$, and I am generating a random map from edges of $K_\ell$ to $\mathbb{R}$. The covariance matrix has, as its two eigenspaces, the cut space and the edge space of $K_\ell$.)

Answer 1

Surprisingly, the conjecture is false. Orthant probabilities are not always ordered by the vector norm $\lVert v_1 \rVert$.

For a counterexample, take the 4-dimensional normal distribution with the following covariance matrix:

  75.990348312987877 -14.891382893378880   9.019694111954474   8.518647696984502
 -14.891382893378880  21.973634189398940  31.423796140503811  -5.365552691385230
   9.019694111954474  31.423796140503811  63.319136031398756  -5.748592963390807
   8.518647696984502  -5.365552691385230  -5.748592963390807   2.716881466214341

The two distinct eigenvalues are $\lambda_1 = 81$ and $\lambda_2 = 1$, each one twice. (It is not a coincidence that they are integers: I created the example by stretching a standard normal distribution $9$-fold in two directions, and rotating randomly. If you re-calculate the eigenvalues from the matrix listed above, they match to more than 12 decimal places.)

Looking into the orthants (only showing first half due to symmetry):

orthant  0: ++++  |v1|=1.604343  p=0.081526   p_MC=0.081561
orthant  1: +++-  |v1|=1.644290  p=0.103843   p_MC=0.103680
orthant  2: ++-+  |v1|=0.916062  p=0.003139   p_MC=0.003152
orthant  3: ++--  |v1|=0.627763  p=0.002054   p_MC=0.002070
orthant  4: +-++  |v1|=1.419708  p=0.068214   p_MC=0.068349  !!!
orthant  5: +-+-  |v1|=1.268382  p=0.017151   p_MC=0.017054
orthant  6: +--+  |v1|=1.850191  p=0.198087   p_MC=0.198265
orthant  7: +---  |v1|=1.562552  p=0.025999   p_MC=0.025957  !!!

Comparing orthants 4 and 7, we observe the latter has greater $\lVert v_1 \rVert$, but clearly smaller probability. I calculated the probabilities with two methods: p is from MATLAB mvncdf which claims absolute error tolerance $10^{-4}$, and p_MC is from Monte-Carlo integration with $10^7$ points. The probabilities differ already in the second decimal, so I'm pretty confident that it is not just an artefact of numerical calculation.

More counterexamples can be generated by making random instances and checking the orthant probabilities: the conjecture fails every now and then, although it holds in most of the cases. Interestingly, in 3 dimensions I have not found any counterexample in thousands of random instances (either in the prolate case = large eigenvalue is single and small eigenvalue is double, or in the oblate case = large is double and small is single).

Figure

Here is an attempt to visualize the 4-dimensional distribution (if anyone has better suggestions, I'd like to hear). We are showing $100\,000$ random points from the distribution; points in the 4th orthant +-++ are shown red, points in the 7th orthant +--- are shown blue, and all other points cyan. Perhaps one can see that the red points are more than twice as many as the blue points.

Each subplot shows a 2D projection to two coordinate axes, in the same scale: each box ranges from $-20$ to $+20$ in both directions. Note that the two contending orthants are in the same quadrant of $X_1,X_2$, and in opposite quadrants of $X_3,X_4$.

Code

For convenience, here is Matlab code to reproduce and verify this particular instance. You can adapt it to search for other counterexamples in four or other dimensions.

Sigma = [
  75.990348312987877 -14.891382893378880   9.019694111954474   8.518647696984502
 -14.891382893378880  21.973634189398940  31.423796140503811  -5.365552691385230
   9.019694111954474  31.423796140503811  63.319136031398756  -5.748592963390807
   8.518647696984502  -5.365552691385230  -5.748592963390807   2.716881466214341
   ];

dim   = size(Sigma,1);
[V,D] = eigs(Sigma);
Dval  = diag(D);

% Find the big and small eigenvalues
Ibig   = abs(Dval - max(Dval)) < 1e-3;
Ismall = abs(Dval - min(Dval)) < 1e-3;
fprintf('Big eigenvalues:\n');
Dval(Ibig)
fprintf('Small eigenvalues:\n');
Dval(Ismall)

% Generate random points from the distribution
NN = 1e7;
R  = chol(Sigma);
X = R' * randn(dim,NN);

V1MAG = [];
P     = [];
for i=0:2^(dim-1)-1
    neg  = bitget(i, dim:-1:1)';
    e = ones(dim,1) - 2*neg;
    signvec = repmat('+',1,dim);
    signvec(e<0) = '-';
    
    % Orthant probability from Matlab
    Sigmaflip = e .* Sigma .* e';
    p = mvncdf(zeros(dim,1), zeros(dim,1), Sigmaflip);
    P = [P p];
    
    % Monte Carlo orthant probability
    pmc = mean(all(X .* e > 0, 1));
    
    % Projections of e onto each eigenvector (column of V)
    edots = e' * V;
    
    v1mag = sqrt(sum(edots(Ibig).^2));
    v2mag = sqrt(sum(edots(Ismall).^2));
    V1MAG = [V1MAG v1mag];
    
    fprintf('orthant %2d: %s  |v1|=%.6f  p=%.6f   p_MC=%.6f\n', ...
        i, signvec, v1mag, p, pmc);
end

% Check the order
[~,S] = sort(V1MAG);
pprev = -inf;
fprintf('\nOrthants sorted by |v1|:\n');
for i=1:2^(dim-1)
    fprintf('orthant %2d : |v1|=%.6f  p=%.6f', ...
        S(i)-1, V1MAG(S(i)), P(S(i)));
    if P(S(i)) < pprev - 0.01
        fprintf('   TROUBLE!');
    end
    fprintf('\n');
    pprev = P(S(i));
end