在 MATLAB 中使用双积分 (dblquad) 查找重叠的正态分布。奇怪的行为答案

【问题标题】：Finding normal distributions overlap using double Integral (dblquad) in MATLAB. Strange behaviour在 MATLAB 中使用双积分 (dblquad) 查找重叠的正态分布。奇怪的行为
【发布时间】：2013-04-16 15:36:24
【问题描述】：

我正在使用以下函数计算两个正态双变量分布的重叠

function [ oa ] = bivariate_overlap_integral(mu_x1,mu_y1,mu_x2,mu_y2)
    %calculating pdf. Using x as vector because of MATLAB requirements for integration
    bpdf_vec1=@(x,y,mu_x,mu_y)(exp(-((x-mu_x).^2)./2.-((y-mu_y)^2)/2)./(2*pi));

    %calcualting overlap of two distributions at the point x,y
    overlap_point = @(x,y) min(bpdf_vec1(x,y,mu_x1,mu_y1),bpdf_vec1(x,y,mu_x2,mu_y2));

    %calculating overall overlap area
    oa=dblquad(overlap_point,-100,100,-100,100);

您可以看到这涉及从最小为 2 pdf- s 由点 x,y 的两个分布的函数 bpdf_vec1 给出。

现在，PDF 永远不会为 0，所以我希望区间的面积越大，最终结果就会越大，显然在某个点之后差异可以忽略不计。但是，有时，当我减小间隔的大小时，结果似乎会增加。例如：

>> mu_x1=0;mu_y1=0;mu_x2=5;mu_y2=0;
>> bpdf_vec1=@(x,y,mu_x,mu_y)(exp(-((x-mu_x).^2)./2.-((y-mu_y)^2)/2)./(2*pi));
>>  overlap_point = @(x,y) min(bpdf_vec1(x,y,mu_x1,mu_y1),bpdf_vec1(x,y,mu_x2,mu_y2));
>> dblquad(overlap_point,-10,10,-10,10)
ans =
    0.0124
>> dblquad(overlap_point,-100,100,-100,100)
ans =
    1.4976e-005  -----> strange, as theoretically cannot be smaller then the first answer
>> dblquad(overlap_point,-3,3,-3,3)
ans =
    0.0110  -----> makes sense that the result is less than the first answer as the           
interval is decreased

在这里我们可以检查在区间的边界点重叠是否（接近）0。

>> overlap_point (100,100)
ans =
 0
>> overlap_point (-100,100)
ans =
 0
>> overlap_point (-100,-100)
ans =
 0
>> overlap_point (100,-100)
ans =
 0

这可能与 dblquad 的实现有关，还是我在某处犯了错误？我使用 MATLAB R2011a。

谢谢

【问题讨论】：

我检查了你的代码，没有发现任何错误，所以我认为这是 matlab 的精度问题（精度与限制的范围有关）。我会用八度音阶测试你的代码，看看它是否有同样的问题。顺便说一句，您可以尝试使用可选的精度参数，例如 dblquad(-100,100,-100,100,1e-8)。
有趣的是，在 gnu-Octave 中运行完全相同的代码时并没有出现问题。在默认精度下，[-10,10,-10,10] 上的积分结果与 [-100,100,-100,100] 上的积分结果几乎没有区别。两个结果都约为 0.01241933 btw。
你能不能通过逐渐增加限制来测试，看看它首先开始失败的地方。像 [-30,30,-30,30] 然后 [-40,40,-40,40] 等等，直到结果开始发散。然后尝试添加可选的精度参数 1e-7 或 1e-8 等，看看是否能解决问题。顺便提一句。我的 Octave 版本是 3.6.1
感谢 cmets。结果从 [-10,10,-10,10] >> dblquad(overlap_point,-9,9,-9,9) ans = 0.012434302954144 >> dblquad(overlap_point,-10,10,-10,10) ans = 0.012426083045525 开始递减。指定（一个相当大的）精度参数似乎确实有积极的效果，但并不能完全解决问题，例如：>> dblquad(-9,9,-9,9,1e-10) ans = 0.012419331486452、dblquad(-10,10,-10,10,1e-10) ans = 0.012419330851169、dblquad(-100,100,-100,100,1e-10) ans = 0.012415587111195 结果彼此更接近，但仍然从[-10,10,-10,10]。
Octave 中的这些值是略微减少还是增加？

标签： matlab statistics normal-distribution integral

【解决方案1】：

恭喜！您因成为第 12 百万个提出这个问题的人而获奖。 :) 我想说的是，这是一个每个人一开始都会遇到的问题。老实说，这个问题被问了一遍又一遍，所以说真的，这个问题应该被标记为 dup。

当从足够远的地方观察时，这些事情发生的情况是双变量法线本质上是一个 delta 函数。而且您不需要真正将该区域分散得太远，因为正常密度会迅速下降。在您尝试整合的大部分领域，它基本上为零，至少在所采用的容差范围内。

因此，如果求积恰好碰到有质量的区域附近的一些样本点，那么您可能会得到对积分的一些实际估计。但是，如果该工具看到的所有数字都是在整个域上基本上为零的数字，则它得出的结论是积分为零。请记住，自适应集成工具并非无所不知。他们对您的功能一无所知。对他们来说，这是一个黑匣子。这些不是象征性的工具。

顺便说一句，这不是我希望看到 Octave 与 MATLAB 始终不同的东西。这只是自适应积分器的问题，以及它选择在哪里设置其采样点。

【讨论】：

所以，本质上唯一的解决方案是让积分的界限更小？
我并不是真的在抱怨，只是说你不是唯一一个被这个问题绊倒的人。我想我们都曾经这样做过。谨慎、明智地选择边界将非常有帮助。然而，太小了，你会失去一些质量。一些正交工具允许您指定兴趣点，因此如果您可以指示法线的峰值，那么您将遥遥领先。我不记得 dblquad 有这种能力。我已经编写了二维自适应集成工具，而自动完成这个问题绝非易事。

【解决方案2】：

好的，这是我的八度音阶结果。

>fomat long
>z10 = dblquad(overlap_point,-10,10,-10,10)
z10 =  0.0124193303284589
>z100 = dblquad(overlap_point,-100,100,-100,100)
z100 =  0.0124193303245106
>z100 - z10
ans = -3.94835379669001e-012

>z10a = dblquad(overlap_point,-10,10,-10,10,1e-8)
z10a =  0.0124193306514996
>z100a = dblquad(overlap_point,-100,100,-100,100,1e-8)
z100a =  0.0124193306527155
>z100a-z10a
ans = 1.21588676627038e-012

顺便说一句。我以前用数值解法就注意到了这种类型的问题。有时您进行的更改会提高结果的准确性（在这种情况下，通过使您的限制更接近完整平面的理想情况），但您会得到相反的效果，结果会变得不那么准确。在这种情况下发生的情况是，通过“扩大”到 -100..100，您将焦点从函数中真正重要的动作发生的地方转移，这接近原点。在某些时候，您正在使用的 dblquad 的实现必须随着您增加限制而开始增加样本间距离，然后它开始丢失靠近原点的一些重要内容。

也许运行更高版本的 matlab 的人可以检查一下，看看它是否得到了改进。

【讨论】：