python - EM 算法的缓冲区溢出

Question

我正在尝试为高斯混合实现 EM 算法。我基于以下教程:

http://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html

该算法适用于某些值，但当数字较大时，我通常会在矩阵上得到无穷大。

这是我的代码:

from scipy.stats import multivariate_normal as mvn
import numpy as np
import matplotlib.pyplot as plt

def em_gmm_orig(xs, pis, mus, sigmas, tol=0.01, max_iter=100):

    n, p = xs.shape
    k = len(pis)

    ll_old = 0
    for i in range(max_iter):
        exp_A = []
        exp_B = []
        ll_new = 0

        # E-step
        ws = np.zeros((k, n))
        for j in range(len(mus)):
            for i in range(n):
                ws[j, i] = pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
        ws /= ws.sum(0)

        # M-step
        pis = np.zeros(k)
        for j in range(len(mus)):
            for i in range(n):
                pis[j] += ws[j, i]
        pis /= n

        mus = np.zeros((k, p))
        for j in range(k):
            for i in range(n):
                mus[j] += ws[j, i] * xs[i]
            mus[j] /= ws[j, :].sum()

        sigmas = np.zeros((k, p, p))
        for j in range(k):
            for i in range(n):
                ys = np.reshape(xs[i]- mus[j], (2,1))
                sigmas[j] += ws[j, i] * np.dot(ys, ys.T)
            sigmas[j] /= ws[j,:].sum()

        # update complete log likelihoood
        ll_new = 0.0
        for i in range(n):
            s = 0
            for j in range(k):
                s += pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
            ll_new += np.log(s)

        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new

    return ll_new, pis, mus, sigmas


np.random.seed(123)

# create data set
n = 100
for i in range(0,4):
    num1 = np.random.uniform(-500,500)
    num2 = np.random.uniform(-500,500)
    num3 = np.random.uniform(-500,500)
    num4 = np.random.uniform(-500,500)
_mus = np.array([[num1,num2], [num3,num4]])

_sigmas = np.array([[[3, 0], [0, 0.5]], [[1,0],[0,2]]])
_pis = np.array([0.6, 0.4])
xs = np.concatenate([np.random.multivariate_normal(mu, sigma, int(pi*n))
                    for pi, mu, sigma in zip(_pis, _mus, _sigmas)])

# initial guesses for parameters
num_clusters = 4
pis = np.random.random(num_clusters)
pis /= pis.sum()
mus = np.random.random((num_clusters,2))
print(mus)
sigmas = np.array([np.eye(2)] * num_clusters)

ll1, pis1, mus1, sigmas1 = em_gmm_orig(xs, pis, mus, sigmas)

#Ploy
intervals = 101
ys = np.linspace(-8,8,intervals)
X, Y = np.meshgrid(ys, ys)
_ys = np.vstack([X.ravel(), Y.ravel()]).T

z = np.zeros(len(_ys))
for pi, mu, sigma in zip(pis1, mus1, sigmas1):
    z += pi*mvn(mu, sigma).pdf(_ys)
z = z.reshape((intervals, intervals))

ax = plt.subplot(111)
plt.scatter(xs[:,0], xs[:,1], alpha=0.2)
plt.contour(X, Y, z, N=10)
plt.axis([-8,6,-6,8])
ax.axes.set_aspect('equal')
plt.tight_layout()    

我通常会收到以下错误:

    Traceback (most recent call last):
      File "C:\Usy", line 82, in <module>
        ll1, pis1, mus1, sigmas1 = em_gmm_orig(xs, pis, mus, sigmas)
      File "C:\Usy", line 48, in em_gmm_orig
        s += pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
      File "y", line 354, in __call__
        allow_singular=allow_singular)
...
    ValueError: array must not contain infs or NaNs

有什么办法可以解决这个问题吗？

Answer 1

最有可能的是，您的 sigma 值在 mvn(mus[j], sigmas[j]).pdf(xs[i]) 中为零，这给出了奇点。您可能想要检查您的西格玛值。