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)にある。
　　〈p_N^k(i|s)〉= p^k(i|s) + 1/N_{s k}C₂p^k-2(i|s) [q(i|s) - p²(i|s)] + higher order
（オリジナルの式ではp, qに皆チルダがついているけど、webでは出せないから略した。またp(i|s) ≡〈p_N(i|s)〉。）
ここでqの定義は(18)式にσ=0を代入して
　　q(i|s) = ∫_mi-1 ^m_idr P(r|s) = p(i|s)
となる。ここで mi はi番目のbinの上限である。けど、これはおかしい。ここがハマリポイント。本当は一般論としては(A4)の式は
　　〈p_N^k(i|s)〉 = p^k(i|s) + _kC₂p^k-2(i|s) [〈p_N²(i|s)〉- p²(i|s)] + higher order
であるべきはず。

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

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

で、期待度数λとすると私の計算では
　　Σ_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

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なのですからスパイク一つ一つを見ているだけでは特定は出来ません。

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

実際に計算してみます。刺激１には10Hz、刺激２には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個、刺激２はspike２個 ...となりますから、10種類の刺激を弁別出来て log 10/log 2 = 3.32 bitとなりより大きな情報量となるように見えます。

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

このランダム性を取り入れて計算してみます。しかし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個あるとします。刺激１には10Hz、刺激２には20Hz、...刺激10には平均100Hzのポアソン過程で反応するようなニューロンをまた考えます。この発火確率が分かっている状態で、実際に観測されたスパイクからどの刺激が与えられたかデコードするとします。

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

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

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

このようにして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]')

自分のOSXマシンにWikiをインストールしてみた

自分のマシンはOSX 10.4.11でApache 1.3.41が含まれている。
Apacheはコントロールパネルの「Web共有」から起動できる。

1) hiki (CGI型のwiki)

hiki-0.8.7をインストールした。OSX 10.4.11にはRuby 1.8.2が含まれている。

インストールは

