# Speed up Probability

I have a set of inequalities:

(β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y < 1.}}, 0.] >
(β - z - 23/ 2) Piecewise[{{1., -1. < w <  0.},
{0.3333333333333333 (3. - 1. w), 0. <= w < 1.}},  0.] && (β - x -19/2)
Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y <  1.}}, 0.] >
(β - z - 21/2) Piecewise[{{0.3333333333333333 (2. - 1. w), -1. < w < 1.}}, 0.]


and I'm trying to calculate the probability that these inequalities hold given that β is distributed as a uniform between 8.5 and 11.5 and that 1 > z >= w > -1 && 1 > x >= y > -1.

I'm doing this:

Probability[(β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y),
-1. < y <  1.}}, 0.] > (β - y - 23/ 2)
Piecewise[{{1., -1. < w <  0.}, {0.3333333333333333 (3. - 1. w), 0. <= w < 1.}},  0.] &&
(β - z - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y <  1.}}, 0.] >
(β - z - 21/2) Piecewise[{{0.3333333333333333 (2. - 1. w), -1. < w < 1.}}, 0.] &&
1 > z >= w > -1 && 1 > x >= y > -1,
β \[Distributed] UniformDistribution[{8.5, 11.5}]]


This has been running for 26 hours now and mathematica has not found a solution yet. Is there any way I can speed up the calculation?

If the parameter space is just x and y it takes less than a second to get the solution:

Probability[(β - x - 19/2) Piecewise[{{0.3333333333333333 (1. -1. y), -1.    < y <  1.}},
0.] > (β - x - 23/2) Piecewise[{{1.,-1.< y <  0.}, {0.3333333333333333 (3. - 1. y),
0. <= y < 1.}}, 0.] && (β - x -  19/2) Piecewise[{{0.3333333333333333 (1. - 1. y),
-1. < y <  1.}}, 0.] >
(β - x -21/2) Piecewise[{{0.3333333333333333 (2. -1. x), -1. < x < 1.}}, 0.] &&
1 > x > -1 && 1 > y > -1, β \[Distributed]  UniformDistribution[{8.5, 11.5}]]


if I try to use reduce with these two inequalities:

kk = Table[Reduce[(β - x - 19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}}, 0] >
(β - y - 23/2) Piecewise[{{1, -1 < w < 0}, {1/3 (3 - w), 0 <= w < 1}}, 0] &&
(β - x - 19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}}, 0] >
(β - y - 21/2) Piecewise[{{1/3 (2 - w), -1 < w < 1}}, 0] && 1 > z >= w > -1
&& 1 > x >= y > -1 /. β -> a, {x, y, z, w}, Reals], {a, 17/2, 23/2, 1/10}]


I get the solution in about 5 minutes.

If instead I try to reduce these inequalities:

kk = Table[Reduce[(β - x - 19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}}, 0] >
(β - z - 23/2) Piecewise[{{1, -1 < w < 0}, {1/3 (3 - w), 0 <= w < 1}}, 0] &&
(β - x - 19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}}, 0] >
(β - z - 21/2) Piecewise[{{1/3 (2 - w), -1 < w < 1}}, 0] && 1 > z >= w > -1
&& 1 > x >= y > -1 /. β -> a, {x, y, z, w}, Reals], {a, 17/2, 23/2, 1/10}]


I can't get the solution within 2 hours. Note that in the first case the parameters are effectively 3 since z is not used in any inequality, whereas in the second case the parameters are 4.

Edited: the first set of inequalities I posted was

(β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y < 1.}}, 0.] >
(β - y - 23/ 2) Piecewise[{{1., -1. < w <  0.},
{0.3333333333333333 (3. - 1. w), 0. <= w < 1.}},  0.] && (β - x -19/2)
Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y <  1.}}, 0.] >
(β - y - 21/2) Piecewise[{{0.3333333333333333 (2. - 1. w), -1. < w < 1.}}, 0.]


the set of inequalities that I'm actually trying to solve is

(β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y < 1.}}, 0.] >
(β - z - 23/ 2) Piecewise[{{1., -1. < w <  0.},
{0.3333333333333333 (3. - 1. w), 0. <= w < 1.}},  0.] && (β - x -19/2)
Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y <  1.}}, 0.] >
(β - z - 21/2) Piecewise[{{0.3333333333333333 (2. - 1. w), -1. < w < 1.}}, 0.]

