MatLab - 给定条件下的数值积分 - 翻译自 Mathematica

Question

这不想在 Mathematica 中为我工作 - 我不知道为什么不。它通常会给我

0

，但有时也会给我不同的数字。我知道它应该是一个很小的非零数字。

有人可以帮我将其翻译成MatLab吗？我是该程序的新手，如果有经验的人可以帮助我开始，这会节省我的时间。

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

对于那些不熟悉 Mathematica 的人来说，如果满足条件，

Boole

只是给出

1

，否则给出

0

。因此，我希望对所有参数空间进行积分，以找到满足条件的子空间的体积。

Answer 1

这是另一种尝试，这次使用蒙特卡罗积分。积分限制内的值是随机生成的，在这些点处评估函数，积分估计为“命中”的分数乘以积分区域的总体积。

我已经实现了它，以 1000 万个块为单位生成随机数，并且针对每个块迭代更新估计值。该函数会永远运行，一旦估计看起来足够精确就应该停止。

运行几分钟后，积分值似乎接近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

我相信这应该就是您正在寻找的。它不是很优雅，Matlab 并不是专门为此而设计的。我分解了你的条件（Matlab 不支持 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

期望给定的函数被向量化。

我无法告诉你这是否会产生一些明智的结果，因为计算需要很长时间，而且我没有耐心等到它完成。

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