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

Answer 2

I believe that this should be what you are looking for. It's not very elegant, Matlab isn't really made for this. I broke down your conditions (Matlab doesn't support x < y < z), eliminated those that are redundant with the integral limits, and transformed others into changed integral limits. The nested integral is implemented by a series of functions that evaluate one integral and serve as integrand for the next. The for loops are necessary because integral expects the given function to be vectorized.

I can't tell you whether something sensible comes out of this, because it takes ages to compute and I didn't have the patience to wait till it is finished.

function I = LauBo

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 = first(p, e1, e2, e3)
        % integrand
        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));
    end

    function I = second(e1, e2, e3)
        % integrate over  p   from e1 to (1 - e1)
        I = nan(size(e1));
        for i = 1 : numel(e1)
            I(i) = integral(@(p) first(p, e1(i), e2, e3), e1(i), 1 - e1(i));
        end
    end

    function I = third(e2, e3)
        % integrate over  e1  from 0 to e2
        I = nan(size(e2));
        for i = 1 : numel(e2)
            I(i) = integral(@(e1) second(e1, e2(i), e3), 0, e2(i));
        end
    end

    function I = fourth(e3)
        % integrate over  e2  from 0 to 0.5
        I = nan(size(e3));
        for i = 1 : numel(e3)
            I(i) = integral(@(e2) third(e2, e3(i)), 0, 0.5);
        end
    end

% integrate over  e3  from 0 to 0.5
I = integral(@fourth, 0, 0.5);

end