# Calculating volume of a nonparametric probability density

Posted 9 years ago
4316 Views
|
5 Replies
|
0 Total Likes
|
 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>>},{DeveloperPackedArrayForm,{<<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.
5 Replies
Sort By:
Posted 9 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 9 years ago
 Thank you! That gives a much more reasonable answer.
Posted 9 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 9 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 9 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.6346548172696545030293958913630. and 0.030140485748425632536677131015530. for the integral and error estimates. >> `Is there a better choice for some option to reduce the error estimate?