# Numerical differentiation of summation of implicit regions

I'm very new to mathematica and this forum, so I apologize for the poor formatting. I am trying to find a probability density function, given a continuous and piecewise differentiable cumulative distribution function (so taking the derivative). The cumulative distribution function is defined by the 5 regions, which I define to be "TestRegioni" for i ranging from 1 to 5. Note that s1,s2,s3,h1,h2,h3 are all constants, and t is a time parameter with respect to x and y.

TestRegion1[t_]:= ImplicitRegion[(1-(x*s1))/(h1)-2*x*(s1)/(h1)<y&&y<(1-(x*s1))/(h1)&&  y>(1-(2+s2)*x)/(h2)&&h1/(x*(s1*x+h1*y))<t &&  0<x&&x<1,{{x, -2, 2}, {y, -4, 4}}]

TestCDF1[t_]:=NIntegrate[1, {x,y}\[Element]TestRegion1[t]];

TestRegion2[t_]:= ImplicitRegion[(1-(x*s1))/(h1)-2*x*(s1)/(h1)<y&&y<(1-(x*(2+s2)))/(h2)&&  y>(1-(4*s1+s3)*x)/(h3)&&h2/(x*((2+s2)*x+h2*y))<t &&  0<x&&x<1,{{x, -2, 2}, {y, -4, 4}}]

TestCDF2[t_]:=NIntegrate[1, {x,y}\[Element]TestRegion2[t]]

TestRegion3[t_]:= ImplicitRegion[ (1-(x*(4*s1+3*s3)))/(h3)-2*x*(s1)/(h1)<y&&y<(1-(x*(4*s1+s3)))/(h3)&&y>(1-(s1)*x)/(h1) - 2*x*(s1)/(h1)&&h3/(x*((4*s1+s3)*x+h3*y))<t &&  0<x&&x<1,{{x, -2, 2}, {y, -4, 4}}]

TestCDF3[t_]:=NIntegrate[1, {x,y}\[Element]TestRegion3[t]]

TestRegion4[t_]:= ImplicitRegion[ (1-(x*(2 + 3*s2)))/(h2)<y&&y< (1-(x*(4*s1+3*s3)))/h3&&
y>(1-(s1)*x)/(h1)  - 2*x*(s1)/(h1)&&h3/(x*((4*s1+3*s3)*x+h3*y))<t &&  0<x&&x<1,{{x, -2, 2}, {y, -4, 4}}]

TestCDF4[t_]:=NIntegrate[1, {x,y}\[Element]TestRegion4[t]]

TestRegion5[t_]:= ImplicitRegion[ (1-(x*s1))/(h1)-2*x*(s1)/(h1)<y&&y< (1-(x*(2+3*s2)))/h2&&h2/(x*((2+3*s2)*x+h2*y))<t &&  0<x&&x<1,{{x, -2, 2}, {y, -4, 4}}]

TestCDF5[t_]:=NIntegrate[1, {x,y}\[Element]TestRegion5[t]]

TotalCDF[t_]:=TestCDF1[t]+TestCDF2[t]+TestCDF3[t]+TestCDF4[t]+TestCDF5[t];


At this point I am getting no errors. Everytime I evaluate these functions at some value t, I get a numerical value, so my functions make sense. I run into an error when I then define:

TestPdf[s_] := ND[TotalCDF[t],t, s];


This is where I am trying to take the numerical derivative of the cumulative distribution function. The error I get is:

"DiscretizeRegion: DiscretizeRegion was unable to discretize the region ImplicitRegion[<<2>>]."

I'm not sure why I am getting this. Each of the functions evaluate correctly. Maybe it has to do with the fact that I am working with 3 variables (x,y,t)? Maybe it has to do with syntax? My ultimate goal is to Plot the probability density function, but I can't get that far because of this error. Also, I am on version 13.3.1 for Linux x86.

• I'm missing the values of s1,s2,s3,h1,h2,h3 Commented Dec 20, 2023 at 13:59
• The values are s1 = Cos[Pi/7], s2 = Cos[2*Pi/7], s3 = Cos[3*Pi/7], h1 = Sin[Pi/7], h2 = Sin[2*Pi/7], h3 = Sin[3*Pi/7] Commented Dec 21, 2023 at 3:50

Use ?NumericQ and Module to evaluate pdf as follows

