混合ガウス分布とEMアルゴリズムでの手書き文字(MNISTデータセット)のクラスタリングを実行しようとしました。
Jupyter notebookでpythonを書いたところ、エラーメッセージでreturn np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs.T).T * (2 * np.pi) ** self.ndim)
の部分で(34, 'Result too large')と表示され、オーバーフローによるエラーとわかったのですが、対処方法がわかりません。
よろしくお願いいたします。
以下コードです。
#MNISTデータセットをダウンロード import numpy as np from sklearn import datasets from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split import matplotlib.pyplot as plt from sklearn.datasets import fetch_mldata mnist = fetch_mldata('MNIST original') #MNISTの1~255の数字を1にする _, mnist.data, _, y = train_test_split(mnist.data, y,test_size=0.29) mnist_np = np.array(mnist.data) x_0or1 = np.where(mnist_np >= 1,1, 0) #アルゴリズムの設定 class GaussianMixture(object): def __init__(self,n_component): #ガウス分布の個数 self.n_component = n_component #EMアルゴリズムを用いた最尤 def fit(self,X,iter_max=10): #データの次元 self.ndim = np.size(X,1) #混合係数の初期化 self.weights = np.ones(self.n_component) / self.n_component #平均の初期化 self.means = np.random.uniform(X.min(),X.max(),(self.ndim,self.n_component)) #共分散行列の初期化 self.covs = np.repeat(10*np.eye(self.ndim),self.n_component).reshape(self.ndim,self.ndim,self.n_component) #EステップのMステップを繰り返す for i in range(iter_max): params = np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel())) # Eステップ、負担率を計算 resps = self.expectation(X) # Mステップ、パラメータを更新 self.maximization(X, resps) # パラメータが収束したかを確認 if np.allclose(params, np.hstack((self.weights.ravel(), self.means.ravel(), self.covs.ravel()))): break else: print("parameters may not have converged") # ガウス関数 def gauss(self, X): precisions = np.linalg.inv(self.covs.T).T diffs = X[:, :, None] - self.means assert diffs.shape == (len(X), self.ndim, self.n_component) exponents = np.sum(np.einsum('nik,ijk->njk', diffs, precisions) * diffs, axis=1) assert exponents.shape == (len(X), self.n_component) return np.exp(-0.5 * exponents) / np.sqrt(np.linalg.det(self.covs.T).T * (2 * np.pi) ** self.ndim) # Eステップ def expectation(self, X): # PRML式(9.23) resps = self.weights * self.gauss(X) resps /= resps.sum(axis=-1, keepdims=True) return resps # Mステップ def maximization(self, X, resps): # PRML式(9.27) Nk = np.sum(resps, axis=0) # PRML式(9.26) self.weights = Nk / len(X) # PRML式(9.24) self.means = X.T.dot(resps) / Nk diffs = X[:, :, None] - self.means # PRML式(9.25) self.covs = np.einsum('nik,njk->ijk', diffs, diffs * np.expand_dims(resps, 1)) / Nk # 確率分布p(x)を計算 def predict_proba(self, X): # PRML式(9.7) gauss = self.weights * self.gauss(X) return np.sum(gauss, axis=-1) # クラスタリング def classify(self, X): joint_prob = self.weights * self.gauss(X) return np.argmax(joint_prob, axis=1) # MNISTデータセットを実装 def main(): X = x_0or1 n_component=10 model = GaussianMixture(n_component) model.fit(X, iter_max=300) labels = model.classify(X) x_test, y_test = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100)) X_test = np.array([x_test, y_test]).reshape(2, -1).transpose() probs = model.predict_proba(X_test) Probs = probs.reshape(100, 100) cmap = plt.get_cmap("tab10") plt.scatter(X[:, 0], X[:, 1], c=[cmap[int(label)] for label in labels]) plt.contour(x_test, y_test, Probs) plt.show() if __name__ == '__main__': main()
あなたの回答
tips
プレビュー