2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ I'm missing the values of s1,s2,s3,h1,h2,h3 $\endgroup$ Commented Dec 20, 2023 at 13:59
  • $\begingroup$ 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] $\endgroup$ Commented Dec 21, 2023 at 3:50

1 Answer 1

3
$\begingroup$

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", ""}]

Figure 1

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", ""}]

Figure 2

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Commented Dec 21, 2023 at 3:26
  • $\begingroup$ 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. $\endgroup$ Commented Dec 21, 2023 at 4:39
  • $\begingroup$ @CarsonNewman Please, see Update 1 to my post. $\endgroup$ Commented Dec 21, 2023 at 8:04
  • $\begingroup$ Fantastic! And to make the orange line sharper I just have to increase the interpolation and Precision goals? $\endgroup$ Commented Dec 21, 2023 at 18:08
  • $\begingroup$ @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} . $\endgroup$ Commented Dec 21, 2023 at 18:30

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.