contour 在尝试绘制多元高斯时抛出错误
contour throws an error trying to plot multivariate gaussian
我编写了一个小脚本来模拟 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
开始作为评论，但变得太长了...
你怀疑你的协变矩阵是对的。您的 bivariate_normal 分布，z1 和 z2，在第 5 次或第 6 次迭代后始终变为 nan 数组。但是，我认为这并不一定意味着您的代码有缺陷，它可能只是一个不方便的数字问题（我对这里的理论一无所知）。
无论如何，您可以通过设置等高线图的级别
，使等高线在分解前达到第 i 次迭代
levels = np.linspace(z.min(), z.max(), 20)
plt.contour(xmesh, ymesh, z, levels=levels)
据我所知，ValueError 是由于等高线无法设置其最小值所致。
是时候感觉自己像个骗子了。
找到答案here。
简而言之，bivariate_normal() 需要标准偏差值，而不是方差。
改正错误后，一切正常。
我编写了一个小脚本来模拟 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
开始作为评论，但变得太长了...
你怀疑你的协变矩阵是对的。您的 bivariate_normal 分布，z1 和 z2，在第 5 次或第 6 次迭代后始终变为 nan 数组。但是，我认为这并不一定意味着您的代码有缺陷，它可能只是一个不方便的数字问题（我对这里的理论一无所知）。
无论如何，您可以通过设置等高线图的级别
，使等高线在分解前达到第i 次迭代
levels = np.linspace(z.min(), z.max(), 20)
plt.contour(xmesh, ymesh, z, levels=levels)
据我所知，ValueError 是由于等高线无法设置其最小值所致。
是时候感觉自己像个骗子了。
找到答案here。
简而言之，bivariate_normal() 需要标准偏差值，而不是方差。
改正错误后，一切正常。