Convergent Cross Mapping

# お詫び
間に合わなかったので公開後にちょくちょく書き足すことにします…
2018/07/03　情報論的アプローチについていろいろと追加，まだ編集中の部分もあるがとりあえず公開したままにしておきます．とくに肝心のCCMの章がほとんど書きかけ．

# モチベーション
相関関係は因果関係を意味しない。よって因果関係を発見するには相関関係のチェックとは別の手法が必要となる.統計的因果推論の文脈では，因果関係は「要因Xを変化させたときに要因Yも変化する場合に因果関係があるという」(岩波データサイエンスvol.3 p.10より引用)などと定義されている．
因果関係とは何かということについては科学哲学的な分野でそこそこ議論されていた話のようだ．この方向の議論については[このスライド]
(https://www.slideshare.net/JunOtsuka/2015-52074249)が参考になる。
歴史的な経緯はさておき，因果関係を観測・データなどから推定する枠組みを因果推定と呼ぶ．

# 因果推定の手法あれこれ
## 統計的因果推定
因果推定といったときおそらくこれが本流の研究である。公衆衛生であったり、医学の分野においてある介入を行ったときそれが効果があるか？というのを統計的に判断する指針を与えるのが統計的因果推論という分野である.
[機械学習プロフェッショナルシリーズ](http://www.kspub.co.jp/book/series/S043.html)にも[統計的因果探索](http://www.kspub.co.jp/book/detail/1529250.html)なるものがでており，人工知能研究の一つの中心的なトピックといえるだろう．
本稿では（詳しくないので）詳しく触れないが，かなりベイズっぽい印象を持っている．ここら辺の分野では，人工知能の巨人の一人ともいわれるJudea PearlによるCausality-Models, Reasoning and Inference-（邦題：統計的因果推論 -モデル・推論・推測-）という本が有名らしい．

## 時系列因果推定
時系列データ間の因果推定は時間は原因が結果に先行することを利用して因果性の判断を行うので、上記に述べた手法とはやや異なる枠組みを採用する。


### Granger Causality(Granger et al.,1969)
おそらく時系列データの因果関係の推定といったときに代表的な手法はこれであろう。
これのコアのアイデアとなるのは**YがXから影響を受けているならば、Xの情報を使うことでYの予測精度が向上するはずである**というアイデアである。
Grangerの原論文[^1]はこのうち**予測**に対して**VARモデル**を用いてる。
VARモデルとは以下のような形で時系列予測を行うモデルである
$y_t = c + \phi_1 y_(t-1) + \ldots + \phi_py_(t-p) + \epsilon_t, \epsilon_t \sim W.N.(\Sigma)$
要するに次の時間の値は、p時間前までの値の線形な和＋ノイズによって書くことができるというモデルである。
具体的なアルゴリズムとしては沖本(2010)pp.80などを参照するとよい。(http://tjo.hatenablog.com/entry/2013/07/30/191853などにも引用されている)

## Transfer Entropy
情報を使って予測精度を上げられるか，というと情報論的なアプローチが可能ではないかと思う人もいるはずである．実は定式化の一つに情報論的なものもある．まず平均情報量の定義は以下の通り(面倒なので離散的な定義)

### 平均情報量(Entropy)

```math
H(X) = -\sum_X P(X)\log_2{P(X)}
```
これはXという変数の値を知った時に得られる情報量の平均を表す．

続いて，二変数間の情報量を考えるために相対エントロピーと相互情報量の定義を復習する．
相対エントロピーはYを知ったうえでのXの情報量を表す．

### 条件付エントロピー(Conditional Entropy)

```math
\begin{eqnarray}
H(X|Y) &=& -\sum_Y P(Y)\sum_X P(X|Y) \log_2 P(X|Y)\\
&=&-\sum_Y \sum_X P(X,Y) \log_2 P(X|Y)
\end{eqnarray}
```

### 相互情報量(Mutual Information)
相互情報量とはXを知る上でYという情報がどれほど役に立つかを測る尺度．Xの情報量から，Yを知った上でのXの情報量を引くことで測ることができる．相互情報量は以下のように計算できる．

```math
\begin{eqnarray}
I(X;Y) &=& H(X)- H(X|Y)\\
&=& -\sum_Y\sum_X P(X,Y)\log_2{P(X)} + \sum_Y \sum_X P(X,Y) \log_2 P(X|Y)\\
 &=&  \sum_X \sum_Y P(X,Y)\log{\frac{P(X,Y)}{P(X)P(Y)}}
\end{eqnarray}
```
式を見ればわかるように$X,Y$が独立ならば相互情報量は0である．相互情報量は式からもわかるように対称な式であることに注意しよう．このままでは因果の方向性を考えることはできない．

### 条件付相互情報量(Conditional Mutual Information)
Zを知った上でYを知ったときに、Xについて得られる情報量。

```math
\begin{eqnarray}
I(X;Y|Z) &=& H(X|Z) - H(X|Y,Z) \\
&=& H(X|Z) + H(Y|Z) - H(X,Y|Z)
\end{eqnarray}
```
例えば，yがXの未来予測にどれぐらい役に立つか(yがxへの因果を持つか)を知りたいときは，$I(X(t+1);Y(t)|X(t-1))$などを考えればよいことがわかる．

### Transfer Entropy
そしてTransfer Entropyは条件付相互情報量を用いて以下のように定義される．

```math
T_{X\to Y} = I(Y_{t+1};X_t^{(k)}|Y_t^{(l)})\\
```
ただし$X_t^{(k)} = (X_t,X_{t-1},\ldots,X_{t-k+1})$などを意味するとする．
これが有意に0より大きいのであれば，ある種の因果関係があると考えれらる．この手法はShreiber[^3] によって2000年に提唱された．
GCの自然の拡張になっているが，実際に計算するには，同時確率分布を推定する必要があり，そこに難しさがあるだろう．

### Causation Entropy
Transfer Entropyは直接の影響と間接的な影響を区別できない．そこでCausation Entropy[^4]というのも改善策の一つとして挙げられているようだ．Causation Entropyは

```math
C_{X\to Y} = I(Y_{t+1};X_t^{(k)}|Y_t^{(l)},Z_t^{(l)})\\
```
のように書ける量で，外部変数の影響を考えることができる．

### 幾何学的統合情報量
Transfer Entropyなども内包する．系の因果的影響を導出するための手法[^8]. Transfer Entropyなどをこの量から導出できる．詳細については[理研のプレスリリース](http://www.riken.jp/pr/press/2016/20161207_1/)などが参考になる．

```math
\Phi(x,y) = \min_{q\in\mathcal{M}_D} D_{KL}(p(x,y)||q(x,y))
```



## Transfer Entropyに対する批判
Transfer Entropyは2000年に提案された手法でありカオスなどの複雑系を扱う分野において広く用いられつつあるようである．しかしこの手法に対する改善策，批判も存在する．
Transfer Entropyに対する批判も上がっていたので紹介する．[^5][^7]同時確率分布を経験分布であらわす難しさ，直接の影響と間接的な影響を区別できないこと，情報源を特定できないことなどが批判の理由であるようだ．
例えば以下のようなサンプルはTransfer Entropyで正しく因果関係をとらえられない．[^5]
$X_t:0,1$をランダムにとる変数．$Y_t:0,1$をランダムにとる変数としたときに，$Z_t = X_{t-1} \mathrm{XOR} Y_{t-1} $と定めると，$X_t,Y_t$はともに$Z_t$を定めるのに必要な情報である(=因果性を持つ)のにもかかわらず$T_{X\to Z} = T_{Y\to Z} = 0$になる．




# 決定論的な系における因果性
今まではある種の確率論的な系を扱ってきた
以下のような系を考えてみよう．

```math
X(t+1) = r_X X(t)[1-X(t) - \beta_{XY} Y(t)]\\
Y(t+1) = r_Y Y(t)[1-Y(t) - \beta_{YX} Xt)]

```
これは漸化式として解くことができて,

```math
X(t) = \frac{r_X}{\beta_{YX}} \left\{(1- \beta_{XY} Y(t-1))\left(1 - Y(t-1)- \frac{Y(t)}{r_YY(t-1)}\right) -\frac{1}{\beta_{YX}}\left(1-Y(t-1)-\frac{Y(t)}{r_Y Y(t-1)}\right) \right\} \\
Y(t) = \frac{r_Y}{\beta_{XY}} \left\{(1- \beta_{YX} X(t-1))\left(1 - X(t-1)- \frac{X(t)}{r_X X(t-1)}\right) -\frac{1}{\beta_{XY}}\left(1-X(t-1)-\frac{X(t)}{r_X X(t-1)}\right) \right\} \\
```
のようになる．すなわち，$X(t)$の情報のみから$Y(t)$を推定することは可能であり，予測改善性による因果推定が決定論的な系では成り立たないのではないかということを示唆している.これは，Starkの埋め込み定理からも示唆される．

# Convergent Cross Mapping
以上のような決定論的な系における因果関係の特殊性を基に，SugiharaらはConvergent Cross Mapping(CCM)という手法を提案した．[^2]そして、$X\to Y$の関係がある時、むしろXの情報がYに蓄積されるので、**Yを用いることでXを推定できる**と考える。これは以下の
## 遅れ座標
CCMはカオス時系列解析の流れを汲んでいる手法である．カオス時系列解析の重要な解析手法の一つとして遅れ座標がある．遅れ座標とは時系列データ
$X(1),X(2),\ldots,X(N)$があった時に，$x(t) = [X(t),X(t+\tau),X(t+2\tau),\ldots,X(t+(E-1)\tau)]$のように生成されるベクトルの列である．これは実は力学系のアトラクターに対して色々と良い性質を持っており，カオス時系列解析によく使われる．ここでは以下のTakensの埋め込み定理とStarkの埋め込み定理を紹介する．
## Takensの埋め込み定理

### delay embedding map

```math
\Phi_{f,\phi}: M\to\mathbb{R}^dを以下のように定める.\\
\Phi_{f,\phi}(x) = (\phi(x),\phi(f(x)),\ldots,\phi(f^{d-1}(x)))
```
ただし

### Takensの埋め込み定理

```math
Mをm次元のコンパクトな多様体とする．d\geq 2m+1のとき，関数の組(f,\phi)は開で稠密な集合である．
```

これは遅れ座標が埋め込みとなっていることを示す定理である．
Takensの埋め込み定理を視覚的によく表すのが良く知られているローレンツアトラクタの遅れ座標による再構成であろう．

## Starkの埋め込み定理
以下のような"歪んだ系"を考える．

```math
x_{t+1} = f(x_t,y_t)\\
y_{t+1} = g(y_t)
```
Starkの埋め込み定理は，大雑把に言えば，xの観測データから，力学系全体の情報を復元できることを意味している．

### delay embedding map(2)
```math
\Phi_{f,\phi}: M\times N\to\mathbb{R}^dを以下のように定める.\\
\Phi_{f,\phi}(x) = (\phi(f^{(0)}(x,y)),\phi(f^{(1)}(x)),\ldots,\phi(f^{(d-1)}(x,y)))
```
このような観測時系列の系列が，元の力学系の埋め込みになっていることがStarkの埋め込み定理の主張である．
### Starkの埋め込み定理

## アルゴリズム
このアルゴリズムは非線形時系列予測法,シンプレックス投影法[^6]　に基づいている．
### simplex projection
シンプレックス投影法のアイデアは，アトラクタの次の点を，遅れ座標において近傍にある点から
なおこの手法は予測を目的としているのではなく，予測を通じてカオス時系列とランダムな時系列を区別するために考案されていることに注意したい．



### CCMのアルゴリズム
要点をまとめると，
1. 使うアトラクタの点Lを決める．
2. シンプレックス投影法に基づいてXからYを推定
3. 推定した時系列について本来の時系列との相関係数をとる
4. 1~3を繰り返すLを増やしながら繰り返す
5. 相関係数が上昇していけばYがXの原因となっていると考える．
という手法である．

一応Githubで実装は後悔しているが，自前の実験がらみのモジュールがごちゃ混ぜになっているので使うのが大変かも．
https://github.com/caprest/nonlinear_lab


このアルゴリズムをナイーブに実装すると計算量が$O(n^2\log{n})$となる、$n\sim10^3 $程度なので普通に重い。
なお検定方法は，IAAFTサロゲートを使うようだ．

### 具体例
もう一度下の系を考える.

```math
X(t+1) = r_X X(t)[1-X(t) - \beta_{XY} Y(t)]\\
Y(t+1) = r_Y Y(t)[1-Y(t) - \beta_{YX} X(t)]
```
この系をcoupled Logistic Mapと呼ぶ．


# CCMの発展手法について
* 検定についてはtwin-surrogateを用いる方法が中山(2015)に提案されている。因果関係を断ち切った時期列データを大量に生成することで、相関係数が優位に高くなっていることを示す。


# Misc
* ぶっちゃけ中山ら(2015)を読んだ方が良く理解できると思います。
* ここら辺のトピックを幅広く扱った[因果フェス](http://takehiko-i-hayashi.hatenablog.com/entry/2015/06/12/130826)というのがあったらしい。名前がかっこいい。
* 最近GCはDeepと絡めた手法が提案されたらしい．応用の余地はありそうだ．
* http://toyoizumilab.brain.riken.jp/taro/papers/tajima14biolsci.pdf 何かもわかりやすい


# 参考資料
沖本竜義. 経済・ファイナンスデータの計量時系列分析. 朝倉書店, 2010.
[^1]: Granger, Clive WJ. "Investigating causal relations by econometric models and cross-spectral methods." Econometrica: Journal of the Econometric Society (1969): 424-438.
[^2]: Sugihara, George, et al. "Detecting causality in complex ecosystems." science 338.6106 (2012): 496-500.
[^3]: Schreiber, T. (2000). Measuring information transfer. Physical review letters, 85(2), 461.
[^4]: Sun, J., & Bollt, E. M. (2014). Causation entropy identifies indirect influences, dominance of neighbors and anticipatory couplings. Physica D: Nonlinear Phenomena, 267, 49-57.
[^5]: James, R. G., Barnett, N., & Crutchfield, J. P. (2016). Information flows? A critique of transfer entropies. Physical review letters, 116(23), 238701.
[^6]: Sugihara, G., & May, R. M. (1990). Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature, 344(6268), 734.
[^7]: Smirnov, D. A. (2013). Spurious causalities with transfer entropy. Physical Review E, 87(4), 042917.
[^8]: Oizumi, M., Tsuchiya, N., & Amari, S. I. (2016). Unified framework for information integration based on information geometry. Proceedings of the National Academy of Sciences, 113(51), 14817-14822.