# How to represent in a 3D graph the set of solutions to this system of inequalities?

I am studying this piecewise function:

Piecewise[{{f1[x, h], x < a}, {f2[x, a, h], a <= x < b},
{f3[x, a, b, h], b <= x < c}, {f4[x, a, b, c, h],  x >= c}}]


where

f1[x_,h_] := pdf1[x]/pdf[x] (h + (1 - h) 2 cdf1[x])/(h + (1 - h) 2 cdf[x])
f2[x_,a_,h_] := pdf1[x]/pdf[x] (h + (1 - h) 2 (1 - cdf1[x] + cdf1[a]))/
(h + (1 - h) 2 (1 - cdf[x] + cdf[a]))
f3[x_,a_,b_,h_] := pdf1[x]/pdf[x] (h + (1 - h) 2 (cdf1[x] - cdf1[b] + cdf1[a]))/
(h + (1 - h) 2 (cdf[x] - cdf[b] + cdf[a]))
f4[x_,a_,b_,c_] := pdf1[x]/pdf[x] (h + (1 - h) 2 (1 - cdf1[x] + cdf1[c] - cdf1[b] +
cdf1[a]))/(h + (1 - h) 2 (1 - cdf[x] + cdf[c] - cdf[b] + cdf[a]))


where, in turn

pdf1[x_] := PDF[NormalDistribution[1, 1], x];
pdf[x_] := PDF[NormalDistribution[0, 1], x];
cdf1[x_] := CDF[NormalDistribution[1, 1], x];
cdf[x_] := CDF[NormalDistribution[0, 1], x];


I want the first and the third piece of the piecewise function to be smaller than a certain value $$l$$, say $$l=0.48$$; and the second and fourth piece of the function to be greater than the same $$l$$. The following code represents the intervals of $$a$$, $$b$$ and $$c$$ such that the four inequalities are satisfied:

{amin = NArgMin[{f2[a], f2[a] > l}, a], amax = NArgMax[{f1[a], f1[a] < l}, a],
a = RandomReal[{amin, amax}]}
{bmin = a, bmax = NArgMax[{f3[b], f3[b] < l}, b], b = RandomReal[{bmin, bmax}]}
{cmin = b, cmax = NArgMax[{f3[c], f3[c] < l}, c], c = RandomReal[{cmin, cmax}]}


Note however that, for example, you have to give a random value to $$a$$ since the lower bound of the interval of the $$b$$'s depends on $$a$$.

I would like to plot in a 3D Graph all the possible combinations $$(a,b,c)$$ such that the system of inequalities is satisfied. Moreover, I would like to implement the Manipulate command so that I can study how the set of all possible combinations $$(a,b,c)$$ changes as $$h$$ varies between 0 and 1. If a 3D graph was not feasible, actually any graphical representation of the set of solutions of $$(a,b,c)$$ would be fine anyway. Thank you if you can help.

• – kglr Oct 4 '18 at 23:11

The requested plot can be obtained as follows. First, compute the limits on a, using the expressions in the question.

{amin = NArgMin[{f2[a], f2[a] > l}, a], amax = NArgMax[{f1[a], f1[a] < l}, a]}
(* {-0.233969, 0.544866} *)


Next, compute the corresponding limits in b.

btab = Interpolation@Table[{a, NArgMax[{f3[b], f3[b] < l}, b]}, {a, amin, amax, .02}];


and finally in c.

ctab = Interpolation[Flatten[Table[{{a, b}, Quiet@NArgMax[{f3[c], f3[c] < l, b < c},
c]}, {a, amin, amax, .02}, {b, a, btab[a], .02}], 1], InterpolationOrder -> 1];


Since a < b < c, the plot is given by

Plot3D[{b, ctab[a, b]}, {a, amin, amax}, {b, a, btab[a]}, AxesLabel -> {a, b, c},
LabelStyle -> {Black, Bold, Medium}, BoxRatios -> Automatic, ImageSize -> Large,
Mesh -> None, PlotStyle -> Opacity[.5]] where the allowed values of c lie between the upper and lower curves (both of which have the same projection on the {a, b} plane). The limits on a, b, and c were obtained here.

As noted by kgir in a comment above, RegionPlot3D also can be used once the logical structure and interpolation functions have been determined. RegionPlot3D is much slower here than Plot3D but produces a figure that some may prefer. The choice is a matter of taste.

Quiet@RegionPlot3D[amin < a < amax && a < b < btab[a] && b < c < ctab[a, b],
{a, amin, amax}, {b, amin, btab[amin]}, {c, amin, ctab[amin, btab[amin]]},
AxesLabel -> {a, b, c}, LabelStyle -> {Black, Bold, Medium},
BoxRatios -> Automatic, ImageSize -> Large, PlotPoints -> 500, Mesh -> None] As noted by the OP in a comment below, the computation of ctab given here often fails with the error

Interpolation::femimq The element mesh has insufficient quality of 0..

as h is increased. Basically, the triangular grid of Interpolation becomes tangled near the edge at which the two curves in the first plot come together. A simply way of eliminating the error is to replace 0.02 by 0.06 in the computation of ctab. However, doing so decreases the accuracy of the plot significantly for 1/3 < h < 1. A much more accurate but slower computation involves replacing the three lines of code immediately preceding the first plot by

btab = Interpolation@Table[{a, NArgMax[{f3[b], f3[b] < l}, b]}, {a, amin, amax},
(amax - amin)/20}];
fc[a0_, b0_] := Quiet@NArgMax[{f3[c], f3[c] < l, b < c} /. {a -> a0, b -> b0}, c]
Plot3D[{b, fc[a, b]}, {a, amin, amax}, {b, a, btab[a]}, AxesLabel -> {a, b, c},
LabelStyle -> {Black, Bold, Medium}, BoxRatios -> Automatic, ImageSize -> Large,
Mesh -> None, PlotStyle -> Opacity[.5], PlotPoints -> 10, MaxRecursion -> 0]


The resulting plot for h = 1/2 is • Thank you very much! I have a couple of questions: which value of $h$ did you use to get that plot? The limits of $a$, $b$ and $c$ depends on the $h$ you use. Is there a way to implement something like a Manipulate command so that I can study how the graph changes as $h$ varies? – Api Oct 5 '18 at 22:40
• I used h = 0, but other values of h undoubtedly would work. Manipulate would be impractical in this case, because the computation requires several minutes, especially for the second plot. Please note that I added b < c to NArgMax[{f3[c], f3[c] < l, b < c}, c] from the other question, which I also updated, to treat unusual cases where f3 is not continuous. Incidentally, I soon will replace the two figures with ones of higher quality. – bbgodfrey Oct 5 '18 at 22:54
• Ok, thank you for your clarification! I will just produce several graphs for different values of h and then try to come up with some nice way of presenting them. I will check for the higher quality figures then. Thank you for the time you are spending helping me. – Api Oct 6 '18 at 9:43
• @Api The problem is with the ctab Interpolation grid near the edge where the two surfaces meet. A quick fix is to replace 0.02 by 0.06 in the definitions of btab and ctab, but the plot less accurate for larger h`.. I would like to find a better solution, however. – bbgodfrey Oct 20 '18 at 2:10
• @Api I have provided the code, but it is slow. – bbgodfrey Dec 9 '18 at 22:03