• TOP
  • ブログ
  • 【前編】初心者のための時系列間の因果関係検出~らくらくノンパラメトリック・非定常・非線形・多変量に対...

ホットリンク 公式ブログ

データから見えてくる注目のトピックスをとりあげます。

テック

【前編】初心者のための時系列間の因果関係検出~らくらくノンパラメトリック・非定常・非線形・多変量に対応~

開発本部研究開発グループR&Dチームのサンティです。

今回はエンジニア向けに時系列間の因果関係検出に関する論文「A Nonlinear Causality Estimator Based on Non-Parametric Multiplicative Regression」(以降[参考1])について紹介したいと思います。数学は四則計算レベルで、高度な知識は不要です。サンプルコードは、githubにアップロードしてあります。

 

概要

本来、時系列間の因果関係検出は時系列解析のエキスパートの聖域でした。初心者が安易に同じ因果関係検出ツールを使っても、見せかけや偽の因果関係を得てしまうケースが多いのではないでしょうか。なぜなら、因果関係を検証する前に予測モデル(多くの場合は様々な自己回帰)を構築する必要があるからです。この作業はエキスパートと初心者とではモデルの質がかなり異なります。因果関係の結論はモデルの質に依存するため、従来の因果関係検出タスクは初心者には優しくありません。かなり鬼門であると筆者は思います。

私は最近、[参考1]のような珍しい方法に出会い、工学的な問題解決の思考回路は大変素晴らしいと感じました。因果関係パターンを検出するために必要作業を十分行っており、省エネを良しとするエンジニアリング精神に満ちています。ただし、[参考1]の引用、支持、批判するような意見はウェブ上でまだ見受けられません。[参考1]の著者はプロフィールによると、医療用分野のコンテストで受賞した実績もお持ちの方です。今回、ご紹介する技術は学術的な価値の根拠がありません。しかし、このような原石的な技術を知っておいても損はないと思います。

また、この因果関係検出方法の最も印象に残るポイントは、初心者にも優しい技術であることです。

そこで今回、前編後編の2回にわたり「初心者のための時系列間の因果関係検出」についてお伝えしたいと思います。

 

貢献ポイント

[参考1]の論文は因果関係検出に関する下記の課題を解決しようとして報告されたものです。(多少筆者の意見が入りますのでご了承ください。)

  1. パラメトリック手法: データの性質をパラメトリック的に説明できると仮定し、解析的にかつエレガントに解くような手法では、利用可能な条件に一致しない場合は無力です。ここはノンパラメトリック的に、データはデータなりに、本質を語ってもらいたいものです。
  2. 定常性・線形性: グレンジャー因果関係とその考え方を継承した研究、特に自己回帰モデルベースでは、大抵、時系列の定常性や線形性が必要です。また、非線形や非定常なデータに対して、なんらかの変換やカーネルの導入を要する手法も多々あります。何れも初心に優しくありません。
  3. 多変量対応: 元々のグレンジャー因果関係は2つの時系列データの因果関係を検証するもので3つ以上の時系列データの時系列を考慮できるように一苦労して対応する必要があります。
  4. 予測モデル構築: 専門的な知識が要してあるため初心者には難しい又は無理です。また、モデルの過学習問題を防ぐ対策が必要としています。
  5. 予測モデルなし: 移動エントロピー(Transfer Entropy: TE)や相互情報量系は確か、明確なモデルが不要だったと記憶していますが(誤っていたらすみません)、因果関係としての用途は[参考2]で様々な問題が指摘されています。この[参考2]では、移動エントロピーに公開処刑するかのように強くダメ出しをしていますので、ぜひ一読することをおすすめします。少なくとも、基本的なTEは定常性を必要としています。非定常データに対応可能なTE系手法もありますが、かなり苦労します。

[参考1]はグレンジャー因果関係の概念を引き継ぎながらこれらの問題を克服する方法を提案しています。アプローチとしては(正確な)モデル構築とモデルなしのちょうど真ん中くらいです。はたしてそのような効果があるのでしょうか、ご一緒に検証していきましょう。

※注意
サンプルコードはJulia言語バージョン0.5.2でこの原稿を作成しました。また、実験のため、ソースコードの最適化を考慮しておりませんのでご容赦下さい。前述の通り、サンプルコードはiPython形式でgithubにアップロードしてあります。

※おことわり
ここでご紹介する技術は弊社で提供している商品に組み込まれたり、参考にしたりはしていません。 すべては筆者自身の意見であり、内容の正確さの保証はありません。

 

Granger 因果関係のおさらい

