1
$\begingroup$

This doesn't want to work for me - I don't know why not. It usually gives me 0, but sometimes also a different number. I know it should be a small, non-zero number.

myO[e1_, e2_] := (e1)/((e1) + (e2));
myU[e1_, e2_] := (1 - e2)/((1 - e1) + (1 - e2));
myV[p_, e1_, e2_] := p (e1)/(p (e1) + (1 - p) (1 - e2));
myW[p_, e1_, e2_] := p (1 - e1)/(p (1 - e1) + (1 - p) e2);

NIntegrate[ Boole[0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1 && 
                  e3 < myV[p, e1, e2] < 1 - e3 < myW[p, e1, e2] && 
                  myO[e1, e2] < e3 < myU[e1, e2]], 
                 {p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}]
$\endgroup$
4
  • $\begingroup$ Do you mean that sometimes the same command gives different outputs with everything else the same ? $\endgroup$ Dec 6 '13 at 15:13
  • $\begingroup$ Yes, it does. But more importantly, it doesn't give the correct value. It mostly gives zero, which cannot be true. $\endgroup$
    – LBogaardt
    Dec 6 '13 at 15:30
  • 1
    $\begingroup$ In version 9 I get the same result every time. If I change the method to MonteCarlo, I get a slightly different number every time, as expected. Raising MaxPoints makes it converge to the same result the default method gives: 0.00177529. In version 8 I get the same. What version are you using? $\endgroup$
    – Szabolcs
    Dec 7 '13 at 1:06
  • $\begingroup$ I am using Version 9. $\endgroup$
    – LBogaardt
    Dec 9 '13 at 15:59
1
$\begingroup$

You can use Reduce or CylindricalDecomposition to break down the inequalities into forms that can be used as limits of integration. (Note: One might try BooleanMinimize instead of distr below to distribute And over Or, but it changes the order of the inequalities. Consequently, the intervals of integration end up out of order and NIntegrate does not work. One could sort them, but I opted for writing distr, which distributes from the bottom up.)

 (* distribute And over Or *)
distr[e_And] := Distribute[distr /@ e, Or];
distr[e_Or] := distr /@ e;
distr[e_] := e;
 (* convert ineqs to limits of integration *)
sub = {(Less | LessEqual)[a_, b_, c_] :> {b, a, c}, 
   HoldPattern[Inequality[a_, _, b_, _, c_]] :> {b, a, c}};

simp = Reduce[
    0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1 && 
      e3 < myV[p, e1, e2] < 1 - e3 < myW[p, e1, e2] && 
      myO[e1, e2] < e3 < myU[e1, e2],
    {e1, e2, e3, p}];
NIntegrate[1, ##] & @@@ (distr[simp] /. sub /. {And | Or -> List}) // Total

(* 0.00177529 *)
$\endgroup$
1
  • $\begingroup$ Thanks! That's great. $\endgroup$
    – LBogaardt
    Dec 9 '13 at 16:01
1
$\begingroup$

One way to proceed is to try out portions and see where the problem lies. For instance:

NIntegrate[Boole[0 < e1 < e2 < 1/2 && 0 < e3 < 1/2 && e1 < p < 1 - e1], 
     {p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}]
0.0416667

is fine. But the inequalities

NIntegrate[Boole[e3 < myV[p, e1, e2] < 1 - e3 < myW[p, e1, e2]], 
     {p, 0, 1}, {e1, 0, 1/2}, {e2, 0, 1/2}, {e3, 0, 1/2}]
0.0239936 - 3.06317*10^-29 I

are giving (small) complex values. Most likely, this is because there are points where your function is trying to divide by zero. One culprit is the function myW, which, at 0,0,0 gives

myW[0, 0, 0]
Power::infy: "Infinite expression 1/0 encountered."
$\endgroup$

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .