しげ日記

« セミナ | メイン | 図書館 »

2006年01月20日

JNS 1/19

Spike Count Reliability and the Poisson Hypothesis
Asohan Amarasingham, Ting-Li Chen, Stuart Geman, Matthew T. Harrison, and David L. Sheinberg

First author のアソハンは、今はうちのラボのポスドクで、実は席が僕のとなり。数学出身で、専門は確率統計みたいです。僕は統計数理は苦手なので普段はこの手の論文はパスなのですが、pooneil さんのとこでも紹介されてますし、何より席がとなりなのでこの論文はさすがに読まないとまずいですね。。（笑）

ストーリーを要約すると、

(1) ニューロン発火のスパイク列は、統計的に分散が大きいので、Poisson process ではないかと考えられていた。
(2) しかし、今までの統計方法では、分散を大きく見積もりすぎている。
(3) そこで、Poisson 過程であるかどうかを判定する新しい統計手法を開発した。
(4) この統計手法で、実際のサル電気生理から得られたスパイク列のデータを解析すると、Poisson 過程よりも reliable なものであることが分かった

という内容。

では、まず Poisson 分布について。

たとえば、あるニューロンの発火について考えてみます。

ある短い time bin ΔT（例えば 1ms）で、ニューロンが発火する確率 p_ΔT は非常に小さい（例えば p_ΔT = 0.002）が、観測時間 n×ΔT が長い（例えば n = 1000）と、期待値 p×n はあるそこそこの値をもつ。（ここでは p×n＝0.002×1000=2、つまり 2 Hz）。このニューロンの発火を観測した場合に、観測された発火数 が例えば 3 である確率は、p_ΔT が一定とするならば二項分布に従うので、
₁₀₀₀C₃ 0.002³(1-0.002)⁹⁹⁷。この計算を直接実行するのは大変である。しかし、p が非常に小さく、n が非常に大きい場合、観測される発火数 m の確率を求めるのには以下の定理が使える（ポアソンの小数の法則）；

　_nC_m p^m(1-p)^n-m → e^-λλ^m/m!

ここでλは期待値（λ=np）。これが Poisson 分布で、つまり Poisson 分布は二項分布の n → ∞、p → 0 の極限での分布のこと。（また、λを大きくしていくと、ポアソン分布は正規分布に近づく。）

この式から上の確率を求めると、λ=2,m=3より、e^-λλ^m/m!=0.18。