時系列Yから時系列XへのGranger因果性ある(Y→X)というのは時系列Xの現在までのデータだけでモデルを作成するより時系列Yの現在までのデータも取り込んだ方が、有意により正確なXの予測モデルが構築できるということです。Granger因果性の定量測定は下記の式の通りです。σ^2はモデル残差の分散を表しています。伝統的なGranger因果関係検定(Y→X)は、超簡略化していうと、時系列Y付きモデルの残差の分散が時系列Yなしモデルより有意に小さいかどうかの検定です。Granger因果性には原因と結果の方向があります。逆方向の測定は式のXとYを入れ替えるだけで済むのです。


In [1]:

 [参考1]は「伝統的なGranger因果性の概念を継承したいけれど、データの性質をなんらかの型にはめたくないし、多変数間の関係をも明示的に仮定したくない」という著者の意図があります。そこで鍵となるのが「Non-Parametric Multiplicative Regression: NPMR」という生態学の専門家が提案したモデルです。詳細は[参考3]をご覧ください。[参考1]が伝統のGranger因果性検定の自己回帰モデルをこのNPMRモデルに置き換えることで、一気に、パラメトリック的な手法、定常性・線形性、多変量、正確なモデル構築、様々な課題から解放されます。

さて、NPMRモデルはどんなものか、軽く勉強しましょう。 NPMRモデルは目的の時系列Yを下記のような予測因子である多変量時系列Xで予測したい場合、予測値は式(2)で求めます。 パッと見でちょっと掛け算と足し算で重み付き平均を計算するような感じですね。 重みは予測因子データ同士の距離で決めているようです。


In [2]:
 
変数Xの2つの点の重みは下記の式のように計算します。
 

In [3]:
 
関数のコードはこのように書けます。
 
In [4]:
function calDist(x)
    exp(-0.5*x^2)
end
 
Out[4]:
calDist (generic function with 1 method)

In [5]:
  

式(2)を実装するとこんな感じになります。

 

In [6]:
function pred(X,Y,σ;cross_validate=true)
    T = length(Y)
    upper = zeros(T)
    lower = zeros(T)
    for t=1:T, i=1:T
        if cross_validate && i == t continue end
        prodw = prod(calDist.((X[i,:] - X[t,:]) ./ σ))
        lower[t] += prodw
        upper[t] += prodw*Y[i]
    end
    upper ./ lower
end
 
Out[6]:
pred (generic function with 1 method)

[参考1]では、XiからYへの因果性は式(5)で評価します。

In [7]:
 
 次に条件付き因果性に対応します。例えば第三の時系列Zを式(4)のように遅延エンベディングの集合が得られるとします。実装しますと下記となります。
 
In [8]:
function Cnpmr(Xy,Xxy,Yy;σy=NaN,σxy=NaN)
    Ŷy = pred(Xy,Yy,isnan(σy) ? ones(size(Xy)[2]) : σy)
    Ŷxy = pred(Xxy,Yy,isnan(σxy) ? ones(size(Xxy)[2]) : σxy)
    log(var(Ŷy-Yy)/var(Ŷxy-Yy))
end
Out[8]:
Cnpmr (generic function with 1 method)

In [9]:
 
 
この時、条件付き因果性は式(6)で評価します。

In [10]:
 
 条件付きの対応は、入力データの形成を整理するだけで、特に評価関数を書き換えることなく上記の実装のまま流用できるのは[参考1]の強くアピールしたいポイントでもあると論文を読んでいて感じました。因みにデータを整理するコードは下記のとおりです。なぜか計算部分より長いです。
 
