1

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. *)
OkkesDulgerci
  • 10,716
  • 1
  • 19
  • 38
John Smith
  • 173
  • 5
  • 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. – Carl Woll Oct 05 '19 at 21:49
  • There are "MonteCarlo" and "AdaptiveMonteCarlo" integration strategies for NIntegrate. https://reference.wolfram.com/language/tutorial/NIntegrateIntegrationStrategies.html 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. – yarchik Oct 05 '19 at 22:00
  • @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). – Thies Heidecke Oct 06 '19 at 11:36
  • 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. – John Smith Oct 06 '19 at 23:21

2 Answers2

3

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

Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
2

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

OkkesDulgerci
  • 10,716
  • 1
  • 19
  • 38