pca.py 文件源码

def _em(self, X):
        # Constants
        n_samples, n_features = X.shape
        n_components = self.n_components
        max_iter = self.max_iter
        # tol = self.tol

        mu = X.mean(axis=0)
        X_centered = X - sp.atleast_2d(mu)

        # Initialize parameters
        latent_mean = 0
        sigma2 = 1
        weight = sprd.randn(n_features, n_components)

        # Main loop of EM algorithm
        for i in range(max_iter):
            # E step
            M = sp.dot(weight.T, weight) + sigma2 * sp.eye(n_components)
            inv_M = spla.inv(M)
            latent_mean = sp.dot(inv_M, sp.dot(weight.T, X_centered.T)).T

            # M step
            expectation_zzT = n_samples * sigma2 * inv_M + sp.dot(latent_mean.T, latent_mean)

            # Re-estimate W
            weight = sp.dot(sp.dot(X_centered.T, latent_mean), spla.inv(expectation_zzT))
            weight2 = sp.dot(weight.T, weight)

            # Re-estimate \sigma^2
            sigma2 = ((spla.norm(X_centered)**2 -
                       2 * sp.dot(latent_mean.ravel(), sp.dot(X_centered, weight).ravel()) +
                       sp.trace(sp.dot(expectation_zzT, weight2))) /
                      (n_samples * n_features))

        self.predict_mean = mu
        self.predict_cov = sp.dot(weight, weight.T) + sigma2 * sp.eye(n_features)
        self.latent_mean = latent_mean
        self.latent_cov = sigma2 * inv_M
        self.sigma2 = sigma2
        self.weight = weight
        self.inv_M = inv_M

        return self.latent_mean