寻找高斯分布的交点

Question

如何在 MATLAB 中找到高斯混合分布相交的所有点？

Answer 1

您的问题的一般类别是找到两条曲线的交点，这是一项易于管理但并非微不足道的任务（最难的部分是确保您抓住所有交点）。

但您的问题非常具体：您正在寻找两个高斯的交集。这非常好：我们都为您的函数提供了一个解析公式，并且除非您的某些参数相同，否则可以保证恰好有两个交集。*

假设您的分布通过mu1、mu2 和规模sigma1、sigma2 来表征。然后你在x 位置的高斯函数由函数定义

1/sqrt(2*pi*sigma^2) * exp(-(x-mu)^2/2/sigma^2)

事实证明，我们可以在纸上完全解出x 的这个方程：

1/sqrt(2*pi*sigma1^2)*exp(-(x-mu1)^2/2/sigma1^2) == 1/sqrt(2*pi*sigma2^2)*exp(-(x-mu2)^2/2/sigma2^2)
sigma2/sigma1 == exp((x-mu1)^2/2/sigma1^2) * exp(-(x-mu2)^2/2/sigma2^2)
log(sigma2/sigma1) == (x-mu1)^2/2/sigma1^2) - (x-mu2)^2/2/sigma2^2

导致抛物线方程ax^2 + bx + c == 0 where

a = 1/(2*sigma1^2) - 1/(2*sigma2^2);
b = mu2/(sigma2^2) - mu1/(sigma1^2);
c = mu1^2/(2*sigma1^2) - mu2^2/(2*sigma2^2) - log(sigma2/sigma1);

很容易证明D = b^2 - 4 a c判别式是非负的，所以当参数不退化时，方程确实有两个实根。所以这两个交点，按照上面的定义，

D = b^2 - 4 * a * c;
x1 = (-b + sqrt(D))/(2*a);
x2 = (-b - sqrt(D))/(2*a);

使用两个具有伪随机参数的高斯函数：

% define parameters and Gaussian
mu1=1; sigma1=3; mu2=2; sigma2=4;

% intersections
a = 1/(2*sigma1^2) - 1/(2*sigma2^2);
b = mu2/(sigma2^2) - mu1/(sigma1^2);
c = mu1^2/(2*sigma1^2) - mu2^2/(2*sigma2^2) - log(sigma2/sigma1)
D = b^2 - 4 * a * c;
x1 = (-b + sqrt(D))/(2*a);
x2 = (-b - sqrt(D))/(2*a);

只是为了证明上面的交点是正确的：

>> f = @(x,mu,sigma) 1/sqrt(2*pi*sigma^2) * exp(-(x-mu).^2/2/sigma^2);
>> f(x1,mu1,sigma1) - f(x1,mu2,sigma2)

ans =

   2.7756e-17

>> f(x2,mu1,sigma1) - f(x2,mu2,sigma2)

ans =

   1.0408e-17

以上意味着x1和x2点的两个高斯值在机器精度内是相等的，这与任何数字答案一样好。

---

*我最初声称只要高斯不完全相同，我们总是有两个交点。显然，如果均值和方差都相同，我们就会有两条退化曲线，并且交点变得毫无意义。但就像Cris Luengo pointed out in a comment 一样，可能只有一个交点：当方差相同且均值不同时（即，我们有两条形状完全相同的曲线沿x 移动）。在这种情况下a=0，因此对应的方程是b*x + c == 0，给我们x0 = -c/b作为交集。所以一个更准确（但有点伪代码）的答案（给定a、b 和c）是

if a == 0 % or allow some tolerance... <=> sigma1 == sigma2
   if b == 0 % or allow some tolerance... <=> mu1 == mu2
      % degenerate curves: a == b == c == 0, f1(x)==f2(x) for all x
      disp('curves are degenerate...')
   else
      % single intersection: mu1 ~= mu2
      x1 = -c/b;
   end
else
   % two intersections; both parameters are different
   D = b^2 - 4 * a * c;
   x1 = (-b + sqrt(D))/(2*a);
   x2 = (-b - sqrt(D))/(2*a);