あとでCとopmempで書きなおす

npgraph.py

# -*- coding: utf-8 -*-
import numpy as np
from scipy.special.basic import digamma
import logging
class NpGraph():
    def __init__(self,T=10,C=(1.0,1.0),phi=0.001,thresh=0.001,niter = 500):
        self._T = T#最大クラスタ数
        self._c1 = C[0]
        self._c2 = C[1]
        self._phi0=phi
        self._thresh = thresh
        self._niter = niter
        pass

    def fit(self,mat):
        lastll = -float('inf')
        self._mat = mat
        self._N = mat.shape[0]#頂点数
        self._tau1 = np.random.random(self._T)#クラス数分用意
        self._tau2 = np.random.random(self._T)#クラス数分用意
        self._psi = np.random.random((self._T,self._N))#T*N
        self._lambda1 = np.random.rand()
        self._lambda2 = np.random.rand()
        self._q = np.zeros((self._N,self._T))#頂点ごとの帰属確率
        self._s = np.zeros((self._N,self._T))#si
        self._phi = np.zeros((self._T,self._N))+self._phi0#phi
        for i in xrange(self._niter):
            self._estep(i)
            self._mstep()
            ll = self.loglikelihood()
            logging.info('Log Likelihood: %f',ll)
            if np.abs(ll-lastll)<self._thresh:
                break
            lastll = ll

    def loglikelihood(self):
        return np.sum(self._q*self._s)

    def _estep(self,iter):
        """
        E-step
        ν_r相当が無駄に計算されているので直す
        ここでいうs_irは、Newmanのモデルにおけるπ_r*Π_j(θrj*Aij), i=1...N,r=1..T
        digammaを通しているので期待値がlogスケールになっている
        """
        #calculate s_i_r
        for i in range(self._N):
            for r in range(self._T):
                self._s[i,r]=0.0
                #stick-breakingでυ(Newmanのπに相当)を計算
                if r<(self._T-1):#0オリジンの為1引く
                    self._s[i,r]+=digamma(self._tau1[r])-digamma(self._tau1[r]+self._tau2[r])#最後以外計算
                if r>0:#0オリジンのため
                    for l in range(r):#最初以外自分以下のインデックスの値を計算
                        self._s[i,r]+=digamma(self._tau2[l])-digamma(self._tau1[l]+self._tau2[l])
                #対数スケールでのΠθrj*Aijに相当
                for j in range(self._N):
                    self._s[i,r]+=self._mat[j,i]*(digamma(self._psi[r,j])-digamma(np.sum(self._psi[r])))

        #calculate q_i_r
        for i in range(self._N):
            for r in range(self._T):
                self._q[i,r]=np.exp(self._s[i,r])/np.sum(np.exp(self._s[i]))

    def _mstep(self):
        #update tau1
        self._tau1 = self._q.sum(axis=0)+1.0
        #update tau2
        for r in range(self._T):
            self._tau2[r] = self._lambda1/self._lambda2
            for i in range(self._N):
                self._tau2[r]+=np.sum(self._q[i,r+1:])

        #update psi
        #ここ高速化出来るはず
        for r in range(self._T):
            for j in range(self._N):
                self._psi[r,j] = self._phi[r,j]
                for i in range(self._N):
                    self._psi[r,j]+=self._mat[j,i]*self._q[i,r]

        #update lambda
        self._lambda1 = self._c1+self._T-1
        self._lambda2 = self._c2
        for k in range(self._T-1):#このループ消せるはず
            self._lambda2 -= (digamma(self._tau2[k])-digamma(self._tau1[k]+self._tau2[k]))

    def get_cluster(self):
        return np.argsort(self._q,axis=1).T[-1]

if __name__ == '__main__':
    logging.basicConfig(level=logging.INFO,format='%(asctime)s %(levelname)s %(message)s')
    #隣接行列
    mat = np.array([[0,1,1,1,0,0,0,0,0],
                    [1,0,1,1,0,0,0,0,0],
                    [1,1,0,1,0,0,0,0,0],
                    [1,1,1,0,0,0,1,0,0],
                    [0,0,1,0,0,1,1,1,1],
                    [0,0,0,0,1,0,1,1,1],
                    [0,0,0,0,1,1,0,1,1],
                    [0,0,0,0,1,1,1,0,1],
                    [0,0,0,0,1,1,1,1,0]])
    mat.astype(np.float64)
    g = NpGraph(T=5,C=(1.0,1.0),phi=1.0)
    g.fit(mat)
    print g._q
    print g.get_cluster()