• Please, clearly indicate the changes in your question when editing it as civilized people use to do. Your set of inequalities was substantially changed. – user64494 Apr 28 '18 at 13:04
• I'm sorry, there was a mistake in the first block of code that I edited yesterday, but I don't think it was substantial. In the first part of the addition I reported the (wrong) equation that works with reduce. In the second part of the addition I reported the (right) equation that does not work with reduce. I thought explaining in the changes in title of the edit was enough. My bad. – Rby Apr 28 '18 at 13:14
• Sorry, you still didn't indicate the changes in your edit. – user64494 Apr 28 '18 at 13:17
• Hi @Rby, I could use some feedback. Have you had a look at my answer? Did I understand your question or not? Does it work as expected? Cheers! – AccidentalFourierTransform Apr 30 '18 at 21:12
• @Rby Well, to have a symbolic expression instead of a numerical one, you could replace NIntegrate for Integrate. But the computation takes so long that I aborted it. You can try for yourself if you want to, but I really think numerical is the way to go here -- the symbolic result won't be illuminating. – AccidentalFourierTransform May 1 '18 at 15:23

Unless I misunderstood the question, this should do: Let

ineq = ((β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y < 1.}}, 0.] >
(β - z - 23/2) Piecewise[{{1., -1. < w <  0.},
{0.3333333333333333 (3. - 1. w), 0. <= w < 1.}},  0.] && (β - x -19/2)
Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y <  1.}}, 0.] >
(β - z - 21/2) Piecewise[{{0.3333333333333333 (2. - 1. w), -1. < w < 1.}}, 0.]);


With this, the probability is, by definition,

NIntegrate[Boole[ineq && 1 > z >= w > -1 && 1 > x >= y > -1], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, {w, -1, 1}, {β, 8.5, 11.5}, Method -> "MonteCarlo"]/48


which evaluates to $0.22$, in less than 0.01 seconds on my laptop.

If you repeat the calculation a thousand times (to have better statistics), then you get a probability of $0.226(2)$, in around a minute on my laptop.

--

In a comment to user64494's post, OP says that they want the probability as a function of x, y, z, w. This is trivial to implement: we just drop the integration over such parameters:

probability[x_?NumericQ, y_?NumericQ, z_?NumericQ, w_?NumericQ] :=
NIntegrate[
Boole[(((β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y < 1.}}, 0.] > (β - z - 23/2) Piecewise[{{1., -1. < w < 0.}, {0.3333333333333333 (3. - 1. w), 0. <= w < 1.}}, 0.] && (β - x - 19/2) Piecewise[{{0.3333333333333333 (1. - 1. y), -1. < y < 1.}}, 0.] > (β - z - 21/2) Piecewise[{{0.3333333333333333 (2. - 1. w), -1. < w < 1.}}, 0.])) && 1 > z >= w > -1 && 1 > x >= y > -1]
, {β, 8.5, 11.5}]/48


so that, for example, probability[.3, .3, .3, .3] evaluates to $0.06$.

We cannot plot a function of four variables; but if we integrate out x, y, z, and plot the probability as a function of w, we get

Table[
NIntegrate[
Boole[ineq && 1 > z >= w > -1 && 1 > x >= y > -1]
, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, {β, 8.5, 11.5}, Method -> "MonteCarlo"]/48
, {w, -1, 1, 1/20}] // Quiet // ListPlot


Needless to say, the area under the curve is $0.22$, as before.

• Can't reproduce it: I obtain NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 0. and 0. for the integral and error estimates. – user64494 Apr 29 '18 at 17:11
• @user64494 I was missing a parenthesis, thank you for the heads-up! It should work now; please let me know if it doesn't. Cheers! – AccidentalFourierTransform Apr 29 '18 at 17:29
• However, this does not answer the question: OP writes "But the probability should be a function of the parameters x,y,z,w not a number." in a comment to my answer. – user64494 Apr 29 '18 at 19:35
• @user64494 Hmm I didn't see that comment before -- the OP is slightly confusing anyway. Thank you for letting me know; I'll fix the answer in a minute. – AccidentalFourierTransform Apr 29 '18 at 19:37
• The probability is a function of all four parameters simultaneously. See the addition in my answer. – user64494 Apr 29 '18 at 20:17

If I correctly understand the question, the probability likely equals 1. Rationalizing the coefficients, we produce a table by

kk = Table[FindInstance[(\[Beta] - x -
19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}},
0] > (\[Beta] - y - 23/2)
Piecewise[{{1, -1 < w < 0}, {1/3 (3 - w), 0 <= w < 1}},
0] && (\[Beta] - x -
19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}},
0] > (\[Beta] - y -
21/2) Piecewise[{{1/3 (2 - w), -1 < w < 1}}, 0] &&
1 > z >= w > -1 && 1 > x >= y > -1 /. \[Beta] -> a, {x, y, z, w},
Reals], {a, 17/2, 23/2, 1/10}]


