<jupyter_start><jupyter_text>仕訳データのグラフ表現と勘定科目の相関関係 簿記でいう勘定連絡図のようなものを有効グラフで表現してみます。勘定連絡図についてはググってください。  0. セットアップ 1. matplotlibに日本語のフォントを追加、設定します。2. 今回は有限会社アイソフトの製品「税理士いらず」のサンプルデータ「株式会社サンプルデザイン」をお借りします。以下のサイトから仕訳データをダウンロードし、好きな場所に置きます（以下ではこのスクリプトファイルと同じディレクトリに"testdata.csv"の名称で置いてあることを想定してます）。http://www.z-irazu.jp/sample.html  1. データの読み込み<jupyter_code>import numpy as np
import pandas as pd
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
import networkx as nx
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster
import japanize_matplotlib
filename='../Developing/testdata.csv'
JE_raw=pd.read_csv(filename, \
                   header=None,index_col=False, engine="python", encoding="SHIFT-JIS",\
                   names=['date','debit_account_name','debit_amount','debit_tax','credit_account_name','credit_amount','credit_tax','description_1','description_2'])<jupyter_output><empty_output><jupyter_text>2. 行列データの作成<jupyter_code>Names=np.concatenate([JE_raw.loc[:,'debit_account_name'].values,JE_raw.loc[:,'credit_account_name'].values])
TBindex=np.unique(Names)
transaction_mat=pd.DataFrame(np.zeros((TBindex.shape[0],TBindex.shape[0])),dtype=int)
transaction_mat.index=TBindex
transaction_mat.columns=TBindex
for itr in range(JE_raw.shape[0]):
    rcd=JE_raw.iloc[itr,:]
    mask_d=rcd["debit_account_name"] == TBindex
    mask_c=rcd["credit_account_name"] == TBindex
    transaction_mat.iloc[mask_d,mask_c]=transaction_mat.iloc[mask_d,mask_c]+rcd['debit_amount']
    
transaction_mat_net=pd.DataFrame(np.zeros((TBindex.shape[0],TBindex.shape[0])),dtype=int)
transaction_mat_net.index=TBindex
transaction_mat_net.columns=TBindex
transaction_mat_net.loc[:,:]=(transaction_mat.values+transaction_mat.values.transpose())/2<jupyter_output><empty_output><jupyter_text>3. 無向グラフ<jupyter_code>df=transaction_mat_net
threshold=0.000001
G=nx.Graph()
edge_labels={}
for a, b in itertools.combinations(df.columns, 2):
        val1 = df.loc[a, b]
        if abs(val1) > threshold:
            edge_labels[(a,b)]=np.round(val1/1000000,1).astype(str) + " M"
            weigth_tmp=10 * (abs(np.log10(val1)/7) - threshold) / (1 - threshold)
            G.add_edge(a,b, weight=weigth_tmp)

pos = nx.spring_layout(G) # グラフのポジションを力学モデルで決定します。

#plt.rcParams['font.family'] = 'Osaka'
plt.figure(figsize=(50,20))

node_size = 10 
nx.draw_networkx_nodes(G, pos, node_color="b",alpha=0.6, node_size=node_size)
nx.draw_networkx_labels(G, pos, fontsize=32, font_family='Osaka', font_weight="bold")
edge_width = [ d["weight"] for (u,v,d) in G.edges(data=True)]
nx.draw_networkx_edges(G, pos, alpha=0.5, edge_color="g", width=edge_width)
nx.draw_networkx_edge_labels(G, pos,fontsize=32,edge_labels=edge_labels,font_weight="bold")
plt.figure(figsize=(50,20))<jupyter_output>/anaconda3/lib/python3.6/site-packages/networkx/drawing/nx_pylab.py:611: MatplotlibDeprecationWarning: isinstance(..., numbers.Number)
  if cb.is_numlike(alpha):<jupyter_text>4. 勘定の相関関係 ここから本題です。仕訳の借方と貸方の関係になる勘定科目同士（上の図でエッジでつながれた科目同士）は同じタイミングで増減するため、相関関係がありそうです。そこで、勘定科目間の相関行列（相関グラフ）を作成してみます。<jupyter_code>def ConvertToJEvector(JE_raw):
    Names=np.concatenate([JE_raw.loc[:,'debit_account_name'].values,JE_raw.loc[:,'credit_account_name'].values])
    TBindex=np.unique(Names)
    JE_vec=pd.DataFrame(np.zeros((TBindex.shape[0],JE_raw.shape[0])),dtype=int)
    JE_vec.index=TBindex
    for itr in range(JE_raw.shape[0]):
        rcd=JE_raw.iloc[itr,:]
        mask_d=rcd["debit_account_name"] == TBindex
        mask_c=rcd["credit_account_name"] == TBindex
        JE_vec.iloc[mask_d,itr]=rcd['debit_amount']
        JE_vec.iloc[mask_c,itr]=rcd['credit_amount']*(-1)
    return JE_vec
        
