Message Boards

WOLFRAM COMMUNITY

4768 Views

5 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Mathematica Wolfram Language Statistics and Probability

Calculating volume of a nonparametric probability density

Jim Baldwin

Jim Baldwin, Retired

Posted 10 years ago

I've fit a nonparametric density in the following manner: x = Table[{RandomReal[], RandomReal[], RandomReal[]}, {i, 20}]; d = SmoothKernelDistribution[x]; RegionPlot3D[ PDF[d, {x1, x2, x3}] > 0.5, {x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5, 1.5}, AxesLabel -> {"X1", "X2", "X3"}, SphericalRegion -> True, Mesh -> None] Now I'd like to calculate the volume of a region with density greater than, say, z. I've tried the following: z = 0.5; Volume[ImplicitRegion[ PDF[d, {x1, x2, x3}] > z, {{x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5, 1.5}}]] But I get Volume::nmet: Unable to compute the volume of region ImplicitRegion[(\[Piecewise] 0 Less[<<2>>]\|\|Less[<<2>>]\|\|Less[<<2>>]\|\|Greater[<<2>>]\|\|Greater[<<2>>]\|\|Greater[<<2>>] InterpolatingFunction[{{-0.598479,1.48061},{-0.776221,1.78808},{-0.660492,1.69809}},{4,7,0,{32,32,32},{2,2,2},0,0,0,0,Automatic,{},{},False},{<<1>>},{Developer`PackedArrayForm,{<<1>>},{1.0682910^-18,7.0788610^-17,<<47>>,0.,<<32718>>}},{Automatic,Automatic,Automatic}][<<1>>] <<4>> )>0.5&&-0.5<=x1<=1.5&&-0.5<=x2<=1.5&&-0.5<=x3<=1.5,{x1,x2,x3}]. >> Any suggestions would be greatly appreciated.

I've fit a nonparametric density in the following manner:

x = Table[{RandomReal[], RandomReal[], RandomReal[]}, {i, 20}];
d = SmoothKernelDistribution[x];
RegionPlot3D[
 PDF[d, {x1, x2, x3}] > 0.5, {x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5,
   1.5}, AxesLabel -> {"X1", "X2", "X3"}, SphericalRegion -> True, 
 Mesh -> None]

Now I'd like to calculate the volume of a region with density greater than, say, z. I've tried the following:

z = 0.5;
Volume[ImplicitRegion[
  PDF[d, {x1, x2, x3}] > 
   z, {{x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5, 1.5}}]]

But I get

Volume::nmet: Unable to compute the volume of region ImplicitRegion[(\[Piecewise]   0 Less[<<2>>]||Less[<<2>>]||Less[<<2>>]||Greater[<<2>>]||Greater[<<2>>]||Greater[<<2>>]
InterpolatingFunction[{{-0.598479,1.48061},{-0.776221,1.78808},{-0.660492,1.69809}},{4,7,0,{32,32,32},{2,2,2},0,0,0,0,Automatic,{},{},False},{<<1>>},{Developer`PackedArrayForm,{<<1>>},{1.06829*10^-18,7.07886*10^-17,<<47>>,0.,<<32718>>}},{Automatic,Automatic,Automatic}][<<1>>]    <<4>>   
)>0.5&&-0.5<=x1<=1.5&&-0.5<=x2<=1.5&&-0.5<=x3<=1.5,{x1,x2,x3}]. >>

Any suggestions would be greatly appreciated.

POSTED BY: Jim Baldwin

5 Replies

Sort By:

Jim Baldwin

Jim Baldwin, Retired

Posted 10 years ago

I've done more experimenting and the following sort of works: x = Table[{RandomReal[], RandomReal[], RandomReal[]}, {i, 20}]; d = SmoothKernelDistribution[x]; z = 0.4; RegionPlot3D[ PDF[d, {x1, x2, x3}] > z, {x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5, 1.5}, AxesLabel -> {"X1", "X2", "X3"}, SphericalRegion -> True, Mesh -> None] NProbability[PDF[d, {x1, x2, x3}] > z, {x1, x2, x3} \[Distributed] d, WorkingPrecision -> 30] But I get a warning about the NProbability result not being very accurate/precise: NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >> NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.63465481726965450302939589136`30. and 0.0301404857484256325366771310155`30. for the integral and error estimates. >> Is there a better choice for some option to reduce the error estimate?

I've done more experimenting and the following sort of works:

x = Table[{RandomReal[], RandomReal[], RandomReal[]}, {i, 20}];
d = SmoothKernelDistribution[x];
z = 0.4;
RegionPlot3D[
 PDF[d, {x1, x2, x3}] > z, {x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5, 
  1.5}, AxesLabel -> {"X1", "X2", "X3"}, SphericalRegion -> True, 
 Mesh -> None]
NProbability[PDF[d, {x1, x2, x3}] > z, {x1, x2, x3} \[Distributed] d, 
 WorkingPrecision -> 30]

But I get a warning about the NProbability result not being very accurate/precise:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.63465481726965450302939589136`30. and 0.0301404857484256325366771310155`30. for the integral and error estimates. >>

Is there a better choice for some option to reduce the error estimate?

POSTED BY: Jim Baldwin

Ilian Gachevski

Ilian Gachevski, Wolfram Research

Posted 10 years ago

Perhaps In[8]:= Volume[DiscretizeRegion[ImplicitRegion[PDF[d, {x1, x2, x3}] > z, {{x1, -.5, 1.5}, {x2, -.5, 1.5}, {x3, -.5, 1.5}}]]] Out[8]= 0.870019

POSTED BY: Ilian Gachevski

Jim Baldwin

Jim Baldwin, Retired

Posted 10 years ago

Thanks. I tried your suggestion but the result for z = 0.5 is much higher than the other answer and is very much outside of the estimated error. While that doesn't mean the the use of DiscretizeRegion is wrong, I do get very wrong answers with that when z is very low as in z = 0.01. For very low values of z the result needs to be just under 1.0 but the approach you suggest gives me values bigger than 2.4. (And, of course, the true result can't be bigger than 1.0.) I've also tried NProbability[PDF[d, {x1, x2, x3}] > z, {x1, x2, x3} \[Distributed] d, WorkingPrecision -> 30, Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 10000, Method -> "GaussKronrodRule"}] but get the same answer without the Method option.

POSTED BY: Jim Baldwin

Ilian Gachevski

Ilian Gachevski, Wolfram Research

Posted 10 years ago

Yes, it seems the DiscretizeRegion approach is too crude. This is somewhat more plausible: In[15]:= NProbability[PDF[d, {x1, x2, x3}] > 0.01, {x1, x2, x3} \[Distributed] d, Method -> {"NIntegrate", Method -> "LocalAdaptive"}] Out[15]= 0.995179

POSTED BY: Ilian Gachevski

Jim Baldwin

Jim Baldwin, Retired

Posted 10 years ago

Thank you! That gives a much more reasonable answer.

POSTED BY: Jim Baldwin

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback