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