JE_vec=ConvertToJEvector(JE_raw)<jupyter_output><empty_output><jupyter_text>相関行列ベクトルデータは借方の増加方向をプラス、貸方の増加方向をマイナスにしています。したがって、例えば、以下の仕訳で現金と預金は負の相関関係になります。現金 100 / 預金 100そこで、相関行列は-1倍して借方と貸方の関係になる場合に大きい値をとるようにしています。<jupyter_code>JE_Cor=np.corrcoef(JE_vec)*(-1)
plt.figure(figsize=(12,10))
ax=sns.heatmap(JE_Cor,cmap='jet')<jupyter_output><empty_output><jupyter_text>相関グラフ<jupyter_code>df=pd.DataFrame((-1)*JE_Cor)
df.index=TBindex
df.columns=TBindex
threshold=0.05
G=nx.Graph()
edge_labels={}
for a, b in itertools.combinations(df.columns, 2):
        val1 = df.loc[a, b]*(-1)
        if abs(val1) > threshold:
            edge_labels[(a,b)]=np.round(val1,1).astype(str)
            #weigth_tmp=10 * (abs(np.log10(val1)/7) - threshold) / (1 - threshold)
            weigth_tmp=10 * (val1 - threshold) / (1 - threshold)
            G.add_edge(a,b, weight=weigth_tmp)        

plt.rcParams['font.family'] = 'Osaka'
plt.figure(figsize=(50,20))

node_size = 10 
nx.draw_networkx_nodes(G, pos, node_color="b",alpha=0.6, node_size=node_size)
nx.draw_networkx_labels(G, pos, fontsize=32, font_family='Osaka', font_weight="bold")
edge_width = [ d["weight"] for (u,v,d) in G.edges(data=True)]
nx.draw_networkx_edges(G, pos, alpha=1, edge_color="b", width=edge_width)
nx.draw_networkx_edge_labels(G, pos,fontsize=32,edge_labels=edge_labels,font_weight="bold")
plt.figure(figsize=(50,20))<jupyter_output>/anaconda3/lib/python3.6/site-packages/networkx/drawing/nx_pylab.py:611: MatplotlibDeprecationWarning: isinstance(..., numbers.Number)
  if cb.is_numlike(alpha):<jupyter_text>相関係数が-0.1より小さい科目同士にエッジを描いてます。例えば、売上と売掛金は相関係数が-0.7になっており、借方貸方の関係になりやすいことが分かります。逆に仕訳グラフでは預金や現金はいろんな科目と繋がってましたが、相関グラフでは少なくなっています。預金や現金はいろいろな科目によって変動するため、ある特定の科目の増減だけではあまり説明がつきません。このため、預金や現金周りの相関係数は0に近くなっています。  5. 仕訳グラフに基づく重回帰モデル 上のグラフの相関関係を用いれば、ある勘定から他の勘定科目を予想できます。一方で、現金や預金を予測するにはいろいろな科目を説明変数に用いる必要があります。そこで重回帰モデルを考えます。説明変数の選択は難しい問題ですが、今回は仕訳グラフで結合しているノードの科目を説明変数にします。ちなみに、EYの特許では変数選択にLASSOを用いてますが、仕訳が既知であれば真のモデル（仕訳グラフ）を知っているため、それを用いてパラメータの推定を用います。（前期の仕訳データを入手しており、当期の勘定残高の増減を予想する場合などに応用可能です。逆に仕訳データに依拠できない場合はこの方法は有効ではありません。）<jupyter_code>mask=transaction_mat_net>0
from sklearn.preprocessing import StandardScaler
ss=StandardScaler()
JE_vec_z=pd.DataFrame(ss.fit_transform(JE_vec.transpose()).transpose(),index=TBindex)
from sklearn import linear_model
from sklearn.metrics import mean_squared_error
MSE=pd.Series(np.zeros((len(TBindex))),index=TBindex)
BetaMat=pd.DataFrame(np.zeros((len(TBindex),len(TBindex))),index=TBindex,columns=TBindex)


