« August 2008 | Main | November 2008 »

Panzeri-Treves '96について

Panzeri and Treves, Analytical estimates of limited sampling biases in different information measures, Network: Computation in Neural Systems 7 (1996) 87.

を読むときには気をつけたほうがいい。ちょっとはまりポイントがあって読むのに時間がかかった。それについて書く。これが他の人の助けになればいいな。あとは連続カーネルでのsmoothingに興味が無ければそこはすっ飛ばしても大丈夫というのがアドバイス。

ちなみにこの論文はレートコーディングの情報量[1]のupwardbiasの補正公式を与えるものだ。


問題は式(A4)にある。
  〈pNk(i|s)〉= pk(i|s) + 1/Ns kC2pk-2(i|s) [q(i|s) - p2(i|s)] + higher order
(オリジナルの式ではp, qに皆チルダがついているけど、webでは出せないから略した。またp(i|s) ≡〈pN(i|s)〉。)
ここでqの定義は(18)式にσ=0を代入して
  q(i|s) = ∫mi-1 midr P(r|s) = p(i|s)
となる。ここで mi はi番目のbinの上限である。けど、これはおかしい。ここがハマリポイント。本当は一般論としては(A4)の式は
  〈pNk(i|s)〉 = pk(i|s) + kC2pk-2(i|s) [〈pN2(i|s)〉- p2(i|s)] + higher order
であるべきはず。

じゃ、なんで大体正しい結果が出たかというとたまたまビンの占有数はポアソン分布だから
  〈(n(i,s))2〉-〈n(i,s)〉2 = 分散 =〈n(i,s)〉
   pN(i|s) ≡ n(i,s)/Ns
  〈(pN(i|s))2〉- p(i|s)2 = p(i|s)/Ns
だから。この右辺がq(i|s)にたまたま一致しちゃう。ここでn(i,s)は問題のビンの観測度数。p(i|s)2/Nsに相当する部分はまた別なところから出てくるようだけどその議論は割愛。


あともう一つ要注意なのは式(A1)の和Σの上のカレット。これは(A1)の次の3行目に書いてあるとおり期待確率 p(i|s)がゼロのビンを取り除くことを表している。しかし、補正公式を実際に使っている場面では、これを観測度数がゼロのビンを取り除く意味に取っている人がいる。というか私も最初、そう勘違いしました。例えば期待度数が1ならばポアソン分布なので観測度数がゼロになるのはむしろきわめてよく起こることであるけど、そういう時にこのビンを取り除いてはいけない。

で、期待度数λとすると私の計算では
  Σkpoisson pdf(k, λ) k log k -λlogλ
が補正値なのだけど、期待度数が0.3ですら補正 0.42を与えて、λ= ∞の補正値1/2とさして違わない[2]。

しかし現実問題として期待度数が0.3でビンが空になっているのか、期待度数が0でビンが空になっているのか実験から決めることはほとんど不可能だと思う。この点でかなり不便な公式ではある。

期待度数が0.1でもまだ補正は0.23だからいっそのこと期待度数に関わらず1/2という式の方が良いかも。実際、そこにbinがあると言うことは、ほかの刺激では占有されていると言うことなのだから、その刺激でも0.1くらい期待度数があると考える方が自然だと思う、ニューロンの反応のノイジーさを考えると。実際自分の持っているニューロンの実データでもその方が良く合った。(合ったというのは刺激提示前の発火つまり情報量ゼロのはずのところの結果がどれだけゼロに近いか。)

しかし期待度数に関わらずという事にすると結局文献[3]に一致する。

[1] Optican L.M. and Richmond B.J., Temporal Encoding of Two-Dimensional Patterns by Single Units in Primate Inferior Temporal Cortex. III. Information Theoretic Analysis, J. Neurophys. 57 (1987) p.162

[2] 蛇足だけどこの式はλ=1だと0.57くらいになり1/2より大きい。実際その方が補正が正確だった。

[3] Treves and Panzeri, The Upward bias in Measures of Informaiton Derived from Limited Data Samples, Neural Comput. 7 (1995) 399

| | Comments (0) | TrackBack (0)

Rate codingしても別に情報量は大きくならない、というか小さくなる

ニューロンのスパイク活動がコードしている情報量についていまのところ一般的な定義はレートコーディングによるものです。つまり、ある時間windowを決めたうえでのspike数(=レート)と刺激の間の相互情報量を求めています。この定義では文献[1]がよくciteされていますが最初に言い出したのが誰かはしりません。

一方私は、小さなwindowでspike数(0 or 1)と刺激のあいだの相互情報量を取ったほうが良いのではないかと考えています。そのことについては細胞が運ぶ情報量の計算法について (PDF)に書きました。この考え方をとりあえず「スパイクコーディング」と呼んでおきます。先行文献があるかどうかは分かりませんけどあるんじゃないかな。

両者の違いは、情報量についての考えの違いではなく、コーディングに対する考えの違いです。前者は要するに発火率(レート)や発火数が情報をコードしていると考えているのに対して、後者はスパイクが情報について語っているという考えです。そのために前者には時間windowとbinningというartefactが必要になります。また前者の方が多くの実験トライアルを必要とするという欠点もあります。


さて今回はさらにrate codingは引き出せる情報量がよりすくないということを論じます。