{{{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> 0, y -> 0, z -> 0, w -> 0}}, {{x -> -(1/2), y -> -(1/2), z -> -(1/4), w -> -(1/4)}}}

Length[kk]


31

Putting the step equal to $\frac 1 {100}$, we obtain the length equals $301$.

Probability[(\[Beta] - x -
19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}},
0] > (\[Beta] - y -
23/2) Piecewise[{{1, -1 < w < 0}, {1/3 (3 - w), 0 <= w < 1}},
0] && (\[Beta] - x - 19/2) Piecewise[{{1/3 (1 - y), -1 < y < 1}},
0] > (\[Beta] - y - 21/2) Piecewise[{{1/3 (2 - w), -1 < w < 1}},
0] && 1 > z >= w > -1 &&  1 > x >= y > -1,
\[Beta] \[Distributed]UniformDistribution[{17/2, 23/2}]]


\ [Piecewise] 1 (-1-1&&y>-1&&z<1&&x-y>=0&&w-x>=0&&w-z<=0)||(-1-1&&z<1&&x-y>=0&&w-x==0&&w-z<=0)||(w==0&&-1-1&&0<=z<1&&x-y>=0)||(w==0&&x==0&&-1-1&&y>-1&&z<1&&x-y>=0&&w-z<=0&&w+2 x<=0)||(0=0&&w-x==0&&w-z<=0&&3/(-3+w-x)-(2 w)/(-3+w-x)-x/(-3+w-x)-y<0)||(0=0&&w-x>=0&&w-z<=0&&3/(-3+w-x)-(2 w)/(-3+w-x)-x/(-3+w-x)-y<0&&w+2 x>0)||(00&&w-z<=0&&3/(-3+w-x)-(2 w)/(-3+w-x)-x/(-3+w-x)-y==0&&w+2 x>0) (-3+2 w+x-3 y+w y-x y)/(3 (-1+w-y)) (-1-1&&w-x<0&&z<1&&w-z<=0&&w+2 x<=0&&x-y>=0)||(-10&&x-y>=0&&3/(-3+w-x)-(2 w)/(-3+w-x)-x/(-3+w-x)-y<0&&x<1)||(0=0&&3/(-3+w-x)-(2 w)/(-3+w-x)-x/(-3+w-x)-y<0&&x<1) (3-x+3 y+x y)/(3 (1+y)) w==0&&0=0&&0<=z<1 (w-x-w y+x y)/(3 (-1+w-y)) 00&&w-x>0&&y>-1&&3/(-3+w-x)-(2 w)/(-3+w-x)-x/(-3+w-x)-y>0&&w-z<=0&&z<1 0 True

Notice. It should be noticed that my answer was done to the question under consideration before its edit where the variable z was included in the first three inequalities. That was not indicated in the question as civilized people use to do.

• But the probability should be a function of the parameters x,y,z,w not a number. Moreover if the parameter space is just 2: x and y it takes just a few seconds to get the answers – Rby Apr 27 '18 at 21:19
• @Rby: Sorry, but your question is unclearly formulated. Please, elaborate it. Replacing FindInstance by Reduce in the code from my answer, one obtains a table of systems of inequalities. Hope, that would be useful. – user64494 Apr 28 '18 at 3:45
• Using Reduce instead of FindInstance gives me a set of inequalities however that takes 275 seconds on my computer. If I need to add more constraints I fear it would take even more. I thought Mathematica was better at handling symbols rather than numbers, I don't get why with 2 parameters it takes only a couple of seconds and with 4 the time explodes. With Matlab I had 18 parameters and the calulation took something around 30 seconds. This problem does not look that hard to me – Rby Apr 28 '18 at 9:09
• @Rby When you did it in Matlab with 18 parameters, was that calculation symbolic or numerical? – MarcoB Apr 30 '18 at 2:12
• @Rby It is a lot better, but symbolic calculations are orders of magnitude more complex than the corresponding numerical ones, so they are significantly limited, not by the software you are using, but by the nature of the problem. To give an example, solution of a quadratic equation is straightforward symbolically, that of a cubic is much more complex, and higher powers quickly get you to where symbolic solutions simply do not exist. Yet, obtaining a numerical solution of such equations is very fast, almost independently of the symbolic complexity. – MarcoB Apr 30 '18 at 13:24