Message Boards Message Boards

0
|
5013 Views
|
5 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Calculating volume of a nonparametric probability density

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.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:

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

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
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?
POSTED BY: Jim Baldwin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract