This is just an extended comment showing that the same thing happens with the Probability
function:
pd = ProductDistribution[{NormalDistribution[0, 10], 4}]
Probability[x1 > 0, {x1, x2, x3, x4} \[Distributed] pd]
(* 1/2 *)
Probability[x1 < 0 && x1 + x2 > 0, {x1, x2, x3, x4} \[Distributed] pd]
(* 1/8 *)
FullSimplify[Probability[x1 < 0 && x1 + x2 < 0 && x1 + x2 + x3 > 0,
{x1, x2, x3, x4} \[Distributed] pd]]
(* 1/16 *)
Probability[x1 < 0 && x1 + x2 < 0 && x1 + x2 + x3 < 0 && x1 + x2 + x3 + x4 > 0,
{x1, x2, x3, x4} \[Distributed] pd]
(* Probability[x1<0&&x1+x2<0&&x1+x2+x3<0&&x1+x2+x3+x4>0,
{x1,x2,x3,x4}\[Distributed]ProductDistribution[{NormalDistribution[0,10],4}]] *)
Also any positive standard deviation gives the same set of probabilities.
Addition:
I suspect the answer for f@4
is 5/128. Here's why using a brute force approach. We constuct a new set of variables with y1 = x1
, y2 = x1 + x2
, y3 =
x1 + x2 + x3
, and y4 = x1 + x2 + x3 + x4
.
a = {{1, 0, 0, 0}, {1, 1, 0, 0}, {1, 1, 1, 0}, {1, 1, 1, 1}};
a.{x1, x2, x3, x4}
(* {x1, x1 + x2, x1 + x2 + x3, x1 + x2 + x3 + x4} *)
(Σ = a.DiagonalMatrix[{1, 1, 1, 1}].Transpose[a]);
d = MultinormalDistribution[{0, 0, 0, 0}, Σ];
Integrate[PDF[d, {y1, y2, y3, y4}], {y1, -∞, 0}, {y2, -∞, 0}, {y3, -∞, 0}, {y4, 0, ∞}]
$$\frac{\int _{-\infty }^0\int _{-\infty }^0\left(\text{erf}\left(\frac{\text{y3}}{\sqrt{2}}\right)+1\right) e^{\frac{1}{4} \left(-3 \text{y1}^2+2 \text{y1} \text{y3}-\text{y3}^2\right)} \text{erfc}\left(\frac{\text{y1}+\text{y3}}{2}\right)d\text{y3}d\text{y1}}{8 \sqrt{2} \pi }$$
NIntegrate[PDF[d, {y1, y2, y3, y4}], {y1, -∞, 0}, {y2, -∞, 0}, {y3, -∞, 0}, {y4, 0, ∞}]
(* 0.03906249965228397 *)
5/128 // N
(* 0.0390625 *)
Maybe using simulations is the only possible approach for values larger than 4.
Simulation approach
Here are the results using simulations.
a = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 0, 0, 0, 0, 0,
0}, {1, 1, 1, 0, 0, 0, 0, 0, 0, 0},
{1, 1, 1, 1, 0, 0, 0, 0, 0, 0}, {1, 1, 1, 1, 1, 0, 0, 0, 0, 0}, {1,
1, 1, 1, 1, 1, 0, 0, 0, 0},
{1, 1, 1, 1, 1, 1, 1, 0, 0, 0}, {1, 1, 1, 1, 1, 1, 1, 1, 0, 0}, {1,
1, 1, 1, 1, 1, 1, 1, 1, 0},
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};
(Σ =
a.DiagonalMatrix[{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}].Transpose[a]);
d = MultinormalDistribution[{0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, Σ];
nsim = 10000000;
SeedRandom[12345];
x = RandomVariate[d, nsim];
prob = {Length[Select[x, #[[1]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] < 0 && #[[8]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] < 0 && #[[8]] < 0 && #[[9]] > 0 &]]/nsim // N,
Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] < 0 && #[[8]] < 0 && #[[9]] < 0 && #[[10]] > 0 &]]/nsim // N}
{0.499966, 0.125064, 0.0625243, 0.0390055, 0.0273741, 0.0204875,
0.0161639, 0.0130903, 0.0109223, 0.0092713}