轮廓抛出一个错误试图绘制多元高斯

Question

我编写了一个小脚本来模拟EM算法并可视化其迭代步骤。然而，在第5次迭代之后，当试图绘制更新的估计双变量高斯分布时，它会停止。

我怀疑我的协方差矩阵有什么可疑之处，但我并不确定。如果我评论等高线图，脚本运行良好，并且按照它应有的工作方式工作（但当然，遵循估计分布的演变将会很好）。任何帮助，将不胜感激。

import numpy as np 
import scipy.stats as sp 
import matplotlib.pyplot as plt 
from matplotlib.mlab import bivariate_normal 


def expectationMaximization(): 
    # define multivariate gaussian distributions and generate observations 
    u1 = [-1.5, -1.5] 
    cov1 = [[0.2, 0.4], 
      [0, 0.1]] 

    u2 = [1, 1] 
    cov2 = [[0.3, 0.4], 
      [0, 0.3]] 

    samples = 1000 

    x1, y1 = np.random.multivariate_normal(u1, cov1, samples // 2).T 
    x2, y2 = np.random.multivariate_normal(u2, cov2, samples // 2).T 

    x = np.concatenate([x1, x2]) 
    y = np.concatenate([y1, y2]) 
    points = np.concatenate([np.column_stack((x1, y1)), 
          np.column_stack((x2, y2))]) 

    # initialization of classifier models 
    uk1 = np.array([-1.5, 1]) 
    covk1 = np.array([[1, 0], [0, 1]]) 

    uk2 = np.array([1.5, -1]) 
    covk2 = np.array([[1, 0], [0, 1]]) 

    w = np.array([1., 1.]) 
    gamma = np.zeros((2, samples)) 

    # sim loop 
    for idx in range(9): 
     ########################################################## 
     # expectation # 
     ########################################################## 
     # update gamma 
     gamma[0] = (w[0] * sp.multivariate_normal.pdf(points, uk1, covk1)/
        (w[0] * sp.multivariate_normal.pdf(points, uk1, covk1) + 
        w[1] * sp.multivariate_normal.pdf(points, uk2, covk2))) 
     gamma[1] = (w[1] * sp.multivariate_normal.pdf(points, uk2, covk2)/
        (w[0] * sp.multivariate_normal.pdf(points, uk1, covk1) + 
        w[1] * sp.multivariate_normal.pdf(points, uk2, covk2))) 

     ########################################################## 
     # plot # 
     ########################################################## 
     plt.subplot(3, 3, idx + 1) 
     plt.title('Iteration {}'.format(idx + 1)) 
     axes = plt.gca() 
     axes.set_xlim([-3, 3]) 
     axes.set_ylim([-3, 3]) 

     # setup grid for bivariate gaussian plot (only needed once) 
     if idx < 1: 
      xmin, xmax = axes.get_xlim() 
      ymin, ymax = axes.get_ylim() 
      delta = 0.1 
      xticks = np.arange(xmin, xmax, delta) 
      yticks = np.arange(ymin, ymax, delta) 
      xmesh, ymesh = np.meshgrid(xticks, yticks) 

     # update mesh values 
     z1 = bivariate_normal(xmesh, ymesh, covk1[0, 0], covk1[1, 1], 
           uk1[0], uk1[1], covk1[0, 1]) 
     z2 = bivariate_normal(xmesh, ymesh, covk2[0, 0], covk2[1, 1], 
           uk2[0], uk2[1], covk2[0, 1]) 
     z = (z1 - z2) * 10 

     # plot pdf map and sample points 
     plt.contour(xmesh, ymesh, z) 
     plt.scatter(x, y, c=(gamma[0] - gamma[1]) * 10) 
     plt.grid(True) 

     ########################################################## 
     # maximization # 
     ########################################################## 
     # update means 
     uk1[0] = sum(gamma[0] * x)/sum(gamma[0]) 
     uk1[1] = sum(gamma[0] * y)/sum(gamma[0]) 

     uk2[0] = sum(gamma[1] * x)/sum(gamma[1]) 
     uk2[1] = sum(gamma[1] * y)/sum(gamma[1]) 

     # update covariance matrices 
     # calc all distances 
     dist1 = points - uk1[None, :] 
     dist2 = points - uk2[None, :] 

     # calc all outer products 
     matrixSchaar1 = np.einsum('...i,...j->...ij', dist1, dist1) 
     matrixSchaar2 = np.einsum('...i,...j->...ij', dist2, dist2) 

     # calculate sum product of matrices and gammas 
     covk1 = ((matrixSchaar1 * gamma[0][:, None, None]).sum(axis=0)/
       sum(gamma[0])) 
     covk2 = ((matrixSchaar2 * gamma[1][:, None, None]).sum(axis=0)/
       sum(gamma[1])) 

     # update w 
     w[0] = sum(gamma[0])/len(gamma[0]) 
     w[1] = sum(gamma[1])/len(gamma[1]) 


def main(): 
    expectationMaximization() 
    plt.show() 

if __name__ == "__main__": 
    main() 

回溯：

/usr/lib64/python3.5/site-packages/matplotlib/mlab.py:1926: RuntimeWarning: invalid value encountered in sqrt 
    denom = 2*np.pi*sigmax*sigmay*np.sqrt(1-rho**2) 
Traceback (most recent call last): 
    File "bsp4.py", line 120, in <module> 
    main() 
    File "bsp4.py", line 115, in main 
    expectationMaximization() 
    File "bsp4.py", line 76, in expectationMaximization 
    plt.contour(xmesh, ymesh, z) 
    File "/usr/lib64/python3.5/site-packages/matplotlib/pyplot.py", line 2766, in contour 
    ret = ax.contour(*args, **kwargs) 
    File "/usr/lib64/python3.5/site-packages/matplotlib/__init__.py", line 1815, in inner 
    return func(ax, *args, **kwargs) 
    File "/usr/lib64/python3.5/site-packages/matplotlib/axes/_axes.py", line 5644, in contour 
    return mcontour.QuadContourSet(self, *args, **kwargs) 
    File "/usr/lib64/python3.5/site-packages/matplotlib/contour.py", line 1424, in __init__ 
    ContourSet.__init__(self, ax, *args, **kwargs) 
    File "/usr/lib64/python3.5/site-packages/matplotlib/contour.py", line 864, in __init__ 
    self._process_levels() 
    File "/usr/lib64/python3.5/site-packages/matplotlib/contour.py", line 1202, in _process_levels 
    self.vmin = np.amin(self.levels) 
    File "/usr/lib64/python3.5/site-packages/numpy/core/fromnumeric.py", line 2359, in amin 
    out=out, keepdims=keepdims) 
    File "/usr/lib64/python3.5/site-packages/numpy/core/_methods.py", line 29, in _amin 
    return umr_minimum(a, axis, None, out, keepdims) 
ValueError: zero-size array to reduction operation minimum which has no identity 

Answer 1

时间感觉像一个笨蛋。

找到了答案here。

总之，bivariate_normal()需要标准偏差值，而不是差异。

纠正我的错误后，一切按预期工作。

Answer 2

开始作为一个评论，但变得太长......

你是对怀疑你的协变矩阵。您的bivariate_normal分布，z1和z2，在第5次或第6次迭代后始终变为nan阵列。但是，我认为这不一定意味着你的代码有缺陷，这可能只是一个不便的数字问题（我对此理论一无所知）。

无论如何，你可以通过设置级别的等高线图，让你的轮廓到i次迭代击穿之前

levels = np.linspace(z.min(), z.max(), 20) 
plt.contour(xmesh, ymesh, z, levels=levels) 

，因为从我可以告诉ValueError关于从轮廓不来能够设定其最小值。