Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Need help NIntegrating stepwise 3D probability density functions

• To: mathgroup at smc.vnet.net
• Subject: [mg114102] Re: Need help NIntegrating stepwise 3D probability density functions
• From: Andrew Moylan <amoylan at wolfram.com>
• Date: Tue, 23 Nov 2010 06:02:15 -0500 (EST)

Hi Jason,

The integrand is only non-zero in a small region. Visualize the region:

testDistribution = 
  MultinormalDistribution[{1, 2, 
    4}, {{2, -1/4, 1/3}, {-1/4, 2/3, 1/5}, {1/3, 1/5, 1/2}}];
thresholdDensity = 0.05;

f[x_, y_, z_] = 
  PDF[testDistribution, {x, y, z}]*
   UnitStep[PDF[testDistribution, {x, y, z}] - thresholdDensity];

region = RegionPlot3D[
  f[x, y, z] != 0, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}, 
  PlotPoints -> 35, PlotStyle -> Opacity[0.5], Mesh -> None]


With the Method option you used, NIntegrate's sampling points all land on places where the integrand is zero. Visualize the sampling points:

pts = Last[
   Last[Reap[
     NIntegrate[
      h[x, y, z] f[x, y, z], {x, -Infinity, Infinity}, {y, -Infinity, 
       Infinity}, {z, -Infinity, Infinity}, 
      Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
      EvaluationMonitor :> Sow[{x, y, z}]]]]];

Show[region, Graphics3D[{Red, Point[pts]}]]


The default NIntegrate methods use symbolic methods to first isolate the non-zero region. This works fine: 

NIntegrate[
 f[x, y, z], {x, -Infinity, Infinity}, {y, -Infinity, 
  Infinity}, {z, -Infinity, Infinity}]


Visualize these sampling points:

pts2 = Last[
   Last[Reap[
     Quiet[NIntegrate[
       f[x, y, z], {x, -Infinity, Infinity}, {y, -Infinity, 
        Infinity}, {z, -Infinity, Infinity}, MaxRecursion -> 0, 
       EvaluationMonitor :> Sow[{x, y, z}]]]]]];

Show[region, Graphics3D[{Red, Point[pts2]}]]


If you like you can use this internal function to see the region(s) into which NIntegrate subdivides an integral:

NIntegrate`PiecewiseNIntegrateMultipleRanges[
 f[x, y, z], {{x, -Infinity, Infinity}, {y, -Infinity, 
   Infinity}, {z, -Infinity, Infinity}}, {}, {}]


It is the option "SymbolicProcessing" -> 0 that disabled NIntegrate's detection of the non-zero region. Any reason you need that option?

I also recommend "GlobalAdaptive" (which is the default) instead of "LocalAdaptive".

If symbolic processing is too expensive, you can turn it off---but then you need some other way to tell NIntegrate where the small non-zero region is. For example you could use:

NIntegrate[
 f[x, y, z], {x, -Infinity, 3, Infinity}, {y, -Infinity, 3, 
  Infinity}, {z, -Infinity, 3, Infinity}, 
 Method -> {Automatic, "SymbolicProcessing" -> 0}]

You will get messages and slow slow convergence though, since the integrand is non-smooth over the specified region.

Finally, you can try Monte Carlo integration:

NIntegrate[
 f[x, y, z], {x, -Infinity, 3, Infinity}, {y, -Infinity, 3, 
  Infinity}, {z, -Infinity, 3, Infinity}, Method -> "MonteCarlo", 
 MaxPoints -> 10^6]

Note that MonteCarlo integration uses PrecisionGoal -> 2 instead of PrecisionGoal -> 6 by default.

Regards,

Andrew Moylan
Wolfram Research






----- Original Message -----
From: "Jason Neuswanger" <jrneuswanger at alaska.edu>
To: mathgroup at smc.vnet.net
Sent: Monday, November 22, 2010 11:40:36 PM
Subject: [mg114102] [mg114071] Need help NIntegrating stepwise 3D probability density	functions

I need to integrate 3D probability density functions over the region
in which the probability density exceeds some threshold (the end goal
being to calculate isopleths in 3D animal territories fit by kernel
distributions).

The simplest test case I could come up with is this:

testDistribution =
  MultinormalDistribution[{1, 2,
    4}, {{2, -1/4, 1/3}, {-1/4, 2/3, 1/5}, {1/3, 1/5, 1/2}}];

thresholdDensity = 0.05;

NIntegrate[
 PDF[testDistribution, {x, y, z}]*
  UnitStep[PDF[testDistribution, {x, y, z}] - thresholdDensity]
 , {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity,
  Infinity}
 , Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}
 ]

I've tried quite a few things and just can't get remotely sensible
results.  The analogous case in 1 dimension works just fine, shown
below:

testDistribution2 = NormalDistribution[];

thresholdDensity2 = 0.35;

NIntegrate[
 PDF[testDistribution2, {x}]*
  UnitStep[PDF[testDistribution2, {x}] - thresholdDensity2]
 , {x, -Infinity, Infinity}
 , Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}
 ]

(* check answer using symbolic methods *)

thresholdRoot =
  x /. FindRoot[
    PDF[testDistribution2, {x}] - thresholdDensity2, {x, 1}];
Integrate[
 PDF[testDistribution2, {x}], {x, -thresholdRoot, thresholdRoot}]

Any help here would be VERY much appreciated!