Re: Need help NIntegrating stepwise 3D probability density functions

• To: mathgroup at smc.vnet.net
• Subject: [mg114098] Re: Need help NIntegrating stepwise 3D probability density functions
• From: Jason Neuswanger <j.neuswanger at sfos.uaf.edu>
• Date: Tue, 23 Nov 2010 06:01:24 -0500 (EST)

```Barrie,

Thanks for your help!  Unfortunately, I apparently picked a too-simple
test case before, because the original problem remains when using
KernelMixtureDistribution or SmoothKernelDistribution.  Here's an
example with actual data:

dataPoints = {{-0.053, -0.006, 0.077}, {-0.062, 0.037,
0.068}, {-0.201, 0.141, 0.104}, {-0.164, -0.009, 0.107}, {-0.221,
0.135, 0.147}, {-0.091, 0.004, 0.072}, {-0.063, 0.026,
0.008}, {-0.199, 0.069, 0.153}, {-0.144, -0.046, 0.054}, {-0.128,
0.01, 0.038}, {-0.21, 0.097, 0.172}, {-0.147, 0.147,
0.12}, {-0.098, 0.069, 0.06}, {-0.138, 0.06, 0.05}, {-0.083,
0.061, 0.054}, {-0.083, 0., 0.09}, {-0.187, 0.041,
0.141}, {-0.093, 0.008, 0.061}, {-0.116, -0.007, 0.159}, {-0.06,
0.041, 0.051}, {-0.172, 0.212, 0.164}, {-0.076, 0.029,
0.065}, {-0.079, 0.008, 0.085}, {-0.11, 0.025, 0.077}, {-0.12,
0.085, 0.085}, {-0.244, 0.261, 0.185}, {-0.08, 0.028,
0.087}, {-0.18, 0.191, 0.071}, {-0.13, 0.028,
0.062}, {-0.082, -0.022, 0.064}, {-0.081, 0.009, 0.068}, {-0.231,
0.256, 0.29}, {0.056, 0.092, 0.}, {-0.083, 0.009,
0.045}, {-0.053, -0.02, 0.102}, {-0.117, 0.018, 0.006}, {-0.051,
0.035, 0.036}, {-0.048, 0.07, 0.05}, {-0.177, 0.16,
0.192}, {-0.111, 0.075, 0.048}, {0.01, 0.008, -0.002}, {-0.196,
0.068, 0.137}, {-0.071, 0.01, 0.074}, {-0.079, 0.067,
0.062}, {-0.106, 0.039, 0.04}, {-0.033, 0.055, 0.055}, {-0.249,
0.258, 0.2}, {-0.23, 0.159, 0.171}, {-0.078, 0.023,
0.047}, {-0.099, 0.025, 0.065}, {-0.129, 0.004, 0.094}, {-0.069,
0.006, 0.101}, {-0.098, 0.067, 0.058}, {-0.101, 0.052,
0.059}, {-0.088, 0.031, 0.073}};
kernelTestDistribution = KernelMixtureDistribution[dataPoints, "Scott"];

thresholdDensity = 100; (* the max probability density in this PDF is
around 700 *)

NIntegrate[
PDF[kernelTestDistribution, {x, y, z}]*
UnitStep[PDF[testDistribution, {x, y, z}] - thresholdDensity]
, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity,
Infinity}
]

This gives the value of the integral as 0 anytime the threshold density
is >0 (all the useful cases).  When the threshold density is <=0 and the
function should integrate to 1, it gives 0.74001.  This value does
change slightly when I use different data, but it's obviously all wrong.

I hope someone can figure this out!

Thanks,
Jason

On 11/22/10 6:21 PM, Barrie Stokes wrote:
> Hi Jason
>
> I admit I tried a few things before I tried the following:
>
> thresholdDensity = 5/100;
> NIntegrate[
>   PDF[ testDistribution, {x, y, z} ]*
>    UnitStep[ PDF[ testDistribution, {x, y, z} ] - thresholdDensity ]
>   , {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity,
>    Infinity} , Method ->  {"LocalAdaptive"}
>   ]
>
> thresholdDensity = 4/100;
> NIntegrate[
>   PDF[ testDistribution, {x, y, z} ]*
>    UnitStep[ PDF[ testDistribution, {x, y, z} ] - thresholdDensity ]
>   , {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity,
>    Infinity} , Method ->  {"LocalAdaptive"}
>   ]
>
> thresholdDensity = 6/100;
> NIntegrate[
>   PDF[ testDistribution, {x, y, z} ]*
>    UnitStep[ PDF[ testDistribution, {x, y, z} ] - thresholdDensity ]
>   , {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity,
>    Infinity} , Method ->  {"LocalAdaptive"}
>   ]
>
> thresholdDensity = 0/100;
> NIntegrate[
>   PDF[ testDistribution, {x, y, z} ]*
>    UnitStep[ PDF[ testDistribution, {x, y, z} ] - thresholdDensity ]
>   , {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity,
>    Infinity} , Method ->  {"LocalAdaptive"}
>   ]
>
> The answers look like they are behaving correctly. I too read the documentation in UnitStep about "SymbolicProcessing" ->  0, and I only tried removing it on a whim. Go figure. Maybe Daniel Lichtblau can explain it for us?
>
> Cheers
>
> Barrie
>
>
>
>>>> On 22/11/2010 at 11:40 pm, in message<201011221240.HAA06705 at smc.vnet.net>,
> Jason Neuswanger<jrneuswanger at alaska.edu>  wrote:
>> 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!

```

• Prev by Date: Re: Need help NIntegrating stepwise 3D probability density functions