Hikiのインストールと設定(http://hikiwiki.org/ja/install.html)の通りでOKだった。

ただApacheをCGI向けに設定しなくてはならない。

Max OS X で Apache サーバを立てよう - Apacheの設定 -(futomi's CGI cafe)
を参考にしたのだけど、/private/etc/httpd/httpd.confをいじる必要は結果としてなかった。

/private/etc/httpd/users/nnnn.confを以下のようにしたら充分だった。効果の範囲が自分だけだしファイルも小さいのでやりやすい。

< Directory "/Users/nnnn/Sites/">
    Options Indexes MultiViews
    AllowOverride None
    Order allow,deny
    Allow from all
</Directory>
< Directory "/Users/nnnn/Sites/hiki/">
    AllowOverride All
    Options +ExecCGI
    AddHandler cgi-script .cgi
< /Directory>

nnnnは自分の名前。６行目まではもとのデフォルトのまま。

最初は、先ほどのページをみても検索でかかったほかのOSX Apacheについてのページをみてもうまくいかなくて、はまった。けど、その理由はAllowOverrideの範囲が足りないからのようだ。
AllowOverride Allとしたら動いたけどこれはやりすぎかも知れないので責任が持てない。hikiの.htaccessの中身をみて最小にした方がよいのかもしれない。

あと上のページにOSXのシェルに慣れていない人のためのガイドがいろいろある。

2) pukiwiki (mod_php型のwiki)

pukiwiki-1.3.7をインストールした。OSXに含まれていたのはPHP 4.4.8。

pukiwikiの公式ページからダウンロードして展開したファイルに含まれるreadme.txt通りの手順でインストール。私の所では実質chmod以外はする必要がなかった。

あとは/private/etc/httpd/httpd.confの240行付近

#LoadModule php4_module        libexec/httpd/libphp4.so

と280行付近

#AddModule mod_php4.c

の#をコメントアウトしたらOKだった。あとのことは標準のhttpd.confの後ろの部分が面倒を見てくれる。

ブラックホール解の内部解の非一意性？

ブラックホールに吸い込まれると、主観視点では有限の固有時間ののちspacelike singularityにぶつかると言われている。これはSchwarzschildの外部解をKruskal座標を使って延長して求められた結果である。しかし、このspacelike singularityというのは物理的な解釈としては変、というかありえ無さそうに私には見える。

あと、Kruskal解のペンローズ図はもとの宇宙が右側だとしたら、左側から別の宇宙が繋がっている。しかしこれも物理学としては不健全な感じがする。

しかし微分方程式の解は一意なのでこれで仕方がないと思っていた。でもそうでもないのではと思えてきた。

例えば普通のflat時空の電磁場などでよいけど、一点から出た波の解が一意というのは、他のlightconeからの影響がないという条件の話だ。同じようにKruskal座標のブラックホール内部の部分は、外部からみればこれも一つのlightconeであって、内部に重力のソースがあればこの部分の解は変更される。

言い直せば、事象の地平は内側から見てもlightconeであって、内側にソースがあれば解は変更しうる。実際ブラックホールでは内部にソースがあると考えた方が自然ではないだろうか？なぜならcollapsing starの表面をみながら観測者がブラックホールに飛び込んだとする。観測者の固有時間では有限時間で事象の地平を超える。そのとき、まだ目の前には小さな半径でcollapsing starの表面が見えているはずである。そしてそれは重力のソースであるはずだ。ならば内部解は変更されない方がおかしい。

相対論には詳しくないし、物理ってちょっと間違えることは良くあるんで、びくびくなんだけど。どうだろう。
一般相対論に詳しい人に教えを乞いたいです。

Higgs粒子はいつ捕まるか？

簡単なまとめ。

LEP：最大115GeV。未発見で終了。

Tevatron (Fermilab)：最大180GeV。運転中。
特に114〜130GeVと150〜170GeVに強いらしい。

LHC：最大1TeV。'07 11月にビーム開始。08年春にデータを取りはじめる予定。
ATLASは09-10年中に見つけるという予想。
magnetの設計ミスが発覚。08年春開始に遅れるかも。
最大の性能を発揮するまで数年は掛かるものと見られている。

ググッた情報以外は
SCIENCE 315, p1656
SCIENCE 316 (2007, 4/6) p31
読売新聞４月２９日（埼玉版？）

確率微分方程式超入門

どうも世の中の確率微分方程式の入門はハイブロウすぎると思う。経済学のいくつかの講義ノートを見たけど、ボレル集合族とか確率空間とか伊藤積分とかはやりすぎ。そういうのは微分を学ぶのにいきなりε-δ論法から入門するみたいなものだ[1]。物理屋風の方法論なら以下のように簡単に伊藤の補題も導きだされる。きちんとしたのはその後でよい。

さて、確率微分方程式とは連続時間でランダムな擾乱を受ける系を記述する微分方程式である。

例えばブラウン運動を考えてみよう。ブラウン運動では粒子はランダムに動く。そして試行をくり返せばその統計集団は正規分布に従う。そして分布の標準偏差は時間tの時には √t に比例する。これがブラウン運動の基本的な性質である。

ブラウン運動を微視的に見れば、それぞれの微小時間Δtごとに大きさ b √Δtの正規分布に従う動きをしている。一般の微分方程式であればΔtのあいだの粒子の動きΔXは、Δtの一次に比例する。それに対してブラウン運動では、√Δtに比例する。このブラウン運動の微小過程を
　　ΔX = b ΔZ
と書いてみる。ΔZは標準偏差√Δtの正規分布を表す確率変数である。するとXも確率変数となる。ここで異なる時間のΔZ(t)とΔZ(t')は独立であるのがポイントである。さらに一定速で水が流れているような時のブラウン運動
　　ΔX = a Δt + b ΔZ
なども考えられる。これがgeneralized Wiener過程である。流速 a は空間座標xに依存しても良いし、時間に依存してもよい。
　　ΔX = a(X, t) Δt + b ΔZ
（場所によって温度がちがうかも知れないし）ブラウン運動の大きさ自体も空間、時間に依存してよい。
　　ΔX = a(X, t) Δt + b(X, t) ΔZ
とここまで一般化したのが、伊藤過程である。

さてこのように確率変数Xの微小変化が記述されている時に確率変数 X²の微小変化を考えてみる。普通にチェインルールを使って
　　ΔX² = 2 X ΔX　　　　(1)
だろうか？ そうは行かない。単純なブラウン運動
　　ΔX = b ΔZ
のときtにおける分散は b t だからtに比例して大きくなる。(1)式は増大と減少が同じだから期待値としては0のままであるけど、実際はX²の期待値は t に比例して増大しなければならない。

一般の関数 f(X, t) で考えよう。ある時刻 t に座標x0に粒子が居てΔX = a + b ΔZ で運動するとする。時刻Δt では 粒子の x 座標は x0 + aΔt を中心とした分散 b √Δt の正規分布をとる。このとき、f(X, t)は中心 f(x0, t) + b²Δt f''(x0, t) で分散 (b f'(x0))²Δtの正規分布をとる。つまり非線形性の為に分布の中心がΔt に比例してずれている。ここで f' は f の X による一階微分である。ここでΔtの１次より上の次数の項は無視した。またΔt が充分小さければ、局所的には関数は線形なので、この分布は正規分布に充分近くなる。これ以外の項はΔt より高いオーダになるので効かない。

これを考えに入れると微小時間Δtの間にΔfも以下の式のように一般の伊藤過程で動く。
　　Δf = a (df(X,t)/dX) Δt + (df(X,t)/dt) Δt + b²(d²f(X,t)/dX²) Δt + b (df(x,t)/dX) ΔZ
第３項が普通のchain ruleで出てこない項であり、先ほどの中心のずれに対応する部分である。そしてこの式が有名な伊東の補題である。

以上の議論はΔtを半分に、そのまた半分に・・・とずっと小さくして行っても結果が変わらない。つまり極限が上手くいっているので、有限を示唆するΔではなくdと書いても良いだろう。だから以降は
　　df = a (df/dx) dt + (df/dt) dt + b²(d²f/dx²) dt + b (df/dX) dZ
のように書く。

ではこの確率微分方程式のご利益とはなんだろう。まず反対に確率微分方程式を使わなければ、どうなるか考えよう。ランダム揺動を受ける粒子の運動のような現象は確率分布p(x, t)に対する偏微分方程式で記述し解かなくてはならない。一般に確率微分方程式
　　dX = a(x,t) dt + b(x,t) dZ
に対応する偏微分方程式は
　　∂_tp(x, t) = (-∂_xa(x,t) + 1/2 ∂_x∂_x b²) p(x, t)
である[2]。p(x,t)は確率分布関数である。この形をした方程式をフォッカー・プランク方程式という。

そして確率微分方程式のご利益はいくつかの例を考えてみるとわかる。たとえば簡単な
　　dX = cX dZ
これは
　　d log X = - X^-2/2 (c X)² dt + X^-1c X dZ
　　　　= - c²/2 dt + c dZ
だから
　　log X = x0 - c² t/2 + c Z(t)
である。ここで Z(t) は分散 t の正規分布である。確率分布で書けば z = log x が
　　p(z, t) = (1/c√t) exp(-1)(z - c²t/2)²/c² t
という確率分布を持つ事がわかる。これをxに変数変換すると
　　f(x, t) = (1/c x√t) exp(-1)((log x - c²t/2)²/c²t)
である。一方これに対応するフォッカー・プランク方程式は
　　∂_tf(x, t) = 1/2 ∂_x∂_x(c x)² f(x, t)
で、解は上の f(x, t) なんだけど、解を見つけるのはなかなか面倒でだろう。必要な変数変換を思い付いたら良いけど思い付かないかも知れない。それにくらべて確率微分方程式はほぼ常微分方程式と同じく一本道で解く事ができる。

もう一つの例をやってみよう。
　　d X = - c X dt + a dZ
これはランジュバン方程式と呼ばれる。e^ctX に伊藤の補題を適用すると
　　d e^ctX = - c e^ctXdt + c e^ctX dt + a e^ctdZ = a e^ctdZ
となる。右辺は時間 t にそって刻々と標準偏差 a e^ctΔtのホワイトノイズが加わった過程であるので、時刻 t には
　　a²∫₀^te^2ctdt = a²/2c * (e^2ct-1)
の分散を持つ正規分布になる。だから
　　e^ctX - x0 = Z( a²/2c * (e^2ct-1))
が解である。一方対応するフォッカー・プランク方程式
　　∂_tf(x, t) = (-∂_xx + 1/2 ∂_x∂_xa²) f(x, t)
解は
　　f(x,t) = c^1/2(πa²(1-e^-2ct))^-1/2 exp(-c(x-x0e^-ct)²/a²(1-e^-2ct))
だそうで変数変換がかなりめんどくさそうである。（私も真面目にやる気がしない。）

このように、ある種の問題は偏微分方程式を使うよりも確率微分方程式を使った方が遥かに簡単に解くことができる。（流体力学にもオイラー描像とラグランジュ描像があったことを思い出させる。どっちかの形式で難しい問題も別な形式に移ると簡単だったりする。）

以上は物理学的な（つまり数学に比べると乱暴な）導出法であったけど、これを数学的にきちんとするには、まずは伊藤積分を定義してそれの逆として確率微分を定義する。伊藤積分は測度空間の一種である「確率空間」を使って定義される。確率空間は我々もよく知る「サンプルの空間」を非加算に拡張したようなものである。[3]

私はここは面白いところだと思う。例えばブラウン運動のような確率過程は実際に粒子があって、何度も観測を繰り返した結果を記述している。つまりサンプリングの結果である。しかし、物理で一般的な確率分布 p(x, t) による記述はこの粒子が全く見えなくなっている記述である。それに対して、確率微分方程式はこの粒子の追跡を基本としている。それが解の簡単さにつながり、また数学的厳密さにつながるということは面白い。

ランダム項が非常に小さくなると普通の運動方程式になる。一方、同じ過程を確率分布の偏微分方程式に直してしまうとその解き方は運動方程式とは全く違ってしまう。しかし、確率微分方程式はこの極限が見やすいという利点もある。

このように確率微分方程式は一つの役にたつ方法であり考え方である。

余談。株価もそれぞれ微小なランダム変動を繰り返している。ポートフォリオというのはそのような複数の株価の関数 f(x1, x2, ...) であるから、これに確率微分方程式と伊藤の補題を適用できる。経済ノーベルもどき賞を受賞したブラック=ショールズ微分方程式というのは、ただそれだけのことである。

[1] ε-δ論法を学ばなくて良いと言う意味ではない。最初は微小量の考え方で理解した上で、厳密な方法を知るべきという事。

[2] なんで b² が∂_x∂_xの右にあるかと、なぜ２乗なのか考えてみよう。またここで示したのは「伊東の方法」と呼ばれる流儀であるが、もう一つ「ストラトノビッチの方法」という流儀がある。（培風館物理学辞典を参照のこと）

[3] もう少し真面目な入門は確率微分方程式（PDF,by 千葉逸人）をどうぞ。これを読む前提はルベーグ積分を知っていること。ルベーグ積分の前提は位相空間。

テクスチャ合成

まず、右のイメージを見てほしい。囲みの中と外では違っている。

しかし、私の場合、まず遠目にちらっと見た時は違いが分からなかった。低解像度でも分からないかも知れない。つぎによく見ると、周りがランダムで中心が規則的だと思った。そして、ランダムな方が自然だと思った。最後に高解像度で良くみると、周りのほうが変なボケがある事に気がついて、それでこちらが合成であると気がついた。慣れてくると変な ねじれのような ゆがみが ある事にも気がつく。

この研究は99年くらいだから旧聞なんだけど、PortillaとSimoncelliらによるテクスチャ合成の研究により作られたものだ。つまり現実のテクスチャから特徴をパラメータとして抽出して、そのパラメータから合成されたのがこれらのイメージの周辺部分という訳である。彼ら自身による紹介は以下のページである。イメージもそこからとった。

Parametric Model for Visual Texture Representation and Synthesis

この合成の結果は実に見事だと思う。また、ゆがみ具合も、ちょうど夢で見るような感じ。あるいはこれなんかはゴッホに似ているような感じもする。（真っすぐさが捉えられて居ないようだ。つまり長距離の相関を捉え損なっている。）

まえに「image statisticsの研究はまだ現実のイメージとは何かが分からない」というような事を書いたけど、その見解は修正しなくてはならないようだ。つまり現実のイメージを再現する特徴の抽出はほぼ出来ている。

ただし、抽出されるパラメータは多岐に渡っていて、その総数は706個。そもそも合成されるイメージの一単位のピクセル数が同じくらいだから、意地悪く言えば、特徴を抽出しているというよりはイメージの特殊なストアの仕方に過ぎないと言う見方もできる。

Mind Hacks

Mind Hacks 実験で知る脳と心のシステム
Tom Stafford, Matt Webb, 夏目 大 (訳) (オライリージャパン, 2005)

読了。面白い本だと思う。

じつは認知的脳科学の最近の研究と言うのは、まさにこの本に書かれたようなことをやっている。科学の論文と言うのは正確さと間違いのなさを重視して書かれるから分かりにくくなるけど、正確さを無視して分かりやすく書くとこういう感じになるという本。だからこの分野について入門したい人には面白いと思う。

以下は、自分の備忘の為に：

▼hack #92 コーヒーの味はカップで変わる？

カフェインがドーパミンと競合する神経伝達物質（アデノシン）の働きを妨害するため、その結果としてドーパミンの働きを活発にすると言われている。(p365)

カフェインは脳の「報酬回路」を化学的にハッキングする。(p363)

これは知らなかった。カフェインは覚醒作用と言われるから、単に覚醒系に作用するのだと思っていた。

▼hack #85　偽物の記憶（p334）
この前後の節の偽の記憶に関わる話は面白い。

▼hack #76 共同運命の原理

一説には、情報の「結び付け」の際には、多数のニューロンが同じタイミングで発火するとも言われ、「共通運命の原理」はそのメカニズムの反映ではないかとも言われる。(p299)

bindingのsynchronizationメカニズムによるこの説明は、（私から見ると分の悪い）仮説にすぎない。しかし、それ以上にメカニズムが現象を説明するという考え方は、生物に対してはまずい考え方だ。現象が何らかの適応的意味をもち、それを実現するためのメカニズムがあると考えるべき。

▼hack #75 ゲシュタルト原理

近接
そばにあるものは自然に同じグループに属するとみなす事が多い。
類同
見た目が似ているものもひとまとめにする事が多い。
閉合
「ある図形やパターンに近いが不完全」というものがあった時は、それを完成させようとする
連続
図形が連続する線の上にあると認識しようとする
（たとえば交差する二つの直線を二つの楔とみなさないこと）
(以上p294 を要約)

これらは尤もな話だし、実験をしてもある程度の範囲でその通りの結果が出るだろう。こういうのがゲシュタルト心理学なわけだけど、間違っている訳でも無いのに下火になった理由は「事が多い」と言っても、そうで無い場合もあるし、とかでそれ以上精密化も定量化も出来なかっただろう。

私の考えではこれはみな確率論によって解決される。つまり、近接の原理は（実際だれもの日常経験通り）同じグループの物が近くにある事が多いという確率的事実の反映だし。見た目が似た物が一つのまとまりを形成している事は多いし、軸の合ったX形の線が見えた時それが二つの>形の物体と<形の物体が合わさった物である事は確率が低い。

そう言えば、ビーダーマンのgeonも同じように確率論で説明できるな。次節の「共通運命の原理」も、たとえばsynchronicityが強いほど、同じ原因から起きたことであるという確率的事実の反映に過ぎない。

脳は、このように生まれてからの全ての知覚/認知の結果を確率的に積算してその確率を使って現在起こっている事に対して推定をかけるようになっている。

▼hack #61 独り言

我々は、ある物体が赤い事、コーヒーカップが２つある事、などは何も考えなくとも把握できる。ただ、「赤い」という情報、「コーヒーカップが２つある」という情報を組み合わせて使うには、言語の構文を利用しなくてはならない、というのがCarruthersの主張だ。(p227)

私としては賛成できない。しかし、私なら言語は単に表現の媒体程度にしか考えないのに対して、西洋人と言うのは言語と知能の本質を結び付けて考える傾向が強い。改めてそう思った。

▼hack #34 注意を向ける能力の限界
p131のattention resolutionが面白い。視覚的注意にはspot lightみたいな物であるあるのは内省的観察でもわかる。このspotlightの空間解像度がこんな簡単な実験で分かるとは面白い。

高次視覚と記憶

Memory modulates color appearance. Hansen T, Olkkonen M, Walter S, Gegenfurtner KR., Nat Neurosci. 2006 Nov;9(11):1367-8.

バナナは全くの白黒写真にしても黄色く感じられるし、ミカンはオレンジ色に、レタスは緑に感じられるという話し。

色知覚が記憶によって強くmodulateされるという話は何となく知っていたが、こんなにはっきりした効果であるとはおもしろい。しかも、「バナナは黄色」みたいな典型の記憶になっていると。

私の考えでは、色知覚だけでなくて、高次視覚は基本的な仕組みとして記憶を利用していると思っている。

知覚が記憶に長時間影響を受けると言う事でマッカロー効果に関係があるかと思ったけど、あれは補色になるらしいので反対だな。

確率モデルを観測された確率分布にfittingする時の、規格化の問題

θiでパラメトライズされたモデル確率分布q(x; θ)を観測された確率p(x)にfittingしたい事が良くある。というか多変数確率論や学習の問題はほとんどそれだ。

fitの基準としてはKL-divergence
　　L = ∫dx p(x) log q(x)/p(x)
がもっとも普通である。ここでxは一般の多次元だと言う事に注意。

確率は全ての場合をあわせれば１なので、当然常に∫dx q(x;θ1,θ2, ...) = 1が前提なんだけど、これが困る。一般的には q(x;θ1,θ2, ...) = c Q(x;θ1,θ2, ...)として、Qが何かの関数であるとして、∫dx q(x) = 1 となるように c を選ぶ。しかしQ(x; θ)が初等関数である場合ですら、積分してcを求める事なんかめったに出来ないから困る。

そこで、どうするか。たとえばHyvarinenの解答については前に紹介した。しかし、これもKL-divergenceじゃない根拠の怪しい基準を使っているのが問題だ。

この問題に対して私は以下のようにすれば良いのではないかと考えた。

まずgradientさえ分かればよい。
　　∂i L = ∫dx p(x) ∂i log q(x)　　　　　　　(1)
ここで∂i は∂/∂θiの略記であり、 ∂i p(x) = 0を使った。またθ0 = cとする。微小変化のベクトルΔθが∫dx q(x) = 1を壊さないための拘束条件は
　　∫dx Σi Δθi ∂i q(x) = 0
であり、ベクトル∂i Lから拘束条件に平行な成分を引き去ればよい。

これは本物の神経のネットワークにとっては、出来ない訳では無いけど、難しい仕事だ。∂i L とベクトル∫dx ∂i q(x)を別々に求めてストアしなくてはならないので２倍多くのニューロンが必要になる。
またp(x) ~ q(x)に収束してくるとベクトルがほとんど同じになって、ニューロン上で桁落ちのような状況が起きるので、かなり精密にキャリブレーションする必要がある。（コンピュータでsimulateする時は、それほどの困難では無いけど。）p(x)は実データの有限サンプリングであるからノイズを含んでいるし、q(x)はMetropolis法（別名Markov Chain Monte Carlo法）によるランダムサンプリングであるためにノイズを含んでいるから状況はなお悪い。

しかし、ここで魔法のように上手く行く公式を見つけた。
　　∫dx (p(x)-q(x)) ∂i log Q(x)　　　　　(2)
を勾配にすればよいのだ。これは上の引き去るべき成分がちょうど１である事を意味する。

導出の仕方をふた通り示す。

第一の導出法はKL-divergence
　　∫dx p(x) log q(x)/p(x)
　　=∫dx p log (1+(q/p-1) )
　　~∫dx p { (q/p-1) - (q/p-1)^2/2 + (q/p-1)^3/3 -...}
で、展開の１次の項がqに拘束条件が課されていない時に悪さをすることが分かる。ならばKL-divergenceから∫dx p (q/p-1)を引けばよい。そうすれば、qの総確率が１からずれた時に復元力も働く。
　　L' =∫dx p log q/p -∫dx (q-p)
で、これgradientは
　　∫dx (p(x)-q(x)) (∂i log Q(x))　　　　(3)
と
　　∫dx (p(x)-q(x)) (∂c log c)　　　　　(4)
である。仮に(4)式により完璧にnormalizationが出来ているとした時、(3)式に従えばcを知る事無しに勾配を求める事ができる。それが(2)式。そして、実のところ、∫dx q(x) (∂i log Q(x))の計算はMetropolis法でやる訳だから、Q(x)が与えられた時のq(x)のサンプリングはcを知る事無しに自動的に行う事ができる。

第二の導出法はp(x)~q(x)のときはKL-divergenceは
　　-∫dx 1/p(x) (q(x)-p(x))^2/2
で近似される。これの勾配は
　　∫dx (∂i q(x))(p(x)-q(x))/p(x)
p(x) ~ q(x)だから
　　〜∫dx (p(x)-q(x)) (∂i q(x))/q(x)
これは(3),(4)式と同じ。

註：この勾配法において、cの値をadjustしてから、シータを更新しまたcの値をadjustし・・・とかやるのは上手く行かない。たとえばc exp(-θx^2)で考えてみよ。c =0 に落ち込んでしまう。