import warnings
warnings.filterwarnings('ignore')
for itr in range(len(TBindex)):
    est = linear_model.LinearRegression()
    target=JE_vec_z.loc[TBindex[itr],:].as_matrix()
    exp_data=JE_vec_z.loc[mask.loc[TBindex[itr],:],:]

    exp=exp_data.transpose().as_matrix()
    est.fit(exp,target)
    pred = est.predict(exp)
    BetaMat.loc[TBindex[itr],mask.loc[TBindex[itr],:]]=est.coef_
    MSE.loc[TBindex[itr]]=mean_squared_error(target,pred)<jupyter_output><empty_output><jupyter_text>6. グラフ上の一つ向こうのノードまで使用した重回帰モデル Aを隣接行列とすると、$S_{n}=\sum_{k=1}^{n}A^{k} $ の (i,j)成分$S_{(i,j)}$ はノードiからノードjまでの経路の場合の数（同じ経路を複数回通ることを許す）となる。したがってノードiのグラフ上で一つ向こうまでのノードの範囲はBoolian行列B(1または0で論理和と論理積が定義)を用いて$B^{2} + B - E$ の(i,:)で与えられる。<jupyter_code>from sklearn import linear_model
from sklearn.metrics import mean_squared_error
MSE_2=pd.Series(np.zeros((len(TBindex))),index=TBindex)
BetaMat_2=pd.DataFrame(np.zeros((len(TBindex),len(TBindex))),index=TBindex,columns=TBindex)

mask_2=pd.DataFrame(np.dot(mask,mask)*1+mask*1-np.eye(len(TBindex)),index=TBindex, columns=TBindex)

import warnings
warnings.filterwarnings('ignore')
for itr in range(len(TBindex)):
    est = linear_model.LinearRegression()
    
    target=JE_vec_z.loc[TBindex[itr],:].as_matrix()
    exp_data=JE_vec_z.loc[mask_2.loc[TBindex[itr],:]>0,:]

    exp=exp_data.transpose().as_matrix()
    est.fit(exp,target)
    pred = est.predict(exp)
    BetaMat_2.loc[TBindex[itr],mask_2.loc[TBindex[itr],:]>0]=est.coef_
    MSE_2.loc[TBindex[itr]]=mean_squared_error(target,pred)<jupyter_output><empty_output><jupyter_text>説明変数の数が多くなるのでLASSOを用います。<jupyter_code>MSE_LASSO=pd.Series(np.zeros((len(TBindex))),index=TBindex)
BetaMat_LASSO=pd.DataFrame(np.zeros((len(TBindex),len(TBindex))),index=TBindex,columns=TBindex)
for itr in range(len(TBindex)):
    est = linear_model.LassoCV(alphas=10 ** np.arange(-10, 1, 0.1), cv=5, max_iter=10000)
    
    target=JE_vec_z.loc[TBindex[itr],:].as_matrix()
    exp_data=JE_vec_z.loc[mask_2.loc[TBindex[itr],:]>0,:]

    exp=exp_data.transpose().as_matrix()
    est.fit(exp,target)
    pred = est.predict(exp)
    BetaMat_LASSO.loc[TBindex[itr],mask_2.loc[TBindex[itr],:]>0]=est.coef_
    MSE_LASSO.loc[TBindex[itr]]=mean_squared_error(target,pred)<jupyter_output><empty_output><jupyter_text>分散説明率%（0〜100で0は全く予測できてない）を比較してみます。モデル選択ではAICなどのように説明変数の個数で罰則を設けるのが通常ですが、今回はあてはまれば良しとします。（また、本当はtrainingとtestデータを分離する必要があります。）説明変数を一つ向こうのノードまで増やすことにより、売上高、売掛金、普通預金、現金、給料手当で大きく説明率が上昇してます。その他は説明変数の増加によるわずかな増加です。ほとんど100%になると思っていたのですが（勘定科目の増減は仕訳の相手の勘定科目の増減で説明できるはず）、そううまく行きませんでした。ちなみにLASSOは縮小推定で分散説明率が落ちているだけでした。<jupyter_code>100-round(pd.concat([MSE,MSE_2,MSE_LASSO],axis=1)*100,2)<jupyter_output><empty_output><jupyter_text>ベータ（回帰係数）行列は以下のようになります。<jupyter_code>np.dot(mask,mask)
plt.figure(figsize=(6,5))
ax=sns.heatmap(BetaMat_2,cmap='jet')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14933774834437086, 'fraction_numerical': 0.02041942604856512, 'mean_word_length': 8.227087576374746, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '544751', 'license': 'MIT', 'n_tokens_mistral': 4970, 'n_tokens_neox': 4610, 'n_words': 311, 'path': 'JE_graph_Linear_Modeling.ipynb', 'repo_name': 'nororo/AuditDataAnalysis', 'size': '9060'}