1
$\begingroup$

Let's say I have two functions $f(x,y,z), g(x,y,z)$. I would like to know the probability that $f(x,y,z) < g(x,y,z)$ over some defined, numerical ranges of $x, y, z$.

Any thoughts on how to do it besides nested for-loops?

P.S. I have more than 3 variables, but for simplicity let's assume I only have $x,y,z$.

Edit:

dL = 1; rL = 1; c = 0; 
f[dH_, rH_, α_, β_, γ_, 
   b_] := -(1/
     2) (-1 + β) Sqrt[-(((-1 + α) (dH (-1 + β) - \
β (1 + rH γ)^2))/(b^2 dH (-1 + β)))] + 
   1/2 β Sqrt[(α (dH - 
       dH β + β (1 + rH γ)^2))/(
    b^2 dH^2 β)];
g[dH_, rH_, α_, β_, γ_, b_] := 1/(
  2 (b + b (-1 + dH) β));

reg1 = ImplicitRegion[
   0 <= α <= 1 && 0 <= β <= 1 && 0 <= b <= 1 && 
    0 <= γ <= 1 && 1 <= dH <= 5 && 1 <= rH <= 5, {dH, 
    rH, α, β, γ, b}];

reg2 = ImplicitRegion[
   f[dH, rH, α, β, γ, b] > 
     g[dH, rH, α, β, γ, b] && 0 <= α <= 1 &&
     0 <= β <= 1 && 0 <= b <= 1 && 0 <= γ <= 1 && 
    1 <= dH <= 5 && 1 <= rH <= 5, {dH, 
    rH, α, β, γ, b}];

prob = RegionMeasure[reg2]/RegionMeasure[reg1] // N (* Takes forever to calculate RegionMeasure[reg2] *)

NIntegrate[1, {dH, rH, α, β, γ, b} ∈ 
  reg1] (* 16 *)

NIntegrate[1, {dH, rH, α, β, γ, b} ∈ 
  reg2, Method -> "AdaptiveQuasiMonteCarlo"] (* returns the function call. *)

$\endgroup$
4
  • $\begingroup$ You can either use something like Integrate[Boole[inequalities], {x, ...}, ...] or reg = ImplicitRegion[inequalities, {x, y, z}]; Integrate[1, {x, y, z} \[Element] reg]. An example of your inequalities would help. $\endgroup$
    – Carl Woll
    Commented Oct 5, 2019 at 21:49
  • $\begingroup$ There are "MonteCarlo" and "AdaptiveMonteCarlo" integration strategies for NIntegrate. reference.wolfram.com/language/tutorial/… The Probability is equal to the volume $\Omega_{f<g}$ divided by the total volume $\Omega$, that is $p=\Omega_{f<g}/\Omega$, where $\Omega$ can be expressed as a multi- (3 in your example) dimensional integral. $\endgroup$
    – yarchik
    Commented Oct 5, 2019 at 22:00
  • $\begingroup$ @JohnSmith The best integration scheme depends mostly on the exact dimension and the complexity of your functions. If either the dimension is >=8 or there is no big chance that the integrals can be done analytically (via cylindrical algebraic decomposition), AdaptiveMonteCarlo and AdaptiveQuasiMonteCarlo are good NIntegrate Methods to try (see CarlWolls comment above). $\endgroup$ Commented Oct 6, 2019 at 11:36
  • $\begingroup$ Hi guys, thanks so much for the suggestions! I've tried NIntegrate with AdaptiveQuasiMonteCarlo, AdaptiveMonteCarlo, MonteCarlo but with no luck. Mathematica just returns the function call. I'll update the question with specific inequalities. $\endgroup$
    – John Smith
    Commented Oct 6, 2019 at 23:21

2 Answers 2

2
$\begingroup$

Update: With your function f and g we van calculate probability of f>g in the domain as:

      dL = 1; rL = 1; c = 0;
f[dH_, rH_, α_, β_, γ_, 
   b_] := -(1/
       2) (-1 + β) Sqrt[-(((-1 + α) (dH (-1 + β) - \
β (1 + rH γ)^2))/(b^2 dH (-1 + β)))] + 
   1/2 β Sqrt[(α (dH - 
          dH β + β (1 + 
              rH γ)^2))/(b^2 dH^2 β)];
        g[dH_, rH_, α_, β_, γ_, b_] :=  1/(2 (b + b (-1 + dH) β));

n = 100000;
{dH, rH} = RandomReal[{1, 5}, {2, n}];
{α, β, γ, b} = RandomReal[{0, 1}, {4, n}];
pts = Transpose[{dH, rH, α, β, γ, b}];
Total@UnitStep[f @@@ pts - g @@@ pts]/n // N

0.85568

Original Answer: How about this?

f[x_, y_, z_] := x + y + z
g[x_, y_, z_] := 2 x + 2 y + 2 z

n = 1000000;

{aX, bX} = {-5, 5};
{aY, bY} = {-3, 6};
{aZ, bZ} = {2, 7};
rangeX = RandomReal[{aX, bX}, n];
rangeY = RandomReal[{aY, bY}, n];
rangeZ = RandomReal[{aZ, bZ}, n];

pts = Transpose[{rangeX, rangeY, rangeZ}];

Total@UnitStep[g @@@ pts - f @@@ pts]/n // N

0.920215

$\endgroup$
1
  • $\begingroup$ Thanks! This will work no matter how complex the inequalities are. $\endgroup$
    – John Smith
    Commented Oct 7, 2019 at 0:10
3
$\begingroup$

Using the functions provided by OkkesDulgerci

Clear["Global`*"];

f[x_, y_, z_] := x + y + z;
g[x_, y_, z_] := 2 x + 2 y + 2 z;

{aX, bX} = {-5, 5};
{aY, bY} = {-3, 6};
{aZ, bZ} = {2, 7};

reg1 = ImplicitRegion[
   aX <= x <= bX && aY <= y <= bY && aZ <= z <= bZ, {x, y, z}];

reg2 = ImplicitRegion[
   f[x, y, z] < g[x, y, z] &&
    aX <= x <= bX && aY <= y <= bY && 
    aZ <= z <= bZ, {x, y, z}];

prob = Volume[reg2]/Volume[reg1] // N

(* 0.92037 *)

or

prob = RegionMeasure[reg2]/RegionMeasure[reg1] // N

(* 0.92037 *)

The region is

RegionPlot3D[f[x, y, z] < g[x, y, z],
 {x, aX, bX}, {y, aY, bY}, {z, aZ, bZ},
 BoxRatios -> {bX - aX, bY - aY, bZ - aZ},
 PlotStyle -> Opacity[0.75]]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.