# Sampling an inequality over multi-dimensional space

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. *)


• 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. Oct 5, 2019 at 21:49
• 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. Oct 5, 2019 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). Oct 6, 2019 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. Oct 6, 2019 at 23:21

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

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

• Thanks! This will work no matter how complex the inequalities are. Oct 7, 2019 at 0:10

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