離散的でない実データの相互情報量とtransfer entropyの推定
2つの確率変数が、互いに連動しているかを知る方法として、単純には、相関係数を調べればよいけど、相互情報量を見るという方法もある。一般には、確率変数が独立であることと相互情報量が0であることは同値であるけど、相関係数が0であることは確率変数が独立であるための必要条件でしかない。つまり、
相互情報量=0 <=> 確率変数が独立 => 相関係数が0
まぁ、相関係数が0でも"高次の相互モーメント"(という呼び方があるのか知らないけど)が残ってる可能性はあるけど、相互情報量が0であれば、高次のモーメントが全部消えてることが保証される
相関がないけど、相互情報量が0にならない例として、S1,S2を区間[-1,1]上の独立な一様分布で、X1=S1+S2 , X2=S1-S2とすると、X1,X2は相関はないけど、独立ではない(何の知識もなしに観測すると、X1の絶対値が大きい時は、X2が小さく、逆にX2の絶対値が大きい時は、X1の絶対値は小さいという規則性を発見するはず)。X1,X2から、S1,S2への分解は、独立成分分析で出来る(はず)
別の例としては、Xが適当な分布に従う時系列で、Yは前半はXに等しいけど、後半は-Xに等しい時系列みたいなケースでも、相関は0になるけど、何らかの形で関連しあっていることを疑うはず。人工的に作ったわけでない現実のデータで、相関係数を見ても相関はないけど、相互情報量を見ると、関連性が見える場合があるのかは知らない
相互情報量の定義は、あちこちに書いてあるけど、大体、分布関数を既知とした定義しか書いてなくて、同時確率分布を勝手に仮定したくない場合は多い(ガウス分布とか混合ガウス分布を仮定して推定してしまうような場合もあるかもしれないけど)だろうので、(データが離散的でない時には)何らかの計算方法を考える必要がある。平均や分散のような統計量は、自然な推定方法が確立されてるけど、相互情報量については、多くの場合、確率分布をヒストグラムで近似するという方法が採用されている。幅は均等に分割するとして、ヒストグラムの分割数をどう決めるかが問題だけど、これも色んな方法がある
ヒストグラム
http://ja.wikipedia.org/wiki/%E3%83%92%E3%82%B9%E3%83%88%E3%82%B0%E3%83%A9%E3%83%A0#.E3.83.93.E3.83.B3.E3.81.AE.E6.95.B0.E3.81.A8.E5.B9.85
一番有名なのは、スタージェスの公式というやつだと思うけど、最近は、赤池情報量基準とか使って決めると、それっぽいかもしれない(本当に実用的かどうかは知らないけど)。と、そんな感じで、pythonで実装して見たのが、以下
#!/usr/bin/env python
import sys
#一次元時系列同士の相互情報量
def mic(xseq , yseq):
xmax,xmin = max(xseq),min(xseq)
ymax,ymin = max(yseq),min(yseq)
N = min(len(xseq) , len(yseq))
L, Nx , Ny = sys.float_info.max , 1 , 1
if xmax==xmin or ymax==ymin:
return 0.0
#-- 赤池情報量基準に基づいて分割数を決める
for n in xrange(2,N):
p = dict()
for (x,y) in zip(xseq,yseq):
nx = int( n*(x-xmin)/(xmax-xmin) )
ny = int( n*(y-ymin)/(ymax-ymin) )
p[(nx,ny)] = p.get( (nx,ny) , 0) + 1
AIC = -(sum( [Nxy*log(float(Nxy)/N , 2) for Nxy in p.itervalues()] ) + N*log(n,2)*2) + (n*n-1)
if AIC < L:
L = AIC
Nx = n
Ny = n
#-- 相互情報量の推定
px,py,p = dict(),dict(),dict()
for (x,y) in zip(xseq,yseq):
nx = int( Nx*(x-xmin)/(xmax-xmin) )
ny = int( Ny*(y-ymin)/(ymax-ymin) )
px[nx] = px.get(nx,0)+1
py[ny] = py.get(ny,0)+1
p[(nx,ny)] = p.get( (nx,ny) , 0) + 1
return sum([p0*log(N*float(p0)/(px[nx]*py[ny]) , 2) for ((nx,ny) , p0) in p.iteritems()])/(N*N)
テストデータとして、100日分の日経平均の収益率と個別銘柄の収益率の相互情報量を計算して、大きい順にいくつか並べてみると以下のようになる。
| 証券コード | 銘柄名 | 相互情報量 | 相関係数 | 日経225構成率 |
|---|---|---|---|---|
| 6473 | (株)ジェイテクト | 0.0104 | 0.8327 | 0.33% |
| 8015 | 豊田通商(株) | 0.0101 | 0.8341 | 0.77% |
| 4901 | 富士フィルムホールディングス(株) | 0.0095 | 0.7232 | 0.64% |
| 3402 | 東レ(株) | 0.0094 | 0.7797 | 0.19% |
| 7267 | ホンダ | 0.0093 | 0.8589 | 2.17% |
| 4324 | (株)電通 | 0.0091 | 0.8098 | 1.03% |
| 5901 | 東洋製罐グループホールディングス(株) | 0.0091 | 0.8349 | 0.44% |
| 6502 | (株)東芝 | 0.0089 | 0.7477 | 0.14% |
| 3405 | (株)クラレ | 0.0089 | 0.8112 | 0.41% |
相互情報量は相関係数と違って、0.5くらいあればそこそこ相関がありそうというような絶対的な判定には使いづらいけど、一応日経225採用銘柄ばっかだし、合ってはいそうな気がする。binの分割数は、大体7前後になるけど、これはスタージェスの公式でもlog_2(100)+1≒7~8なので、あんまり変わらない。
ついでに、transfer entropyというものも、実装してみる。transfer entropyは、因果指標の一つと位置付けられていて、一方の変数が他方の変数に先行している、もしくは追従しているというのを見る道具。因果関係を見るのに、異時刻相関や遅延相互情報量を測るのと、アイデアとしては同じようなものだと思う。まぁ、以下のスライドとか見れば意味は分かるはず
Transfer Entropy
http://users.utu.fi/attenka/TEpresentation081128.pdf
transfer entropyを誰が考えたのか知らないけど、これは計量経済学で以前から知られていたGrangerの因果性と同じらしい(論文未読)
Granger causality and transfer entropy are equivalent for Gaussian variables
http://arxiv.org/abs/0910.4514
Equivalence of Granger causality and transfer entropy: a generalization
http://www.m-hikari.com/ams/ams-2011/ams-73-76-2011/schindlerAMS73-76-2011.pdf
やり方としては、相互情報量と同じで、ヒストグラムに分割して計算する。分割数は面倒なのでスタージェスの公式で決めた
#!/usr/bin/env python
#一次元時系列同士のtransfer entropy
#xseq <- yseqを測定
#時系列は、後ろにいくほど古いとする
def tec(xseq , yseq):
xmax,xmin = max(xseq),min(xseq)
ymax,ymin = max(yseq),min(yseq)
N = min(len(xseq) , len(yseq)) - 1
Nx = int( log(N , 2) + 1 )
Ny = Nx
p , px , pxx , pxy = dict() , dict() , dict() , dict()
for (x1 , (x,y)) in zip(xseq[:-1] , zip(xseq[1:] , yseq[1:])):
nx1 = int( Nx*(x1-xmin)/(xmax-xmin) )
nx = int( Nx*(x-xmin)/(xmax-xmin) )
ny = int( Ny*(y-ymin)/(ymax-ymin) )
p[(nx1,nx,ny)] = p.get( (nx1,nx,ny) , 0) + 1
px[nx] = px.get(nx , 0) + 1
pxx[(nx1,nx)] = pxx.get( (nx1,nx) , 0) + 1
pxy[(nx,ny)] = pxy.get( (nx,ny) , 0) + 1
R=[p0*log( float(p0*px[nx])/(pxx[(nx1,nx)]*pxy[(nx,ny)]) , 2) for ((nx1,nx,ny) , p0) in p.iteritems()]
return sum(R)/(N-1)
#-- test
import random
if __name__=="__main__":
xs = [random.uniform(-1,1) for _ in xrange(100)]
ys = [2*x+0.01*random.random() for x in xs[1:]]
xs = xs[:-1] #長さをysに合わせておく
print "T(X<-Y) = %2.5f" % tec(xs , ys)
print "T(X->Y) = %2.5f" % tec(ys , xs)
print "I(X,Y) = %2.5f" % mic(xs , ys)
print "I(X',Y) = %2.5f" % mic(xs[1:] , ys[:-1])
print "I(X'',Y) = %2.5f" % mic(xs[2:] , ys[:-2])で、上記のような人工的なテストデータを作って測定すると
T(X<-Y) = 1.43836 T(X->Y) = 2.53167 I(X,Y) = 0.00225 I(X',Y) = 0.03056 I(X'',Y) = 0.00252
みたいな結果が返ってくるので、xsの値がysの値に先行していると判定できる？(多分、本当はGranger testみたいに、ちゃんと検定しないと、何とも言えないっぽい)。まぁ、この例だと、一方の時系列をシフトした相互情報量を見た方が良い感じもするけど