一見rete codingの方が多くの情報が得られそうに思えます。例えば刺激が10個ある時、binが10個あればrate codingは刺激を完全に特定出来る可能性があります。(rate codingは多段階のgradientがあり中間調を表現できます。)一方spike codingの方はspikeは1/0なのですからスパイク一つ一つを見ているだけでは特定は出来ません。

しかし、実はそうではありません。

実際に計算してみます。刺激1には10Hz、刺激2には20Hz、...刺激10には平均100Hzで反応するようなニューロンを考えます。ここで最初の100msに運ばれる情報量を考えます。

このとき、spike codingでは情報量が12.0 bit/sになります。手で数式をいじるより計算する方が早いのでmatlabコードを示します。このとき100msあたりでは1.20bitとなります。

p2 = 10:10:100;
dt = 1e-6;
p = [p2*dt; 1-p2*dt]/10;
(sum(sum( plog(p) )) - sum( plog(sum(p,1)) ) - sum( plog(sum(p,2)) ))/log(2)/dt
function y = plog(p)
y = p.*log(p + 1e-38);


一方、rate codingの定義だと、100msの間に刺激1はスパイク1個、刺激2はspike2個 ...となりますから、10種類の刺激を弁別出来て log 10/log 2 = 3.32 bitとなりより大きな情報量となるように見えます。

しかし、ニューロンがこのように規則正しく発火することはありません。実際は発火数はランダムにばらつきます。その発火はポアソン過程っぽい感じです。でもポアソン過程であれば Fano factor = (分散/平均) は1ですが、ニューロンのスパイクは実はこれが1よりも大きく2に近い値が本当です。つまりポアソン過程に似ているけどそれよりもランダム性が強い発火をするわけです。

このランダム性を取り入れて計算してみます。しかしrate codingの定義に思いっきり有利に計らってFano factor = 1のポアソン過程にした上で、計算するプログラムが次です。

X = [];
for i=1:10
   X = [X; poisspdf(0:19, i)];
end
p = X/10;
(sum(sum( plog(p) )) - sum( plog(sum(p,1)) ) - sum( plog(sum(p,2)) )) /log(2)

結果は0.72 bitとなり、spike codingの値より大幅に小さくなります。

つまりあの定義が勝つためにはニューロンの揺らぎをうち消すような厳格な符号化が必要なわけです。

おそらくですが、ポアソン過程であるという仮定をおけば常にspike coding定義の方がよいと証明できるのではないかと思います。だとすれば、graded responseを利用できるとか、rateコーディングの可能性とかそういう(これまでこの分野にpersistしてきた)観念は捨た方が良いのではないかと思います。

[1] Optican L.M. and Richmond B.J., Temporal Encoding of Two-Dimensional Patterns by Single Units in Primate Inferior Temporal Cortex. III. Information Theoretic Analysis, J. Neurophys. 57 (1987) p.162


スパイクは 1 or 0 でも中間の刺激の弁別情報を伝えることが出来る

さて、以上で論じたように、スパイクコーディングでは刺激がリッチに10個もあったとしても、神経細胞ができる反応は0か1かの2通りしかありませんから、何段階もの反応を返すと考えるレートコーディングよりも識別力が弱いのではないかと考える人もいるかと思います。しかしそうではありません。

Ex1刺激が10個あるとします。刺激1には10Hz、刺激2には20Hz、...刺激10には平均100Hzのポアソン過程で反応するようなニューロンをまた考えます。この発火確率が分かっている状態で、実際に観測されたスパイクからどの刺激が与えられたかデコードするとします。

実際に一つの合成されたスパイク列から、それぞれの刺激の尤度を計算してそのシェアを示したものが右上の図です。縦軸が確率で横軸が時間[ms]となっています。最初はどの刺激であるか分かりませんからみな確率 1/10としてスタートします。スパイクが来るたびに曲線が不連続になっています。刺激1,2などの低い周波数しか生まない刺激は、スパイクがこないでいる区間でだんだんと尤度を上げていきます。反対に9,10などの高い周波数の刺激は、スパイクが来た瞬間に尤度があがります。

実はこの図は30Hzのポアソン過程でスパイクを生成したものですが、右に行くにしたがって実際に30Hzに対応する刺激3という仮説の尤度が上がってグラフの上に出てきています。これはスパイクがないときの低い刺激の増加傾向とスパイクが来たときの高い刺激の増加傾向の丁度釣り合ったところの刺激の尤度が出てきているわけです。

右下の図は分かりやすいようにポアソン過程でなく規則的な30Hzの過程にして典型が見えるようにしたものです。Ex2

このようにして30Hzのような中間の刺激でも判別しうることが分かります。

これに使ったMATLABプログラムは以下の通りです。

p = (10:10:100)' *0.001;
g = p*0 + 0.1;
h = [];
nspk = 0;
for i= 1:1000
   if(rand(1) < p(3)) g = g.*p/mean(p); nspk = nspk +1; else g = g.*(1-p)/mean(1-p); end
   %if(mod(i,round(1/p(3)))==0) g = g.*p/mean(p); nspk = nspk +1; else g = g.*(1-p)/mean(1-p); end
   g = g/ sum(g);
   h = [h g];
end
plot(h','LineWidth',1)
xlabel('time [ms]')

| | Comments (1) | TrackBack (0)

« August 2008 | Main | November 2008 »