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}
, Method -> {"LocalAdaptive"}
]
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!