GraphLassoによる変数間の関係のグラフ化

import pydot
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib import pyplot as plt
from sklearn.datasets import load_boston
from sklearn.covariance import GraphLasso
import scipy as sp

#X:データ, alpha:L1正則化パラメータ, disp_heatmap:結果を表示するかどうか
def Graphical_Lasso(X, alpha, disp_heatmap=True):
    #データの正規化（必須）
    X=sp.stats.zscore(X,axis=0)
    
    #GraphLasso
    model = GraphLasso(alpha=alpha,verbose=True)
    model.fit(X)
    
    cov=np.cov(X.T) #計算による分散共分散行列(転置を取るかはデータの向きによる)
    cov_ = model.covariance_ #スパース化した分散共分散行列
    pre_ = model.precision_ #スパース化した分散共分散行列の逆行列
    
    #分散共分散行列のヒートマップによる表示
    if disp_heatmap==True:
        f, axes = plt.subplots(1, 3, figsize=(15, 4))
        axes[0].set_title('Sample covariance')
        axes[1].set_title('Estimated covariance')
        axes[2].set_title('Estimated precision')
        
        sns.heatmap(cov, ax=axes[0])
        sns.heatmap(cov_, ax=axes[1])
        sns.heatmap(pre_, ax=axes[2])
    return model

def Disp_Gaussian_Graphical_Model(model,feature_names):
    pre_ = model.precision_ #スパース化した分散共分散行列の逆行列
    
    #グラフ表示のために対角成分は0にする
    pre_zero=pre_-np.diag(np.diag(pre_))
    
    #分散共分散行列の逆行列(≒隣接行列)からグラフを生成する
    g = pydot.Dot(graph_type='graph')
    df = pd.DataFrame(index=feature_names.tolist(),
                      columns=feature_names.tolist(), 
                      data=pre_zero.tolist())
    
    #ノードを追加
    for c in df.columns:
        node = pydot.Node(c)
        g.add_node(node)
    #エッジを追加    
    for i in df.index:
        for c in df.columns:
            if c>=i:
                #エッジの絶対値が0.1より大きいときに表示(ここは好みによる)
                if abs(df.loc[i, c]) > 0.1:
                    edge = pydot.Edge(g.get_node(i)[0], g.get_node(c)[0],
                                      penwidth=5*abs(df.loc[i, c]))
                    edge.set_label('{:.1f}'.format(df.loc[i, c]))
                    g.add_edge(edge)
    
    #結合数によって色を変える
    nodeList = g.get_node_list()
    num_edges=np.sum(abs(pre_zero) > 0.1, axis=0)
    max_num_edges=max(np.amax(num_edges),1)
    for i in range(len(df.index)):        
        node = nodeList[i]
        # H S Vの順 Sを変更
        color=str(0)+" "+str(num_edges[i]/max_num_edges)+" "+str(1)
        node.set_color("black") #nodeの輪郭の色
        node.set_fillcolor(color) #nodeの中身の色
        node.set_style('filled')
    
    # グラフを出力
    #g.write_png('LassoEstimatedGGM.png', prog='neato')
    g.write_png('LassoEstimatedGGM.png')

def Compute_extended_BIC(model,n,p,cov,gamma):
    omega = model.precision_
    E=(np.sum(np.sum(omega != 0, axis=0))-p)*0.5
    #尤度
    MLE=n*(np.log(np.linalg.det(omega))-np.trace(np.dot(cov,omega)))
    
    EBIC=-MLE+E*np.log(n)+4*E*gamma*np.log(p)
    return EBIC

def main():
    #Boston Housing Dataをインポート
    boston = load_boston()
    X=boston.data
    feature_names=boston.feature_names
    
    cov=np.cov(X.T) #標本分散共分散行列
    #普通はn(サンプル数)>p(変数の数)なので
    n=max(X.shape[0],X.shape[1]) #サンプル数
    p=min(X.shape[0],X.shape[1]) #変数の数(ノード数)
    
    #Graphical Lassoの実行
    alpha=0.4 #L1正則化パラメータ
    model=Graphical_Lasso(X, alpha)
    
    #BIC, EBICの計算
    BIC=Compute_extended_BIC(model,n,p,cov,0)
    print("BIC  :"+str(BIC)+" alpha :"+str(alpha))
    EBIC=Compute_extended_BIC(model,n,p,cov,0.5)
    print("EBIC :"+str(EBIC)+" alpha :"+str(alpha))
    
    #グラフ化
    Disp_Gaussian_Graphical_Model(model,feature_names)
    
if __name__ == '__main__':
    main()

def Compute_BEST_GGM_from_extended_BIC(X,n,p,cov,gamma):
    alphas = np.linspace(0, 1, 100)
    #alphas=np.array([0.1,0.2,0.3,0.4])
    num_alphas=alphas.shape[0]
    EBICs=np.zeros(num_alphas) #EBICを記録する配列
    
    for i in range(num_alphas):
        alpha=alphas[i]
        model=Graphical_Lasso(X, alpha,False)    
        #BIC, EBICの計算
        EBIC=Compute_extended_BIC(model,n,p,cov,gamma)
        print("EBIC :"+str(EBIC)+" alpha :"+str(alpha))
        EBICs[i]=EBIC
        
    #正則化パラメータに対するEBICの値のグラフ
    plt.plot(alphas, EBICs)
    plt.xlabel("L1 regularization parameter (alpha)")
    plt.ylabel("extended BIC (EBIC)")
    plt.show()
    
    #最適なmodelのインデックスを取得
    best_idx=np.argmin(EBICs)
    best_alpha=alphas[best_idx]
    best_model=Graphical_Lasso(X, best_alpha)  
    return best_model,best_alpha

def main_select_best_model():
    #Boston Housing Dataをインポート
    boston = load_boston()
    X=boston.data
    feature_names=boston.feature_names
    
    cov=np.cov(X.T) #標本分散共分散行列
    #普通はn(サンプル数)>p(変数の数)なので
    n=max(X.shape[0],X.shape[1]) #サンプル数
    p=min(X.shape[0],X.shape[1]) #変数の数(ノード数)
    model,alpha=Compute_BEST_GGM_from_extended_BIC(X,n,p,cov,0.5)
    print("Best alpha :"+str(alpha))
    #グラフ化
    Disp_Gaussian_Graphical_Model(model,feature_names)
    
if __name__ == '__main__':
    #main()
    main_select_best_model()