（以上の話は「統計学入門」(東大出版会）より改変）。

ポイントは、微小 time bin ΔT における発火確率 p_ΔT が、常に一定であるというのが、Poisson 分布を持つための必要条件である。

ところで、現実のニューロンでは、p_ΔT は、長い時間のスケールで見ると、一定でないように見える。ある event があれば発火率は上昇するだろうし、何もなければ発火率は減少するかもしれない。つまり、発火率 λ は時間によって変化する（すなわち、λ=λ(t)、p_ΔT = p_ΔT(t) ）。ただし、λや p_ΔT が大きなスケール時間的に変動しても、局所的（比較的短い時間スケール）でみれば、 p_ΔT は、すべての time bin に対して等価である、としたのが、inhomogenous Poisson process 。（このへんの話は、Shadlen & Newsome 1999 のあたりで提案された話らしいのですが、僕はあの有名論文、難解なためにいつも途中で挫折してしまっていまだ理解できておりません。。）

つまり、Poisson Process でない、ということは、局所的に見ても、p_ΔT は、すべての time bin に対して等価でないということであります。

　　　　＊　　　　　　　　　＊　　　　　　　　　＊

ここで、本題に入る前に、以後使用する記号をまとめておきます：
　i 番目の trial で観測された発火率　　　　　 m_i
　n 回の trial で観測された発火率の平均　　μ'= (1/n)Σm_i
　n 回の trial で観測された発火率の分散　　σ'²
　真の発火率（発火率の期待値）　　　　　　　　 λ
　真の分散　　　　　　　　　　　　　　　　　　　　　 σ²
（真の値は unknown なので、観測によって得られた値で代用する、ということです）

さて、ニューロンの spike train が Poisson process に従う、という根拠は、「trial毎の発火率の分散が大きい」という実験事実があるためです。

Poisson 分布では、その数学的性質上、期待値と分散は等しくなります。つまり、Poisson 分布では、Fano factor = σ²/λ = 1。

Shadlen 論文での主張は、観測によって得られた Fano Factor（σ'²/μ' ）が 1～1.5 になるので、spike train が Poisson prcess であるとの仮説が成り立つのではないか、というものです。

それで、ここからが本題になるのですが、以上の Shadlen 論文の主張に異議を述べたのが、本論文です。

このような spike 列の統計的な処理を行うときには、「trial毎の観測の揺らぎは存在しない（trail-to-trial statistical stationarity）」という暗黙の仮定をおいています。

つまり、観測による trial毎の揺らぎを考慮していません。しかし、もし、trial毎に揺らいでいるのであれば、Fanofactor（＝分散/平均）が大きい値を持つのは当然、ということになります。

ここで「観測の揺らぎ」と言っているのは、観測自体に内包しているしているような揺らぎのことです。さいころを振るような trial の場合は、tiral毎の観測の揺らぎはないと考えられます。しかし、サルの電気生理のような場合、trail毎の条件が同じであったとしても、サルの集中力や疲れなどの hidden な条件によって揺らぎが生じるので、trail毎の観測の揺らぎが乗ってしまう可能性があります。

つまり、trial毎に揺らぎが考えられる観測では、単に σ'²/μ' が 1～1.5 になるからといってそれが Poisson Process であることの主張とするには根拠が弱い、という主張です。（なぜならσ'にはニューロン発火の純粋な揺らぎ（真の分散）に加えて、サルの集中力といった観測の揺らぎといういらないものが入ってしまうから。）

それで、trail 毎に揺らぎがある場合でも Poisson Process であるかどうかを判定する方法を提案したのが、本論文なわけです。

ではどうやるか。

この論文での戦略として、帰無仮説を立ててそれを棄却できるかどうかで判定するという手法を用います。

帰無仮説（null hypothesis）H₀

　m₁, m₂, ・・・, m_n
　は独立な Poisson random value である

で、この帰無仮説を棄却できれば、Poisson 過程ではないと言える。（棄却できなければ Poisson 過程であるともないとも言えない）。

で、その判定として、以下の不等式を用いる；

　Σ m_i² ≦ f　（論文中の式(2)）

つまり、帰無仮説を棄却できる値 f を見つけることができれば、統計的に Poisson 過程でないかどうかの判定が可能になる、というわけです。それで、この値 f を数学的に見つけてやろう、というのがこの論文の前半での主題。

f は、trail回数 n、観測平均 μ'、および棄却のレベル（どれくらいの確率で棄却できるかということ）αの関数になると考えられる。すなわち、f = f(n,α,μ')。

それで、結局、論文中の式（10）によって f を求めることができる、というのが本論文での結果です。

（式(10)の導出方法の詳細は省略しますが、Trial-to-trial statonarity を仮定した場合との統計的な比較によるみたいです。 ）

数値的には、モンテカルロシミュレーションを使って、f を求めるみたいです。

一応、ホームページに f を求める Matlab プログラムを置いているみたいなので、そのプログラムを使えば、自分のデータを試してみることも可能です。

以上の話をまとめると、
(0) 今までは、Poisson process かどうかの判定に、Fano factor（分散/平均）が～１であるかどうかを使っていた。
(1) しかし、Trial-to-trial に揺らぎ（観測に依る揺らぎ）がある場合は、Fano factor が真の値より大きくなってしまうという不都合がある。
(2) それで、新しい統計手法を開発した。ついでに、簡単に使えるようプログラム化した。
ということです。

それで、論文の後半では、実際のサル電気生理での spike 列のデータを使って、Poisson 過程であるか否かを解析しています。結果を簡単にいうと、stimulus onset から400msぐらいまでの区間（発火率が比較的高い区間）では、スパイク列は Poission 過程ではない、（Poisson よりもっと規則正しい）ことが統計的に示されています（Figure4b）。Stim onset から 500 ms以上離れてくると、Poisson 過程でないと言えなくなってきます（帰無仮説が棄却されないため）。

　　　　＊　　　　　　　　　＊　　　　　　　　　＊

さて、（上で述べたことを繰り返しますが、）Poisson Process でない、ということは、局所的に見ても、p_ΔT は、すべての time bin に対して等価でないということです。

端的に言えば 、ニューロンには、refractory effect（1回発火した直後にはすぐ(～2ms)には発火できない特性）や burst などの biophysical な特性があるが、発火率が上昇すれば上昇するほど、そういう特性に支配されて、ランダム性（Poisson 性）が減少していく。まあ、当たり前といえば当たり前の話ですが、統計数理的に、どのレベルで Poisson 的であるかそうでないかを議論したのが、この論文の新規性ということでしょう。

投稿者 sfujisawa : 2006年01月20日 19:42

コメント

投稿者 pooneil : 2006年01月26日 00:32

投稿者 しげ : 2006年01月28日 00:29

コメントしてください