TestPdf[s_?NumericQ] :=
Module[{s1 = 1, s2 = 2, s3 = 3, h1 = 1, h2 = 2, h3 = 3},
TestRegion1[t_?NumericQ] :=
ImplicitRegion[(1 - (x*s1))/(h1) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*s1))/(h1) && y > (1 - (2 + s2)*x)/(h2) &&
h1/(x*(s1*x + h1*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF1[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion1[t]];

TestRegion2[t_?NumericQ] :=
ImplicitRegion[(1 - (x*s1))/(h1) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*(2 + s2)))/(h2) && y > (1 - (4*s1 + s3)*x)/(h3) &&
h2/(x*((2 + s2)*x + h2*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF2[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion2[t]];

TestRegion3[t_?NumericQ] :=
ImplicitRegion[(1 - (x*(4*s1 + 3*s3)))/(h3) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*(4*s1 + s3)))/(h3) &&
y > (1 - (s1)*x)/(h1) - 2*x*(s1)/(h1) &&
h3/(x*((4*s1 + s3)*x + h3*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF3[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion3[t]];

TestRegion4[t_?NumericQ] :=
ImplicitRegion[(1 - (x*(2 + 3*s2)))/(h2) < y &&
y < (1 - (x*(4*s1 + 3*s3)))/h3 &&
y > (1 - (s1)*x)/(h1) - 2*x*(s1)/(h1) &&
h3/(x*((4*s1 + 3*s3)*x + h3*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF4[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion4[t]];

TestRegion5[t_?NumericQ] :=
ImplicitRegion[(1 - (x*s1))/(h1) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*(2 + 3*s2)))/h2 &&
h2/(x*((2 + 3*s2)*x + h2*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF5[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion5[t]];

TotalCDF[t_?NumericQ] :=
TestCDF1[t] + TestCDF2[t] + TestCDF3[t] + TestCDF4[t] +
TestCDF5[t]; {TotalCDF[s], Derivative[1][TotalCDF][s]}];

lst = TestPdf[#] & /@ Range[0, 2, .1]

T = Range[0, 2, .1];

ListLinePlot[{Transpose[{T, lst[[All, 1]]}],
Transpose[{T, lst[[All, 2]]}]},
PlotLegends -> {"TotalCDF", "TestPdf"}, PlotRange -> All,
Frame -> True, FrameLabel -> {"t", ""}]


Update 1. In a case of parameters like s1 = Cos[Pi/7], s2 = Cos[2*Pi/7], s3 = Cos[3*Pi/7], h1 = Sin[Pi/7], h2 = Sin[2*Pi/7], h3 = Sin[3*Pi/7] we use Rationalize[] and Interpolation[] for simplification as follows

TestCDF[s_?NumericQ] :=
Module[{s1 = Rationalize[Cos[Pi/7.], 10^-16],
s2 = Rationalize[Cos[2*Pi/7.], 10^-16],
s3 = Rationalize[Cos[3*Pi/7.], 10^-16],
h1 = Rationalize[Sin[Pi/7.], 10^-16],
h2 = Rationalize[Sin[2*Pi/7.], 10^-16],
h3 = Rationalize[Sin[3*Pi/7.], 10^-16]},
TestRegion1[t_?NumericQ] :=
ImplicitRegion[(1 - (x*s1))/(h1) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*s1))/(h1) && y > (1 - (2 + s2)*x)/(h2) &&
h1/(x*(s1*x + h1*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF1[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion1[t], AccuracyGoal -> 8,
PrecisionGoal -> 8];

TestRegion2[t_?NumericQ] :=
ImplicitRegion[(1 - (x*s1))/(h1) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*(2 + s2)))/(h2) && y > (1 - (4*s1 + s3)*x)/(h3) &&
h2/(x*((2 + s2)*x + h2*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF2[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion2[t], AccuracyGoal -> 8,
PrecisionGoal -> 8];

TestRegion3[t_?NumericQ] :=
ImplicitRegion[(1 - (x*(4*s1 + 3*s3)))/(h3) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*(4*s1 + s3)))/(h3) &&
y > (1 - (s1)*x)/(h1) - 2*x*(s1)/(h1) &&
h3/(x*((4*s1 + s3)*x + h3*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF3[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion3[t], AccuracyGoal -> 8,
PrecisionGoal -> 8];

TestRegion4[t_?NumericQ] :=
ImplicitRegion[(1 - (x*(2 + 3*s2)))/(h2) < y &&
y < (1 - (x*(4*s1 + 3*s3)))/h3 &&
y > (1 - (s1)*x)/(h1) - 2*x*(s1)/(h1) &&
h3/(x*((4*s1 + 3*s3)*x + h3*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF4[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion4[t], AccuracyGoal -> 8,
PrecisionGoal -> 8];

TestRegion5[t_?NumericQ] :=
ImplicitRegion[(1 - (x*s1))/(h1) - 2*x*(s1)/(h1) < y &&
y < (1 - (x*(2 + 3*s2)))/h2 &&
h2/(x*((2 + 3*s2)*x + h2*y)) < t && 0 < x &&
x < 1, {{x, -2, 2}, {y, -4, 4}}];

TestCDF5[t_?NumericQ] :=
NIntegrate[1, {x, y} \[Element] TestRegion5[t], AccuracyGoal -> 8,
PrecisionGoal -> 8];

TotalCDF[t_?NumericQ] :=
TestCDF1[t] + TestCDF2[t] + TestCDF3[t] + TestCDF4[t] +
TestCDF5[t]; TotalCDF[s]];


Example of usage

lst = Table[{t, TestCDF[t] // Chop}, {t, 0, 5, .025}]
cdf = Interpolation[lst, InterpolationOrder -> 4];

Plot[{cdf[t], cdf'[t]}, {t, 0, 5},
PlotLegends -> {"TotalCDF", "TestPdf"}, PlotRange -> All,
Frame -> True, FrameLabel -> {"t", ""}]


• It still was unable to discretize the region because, in my code, s1 = Cos[Pi/7], s2 = Cos[2*Pi/7], s3 = Cos[3*Pi/7], h1 = Sin[Pi/7], h2 = Sin[2*Pi/7], h3 = Sin[3*Pi/7], and these numbers are more complex than the constants you provided. Commented Dec 21, 2023 at 3:26
• Also, when I estimate these values to the nearest tenth, the integration fails due to the error "Numerical integration converging too slowly; suspect one of the following: singularity, the value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small." How do I take care of this? Especially if I want finer estimations. Commented Dec 21, 2023 at 4:39
• @CarsonNewman Please, see Update 1 to my post. Commented Dec 21, 2023 at 8:04
• Fantastic! And to make the orange line sharper I just have to increase the interpolation and Precision goals? Commented Dec 21, 2023 at 18:08
• @CarsonNewman It could be better to decrease step in lst = Table[{t, TestCDF[t] // Chop}, {t, 0, 5, .025}], the main error coming from here. Use for example, lst = Table[{t, TestCDF[t] // Chop}, {t, 0, 5, .005}]. Don't pay attention to PrecisionGoal up to, {t, 0, 5, .0000001} . Commented Dec 21, 2023 at 18:30