グラフィカルラッソ(Graphical Lasso)変数関係の可視化

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():

    #Dataをインポート

    dataset = sklearn.datasets.load_boston()

    X=dataset.data

    feature_names=dataset.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.1 #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()