In [11]:
function makeDataset(fromX::Array{Float64,2},
                toY::Array{Float64,2},d;Z=Array{Float64}(0,2))
    T = length(toY)-d
    Xy = vcat([toY[i:-1:(i-d+1)]' for i=d:(d+T-1)]...)
    if size(Z)[1] > 0
        Xz = vcat([collect((Z[i:-1:(i-d+1),:]'...))' for i=d:(d+T-1)]...)
        Xy = [Xy Xz]
    end
    Xx = vcat([collect((fromX[i:-1:(i-d+1),:]'...))' for i=d:(d+T-1)]...)
    Xxy = [Xx Xy]
    Yy = toY[(d+1):(d+T)]
    (Xy,Xx,Xxy,Yy[:,:])
end
Out[11]:
makeDataset (generic function with 1 method)
 

次に因果性の強さを求めます。この計算はNPMRモデルの本家のまま引用しています。エッセンスはデータを多少プラスとマイナス両方に変動させて結果への強さを観測するところにあります。あまり深く考えずにこの因果性の強さは絶対値ではなく、相対的なものとして捉えておけばよいと思います。

 

In [12]:
実装しますとnudge、Q、sensitivity関数のようになります。少しめんどくさいですが、論文の添付資料の結果と一致していますので(多分)正しいです。
 
In [13]:
function nudge(X::Array{Float64,2},col,)
    dk = *abs(maximum(X[:,col]) - minimum(X[:,col]))
    nX = repeat(X; inner=[2,1])
    nX[:,col] += repeat([dk; -dk], outer=[size(X)[1],1])
    nX
end
Out[13]:
nudge (generic function with 1 method)

In [14]:
function Q(X,Y,col;σ=NaN,=0.05)
     = pred(X,Y,isnan(σ) ? ones(size(X)[2]) : σ,cross_validate=false)
    nX = nudge(X,col,)
    nY = repeat(Y,inner=[2,1])
    Ŷs = pred(nX,nY,isnan(σ) ? ones(size(X)[2]) : σ,cross_validate=false)
    Qcol = sum(abs(Ŷs - repeat(,inner=[2,1]))) /
            (2*length(Y)**abs(maximum(Y)-minimum(Y)))
end
Out[14]:
Q (generic function with 1 method)

In [15]:
function sensitivity(toY::Array{Float64,2},fromX::Array{Float64,2};d=2)
    (Xy,Xx,Xxy,Yy) = makeDataset(fromX,toY,d)
    xw = size(fromX)[2]
    for col in 1:(d*xw)
        Qcol = Q(Xxy,Yy,col)
        println("Q_x$(col%xw==0?xw:col%xw)(t-$(Int(ceil(col/xw)))): $Qcol")
    end
end
Out[15]:
sensitivity (generic function with 1 method)

 次に統計的検出した因果性が有意かどうかをはっきりさせます。ノンパラメトリックを選んだゆえに解析的な解がありません。 偶然レベルではないと実践で検証しかありません。この時、[参考1]は時差的なSurrogateデータで検証するのが良いと他論文の結論を引用しています。 やり方はとても簡単です。時系列をランダムに長さの3分の1以上シフトすればよいのです。ここでは最大90%シフトと上限を設けます。 下記はデータを適当にぐるぐる回転させるコードです。
 
In [16]:
function surrogate(X;lower=33,upper=90)
    circshift(X,Int(floor(size(X)[1]*rand(lower:upper)/100)))
end
Out[16]:
surrogate (generic function with 1 method)

 次は肝心の有意検証です。まず、算出した因果性の値がゼロ以下の場合、即刻、因果性がないと判断します。そうでない場合はまず、欲しい有意水準αを決め、検証回数を計算します。何回検証すればよいかの計算式がありますので[参考1]の2.4節をご覧ください。 有意の判断はとても簡単です。決まった回数でぐるぐる回転させたデータで因果性の強さを取得します。そのなかに元データの因果性より一つでも因果性が大きいものがあれば、因果性が有意ではないと判断します。なかったら、因果性を認め、trueを返します。実装しますと下記のようなコードとなります。
 
In [17]:
function isCausal(
            fromX::Array{Float64,2},
            toY::Array{Float64,2};d=2,α=0.05,Z=Array{Float64}(0,2))
    (Xy,Xx,Xxy,Yy) = makeDataset(fromX,toY,d,Z=Z)
    Cmain = Cnpmr(Xy,Xxy,Yy)
    println("Cnpmr(X->Y): $Cmain")
    if Cmain < 0 return false end
    surnum = Int(floor(1/α - 1))
    for cnt in 1:surnum
        sfromX = surrogate(fromX)
        (sXy,sXx,sXxy,Yy) = makeDataset(sfromX,toY,d)
        Csurr = Cnpmr(sXy,sXxy,Yy)
        println("$cnt: $Csurr")
        if Csurr > Cmain return false end
    end
    true
end
Out[17]:
isCausal (generic function with 1 method)

 下記は実践編の検証のためのコードです。やっていることは因果性検証を数回(10回)を繰り返し検証するだけです。
 
In [18]:
function batchTest(datfn;d=3,trials = 10)
    x1_x2sigcnt = 0
    x2_x1sigcnt = 0
    (X1,X2,Z) = datfn()
    for i = 1:trials
        if isCausal(X1,X2,d=d) x1_x2sigcnt += 1 end
        if isCausal(X2,X1,d=d) x2_x1sigcnt += 1 end
    end
    println("X1 -> X2: $x1_x2sigcnt / $trials significant count(s)")
    println("X2 -> X1: $x2_x1sigcnt / $trials significant count(s)")
end
Out[18]:
batchTest (generic function with 1 method)

 

前編はここまでとなります。後編は、10月12日に公開予定です。どうぞお楽しみに!

参考

  1. Nicolaou Nicoletta, Constandinou Timothy G.: "A Nonlinear Causality Estimator Based on Non-Parametric Multiplicative Regression", Frontiers in Neuroinformatics vol.10, 2016 https://doi.org/10.3389/fninf.2016.00019
  2. Ryan G. James, Nix Barnett, James P. Crutchfield: "Information Flows? A Critique of Transfer Entropies", Phys. Rev. Lett. 116, 238701 (2016)
  3. McCune, Bruce: "Non-parametric habitat models with automatic interactions", Journal of Vegetation Science 17: 819-830, 2006
 
▼ホットリンク では一緒に働く仲間を募集しています!