MatLab - Numerical Integration given Conditions - Translate from Mathematica

Question

This doesn't want to work for me in Mathematica - 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.

Would someone be able to help me translate this to MatLab? I am new to that program, and it would save me time if someone with experience could help me get started.

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}]

For those not familiar to Mathematica, Boole simply gives 1 if the condition is met, and 0 otherwise. As such, I wish to integrate over all of parameterspace to find the volume of the subspace for which the conditions are met.

Answer 1

Here's another attempt, this time using Monte Carlo integration. Values within the integration limits are randomly generated, the function is evaluated at those points, and the integral is estimated as the fraction of "hits" multiplied with the total volume of the integration area.

I've implemented it such that random numbers are generated in chunks of 10 million, and the estimate is iteratively updated for each chunk. The function runs forever, and should be stopped once the estimate appears to be precise enough.

After a few minutes of running, it appears that the value of the integral is close to 0.001775.

function LauBo2

myO = @(e1, e2) e1 ./ (e1 + e2);
myU = @(e1, e2) (1 - e2) ./ (2 - e1 - 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);

    function I = integrand(p, e1, e2, e3)
        I = double(e3 < myV(p, e1, e2)) .* (myV(p, e1, e2) < 1 - e3) .* ...
            (1 - e3 < myW(p, e1, e2)) .* (myO(e1, e2) < e3) .* ...
            (e3 < myU(e1, e2)) .* (e1 < p) .* ( p < 1 - e1) .* (e1 < e2);
    end

n = 10000000;
S = 0;
N = 0;
while true
    p = rand(n, 1);         % p from 0 1
    e1 = 0.5 * rand(n, 1);  % e1 from 0 to 0.5
    e2 = 0.5 * rand(n, 1);  % e2 from 0 to 0.5
    e3 = 0.5 * rand(n, 1);  % e3 from 0 to 0.5

    S = S + mean(integrand(p, e1, e2, e3));
    N = N + 1;

    I = S / N  * 1 * 0.5 * 0.5 * 0.5;
    fprintf('%.20f\n', I)
end

end