日本神経回路学会誌 Vol. 18, No. 4（2011），194–203

解説

対数線形モデルによるマルチニューロンスパイクデータ解析—定常スパイク相関解析から動的スパイク相関解析への橋渡し—

島 崎

秀 昭

理化学研究所 脳科学総合研究センター∗

Analysis of Multiple Neural Spike Train Data using the Log-linear Model—From Stationary to Time-Varying Spike Correlation—
Hideaki Shimazaki
RIKEN Brain Science Institute∗
概要
神経細胞集団は協調して情報処理を行っています．特に比較的短い時間内に複数の細胞が発火する同期発火現象はセル・アセンブリ仮説やシンファイア・チェーンといった細胞集団による情報処理の仮説との関係から注目を集めてきました．近年実データの同期的発火構造を調べるために対数線形モデルを用いた解析が広く行われています．しかしこれまでの対数線形モデ
ルは定常性を仮定しており行動下の動物から得られるスパイクデータへ適用することはできませんでした．本稿は対数線形モデルの基本事項，パラメータの基本的な推定法（最尤推定・

MAP 推定）を経て，最終的に動的な相関構造の推定を可能にする非線形フィルタの構築に至るまでを紹介します．

1. は じ め にど，この仮定が成り立つと期待できる条件下で記録さ神経細胞は集団で協調して情報を処理していると考えられています．複数の細胞が依存しあって活動を行っ方で行動下の動物の神経細胞の活動はよりダイナミッ
クです．刺激や行動に応じて個々の細胞の発火頻度はているとすると，観測された同時記録スパイク時系列には同期的な集団発火が観測されるかもしれません．

時間的に変動します．同一の試行を繰り返し行い peristimulus time histogram（PSTH）やカーネル密度

Schneidman らの論文1) を皮切りに，同期発火活動の推定法などで個々の神経細胞の発火頻度の刺激依存の統計モデルとして統計物理学のイジングモデル（スピングラスモデル）を援用するデータ解析が広く行われダイナミクスを調べることは電気生理実験の基本的な手法です．さらに発火頻度だけなくスパイク相関構造ています．さらに最近では 2 次のスパイク相関では説もダイナミックに変動する可能性があるとすればどう明のできない，高次のスパイク相関構造にも注目が集
まっています2) ．そこで，本稿では低次から高次まででしょうか．実際 Cross correlogram，Joint PSTH，
Unitary Event 法を用いた解析で皮質や海馬体の 2 つのスパイク相関を統一的に扱うことのできる対数線形の神経細胞の相関が行動に依存して動的に変動するこモデルを用いたマルチニューロンのスパイクデータ解
析を説明します．

とが示されています3〜6) ．こうしたデータに対して通常の対数線形モデルを用いた解析を行うことはできま

多くの対数線形モデルによる解析は定常性，すなわ
ち個々の細胞の発火頻度や相関構造が時間に依存しな

せん．そこで本稿では定常モデルから出発して，こうれたスパイクデータに対して適用されてきました．一いことを仮定しています．そのため，これまでの実デーしたダイナミックなデータを扱うための状態空間モデルによるアプローチ7) を駆け足で紹介します．

タの解析はスライス網膜神経節細胞（Retinal ganglioncell）や皮質培養細胞，麻酔下のサルの皮質神経細胞な2. 高次スパイク相関とはなにか∗ 〒 351–0198

埼玉県和光市広沢 2 番 1 号

議論の始めとして，まず本稿で扱う同時刻スパイク

島崎秀昭：対数線形モデルによるマルチニューロンスパイクデータ解析

195

るなかで，非常にまれに全ての神経細胞が同期発火す
る状況を想像してください．このようなスパイク時系
列では低次のスパイク相関が存在せず高次のスパイク
相関だけが存在する可能性があります．その理由は，
全ニューロンの同期発火という異常事態が部分集団で見た場合には偶然に起こりうる同期発火と見分けがつかないからです．私を含めて発火頻度や 2 次のスパイ
ク相関だけでなく高次のスパイク相関も神経細胞が協図1

秘密の暗号で仲間探し．

相関とはなにかを明らかにしなければいけません．こ調して情報処理を行う様子を明らかにするために重要であると考える研究者がいます．なぜ，高次スパイク
相関に興味があるのでしょうか．このことを説明する
ために，もう一度例を挙げて説明してみましょう．

こではまだ上述のようなダイナミクスは考えないこと

上述のパンチカードの穴を Twitter のようなネットにしましょう．2 次の相互（交互）作用項までを考えた対数線形モデルは多くの文献がありますので，以下上での発信（発言や会話）に置き換えて考えてみましょう．ある適当な時間区間（例えば 10 分程度）の中で一では高次スパイク相関とはなにかを紹介しようと思います．
それでは始めに次のような例から考えて見ましょう．
ある学校のクラスの 5 人の生徒たちに先生が紙切れを配りこう言いました．
「この紙には秘密の暗号で君達の回でも発信をしていれば 1，発信がなければ 0 とする
二値時系列を考えます．二人の人が時間区間内に定期的に会話をしていればその二値時系列は高い 2 次相関を持つことが期待できます．また，各人が勝手に発言

仲間の情報が書かれています．暗号を解読して自分のしていれば独立な二値時系列です．さて，ここで独立に発言していた人々が突然一斉に発言を始めたとしま

仲間を見つけてください．
」紙切れにはパンチカードの．どうやら同ように所々穴が開けられています（図 1）

しょう．この現象をどう解釈するのが適当でしょうか（発言内容は見ることができないとします）．1 つの解じ箇所に穴が開けられていると仲間である確信が持て

釈として，何か特定の事件が起きてそれに各人が独立そうなのですが…．まず生徒は隣同士で紙を見せ合ってみます．すると 1〜2 箇所は穴が一致しています．隣に反応したとみなすことができそうです．すなわち観測できない共通の外部入力による影響が高次相関に現の子は仲間なのでしょうか．紙切れをよく見るとそれぞれの紙には 10 個の升目のうち 4 つに穴が開けられれる可能性があります．

ています．2 つの紙にランダムに穴を開けても同じ箇
所に穴が開いている確率は 0.4 × 0.4 = 0.16 ですから
10 箇所のうち 1〜2 箇所は穴が一致していても不思議
ではなさそうです．それでは次に 3 人が集まって紙を見せ合ってみます．1 箇所ですべての紙に穴が開いています．しかしやはりランダムに穴が開けられていると仮定した場合，(0.4)3 = 0.064 の確率で 3 箇所すべ

どんなに観測技術が発達しても私達が観測できるのは全体のごく一部です．その一部の情報を使って全体の様子を最大限推察したいという希望が高次相関解析の動機の 1 つに挙げられます．こうした高次のスパイ
ク依存関係を簡単に作り出してくれるのが，次に紹介する対数線形モデルです．

3. 対数線形モデルてが一致することがあります．穴が一致した頻度は偶然から考えられるより多いですが，3 人が仲間である

本節では対数線形モデルの定義と良く知られている
性質を簡単にまとめます．本節は多くの優れた文献をと確信は持てません．そこで 5 人が集まってみたとこ
ろ 5 人の全員の紙に同じ箇所で穴が開いていました．

参考に本稿の目的である動的スパイク相関の推定手法の構築に必要となる諸性質のみをまとめたものです．

ランダムに穴を開けた場合その確率は (0.4)5 = 0.0102

より詳しい内容は原著論文8〜10)を，より包括的な文脈ですから，これは故意に同じ場所に穴を開けたと考えても良さそうです．どうやら 5 人全員が仲間だというでの解説は参考文献11〜13)等をご覧下さい．
N 個のニューロンによる同時刻スパイクパターンを結論に達しました．

考えます．ひとつひとつのニューロンの状態は二値の確部分だけを見てもわからない依存関係，それが高次の依存関係です．神経スパイクの話に戻りましょう．

率変数（binary random variable）Xi (i = 1, . . . , N )で表されるとします．ここでは，‘1’ はスパイク，‘0’

複数の神経細胞が低頻度でぱらぱらと独立発火していはスパイクが生成されなかったことを表すとします．

196

日本神経回路学会誌 Vol. 18, No. 4（2011）

N 個のニューロンの集団スパイクパターンはベクトル
X = (X1 , X2 , . . . , XN ) で 表されます．すべての可能なパターンの数は 2N 個です．
N 個のニューロンの同時確率関数（joint probability
mass function）を考えます．ニューロンがとりうる値
として二値変数 xi = 1 or 0 (i = 1, . . . , N ) を使って同時確率関数をP (X1 = x1 , . . . , Xn = xn ) = p(x).

(1)と書くことにしましょう．右辺では簡単のためにベク
トル x = (x1 , x2 , . . . , xN ) を使いました．確率関数と
しての制約



p(x) = 1 がありますから，すべてのス

パイクパターンの確率は 2N − 1 個のパラメータを決め
れば定まります．当然 2N − 1 個のスパイクパターンのそれぞれについて，直接そのパターンが生じる確率

p(x) を決めれば同時確率関数は定まります．しかし，それ以外にも同時確率関数を定める便利なパラメータの組があります．情報幾何の分野では，2N − 1 個のパラメータをすべての分布 {p(x)} をカバーする 2N − 1
次元の ‘座標’ と見なします．座標系の一点・一点が特定の確率分布を表しています．座標系はある確率分布
から別の確率分布へ移動するときの道先案内人の役を担います11) ．ここでよく登場するのが θ-座標系と η座標系です．
3.1 θ-座標系と η-座標系
同時確率関数の対数を次のような展開式で表します．

log p(x) =i

θi xi +



します．これを同時スパイクレートと呼ぶことにします．どのニューロン同士の同期発火の期待値なのかは添え字で表現されています．以下ではこれらをベクト
ル η = (η1 , η2 , . . . , η12 , η13 , . . . , η1···N ) を用いて表します．
式 (2) 及び式 (3) は次のような ‘マルチインデック
ス’ 表記でコンパクトに書くことができます．Ωk を

1 から N までの自然数から k 個を取り出してきた
時のすべての組としましょう：Ω1 = {1, 2, . . . , N },
Ω2 = {12, 13, . . .}, Ω3 = {123, 124, . . .}, etc. この要
素 I ∈ {Ω1 , . . . , ΩN } を使って，自然パラメータ・期
待値パラメータを θI ，ηI と表すことができます．同様に式 (2) の相互（交互）作用項 を
i = 1, · · · , N
i<jFi (x) = xi ,Fij (x) = xi xj ,..
.
F1···N (x) = x1 · · · xN .

(4)とすれば，FI (x) (I ∈ {Ω1 , . . . , ΩN }) を使って式 (2)・
(3) は簡潔に

log p(x) =θI FI (x) − ψ(θ),(5)

I

ηI = E[FI (X)].

(6)と書くことができます．上記の式からわかるように式(5) の対数線形モデルは連結関数（ここでは多変量ロ

θij xi xj + · · ·ジット関数）をとった後にパラメータに対して線形和でi<jあれば良く，特徴関数と呼ばれる FI (x) の部分はデー· · · + θ1···N x1 · · · xN − ψ(θ).

(2)

タを非線形に加工しても良いわけです．そうすると式(θ1 , θ2 , . . . , θ12 , θ13 , . . . , θ1···N ) ．本=
− log p({0, . . . , 0}) はp(x) = 1 を満たすよう
xに決められる（対数）正規化項です．そのため ψ(θ)
は θ に関する複雑な式になっています．式 (2) が対数線形モデル（log-linear model）です．パラメータ
θi , θij , . . . , θ1···N は指数分布族の自然パラメータある
いは正準パラメータと呼ばれ，‘θ-座標系’ を構成します．
もうひとつ重要な座標系が η-座標系です．η-座標系(4) による展開以外にもっと良いデータの加工の仕方

ここで θ



=

稿 で は プ ラ イ ム 記 号 で 転 置 を 表 し ま す ．ψ(θ)

ろですが，今回は式 (4) による展開を用いた相関解析に焦点を当てて先を急ぎましょう．

3.2 対数線形モデルのキュミュラント
先ほど分布の正規化に用いた関数 ψ(θ) は θ に関す
る複雑な関数だと書きました．実際に ψ(θ) を書くと次のようになります．

ψ(θ) = logは次のように構成されます：

ηi = E [Xi ] ,ηij = E [Xi Xj ] ,..
.
η1···N = E [X1 · · · XN ] .

があっても良いように思えます．とても興味深いとこ  θ F (x)
e

I

I

I


.

(7)

x

i = 1, · · · , N
i<j

一般に指数分布族のキュミュラント母関数 K(u) は(3)

これらのパラメータは同時スパイク生成の期待値を表

K(u) = ψ(θ +u)−ψ(θ) で与えられ，確率変数 FI (X)の高次のキュミュラントは ψ(θ) の高階導関数で表されます．θI による変微分を ∂I と書くことにします．
対数正規化項 ψ(θ) を θI に関して一階微分すると

島崎秀昭：対数線形モデルによるマルチニューロンスパイクデータ解析

197∂I ψ(θ) = E[FI (X)]= ηI .

(8)

次に ψ(θ) の二階導関数は二次のキュミュラントを与えます．

∂I ∂J ψ(θ) = E[(FI (X) − ηI )(FJ (X) − ηJ )]= ηI∪J − ηI ηJ .

(9)

ここで記号 ∪ は本稿だけで使う特別な表記です．マルチインデックスの要素を集合とし，和集合を再びマ
ルチインデックス表記にしたものとします．たとえばI = 12 及び J = 234 のとき I ∪ J = 1234 とします．
∂I ∂J ψ(θ) はフィッシャー計量で ∂I ∂J ψ(θ) = gIJ と書くことにします．これから様々なシチュエーション
で対数線形モデルの 2 次近似式を使います．そこで上記の式 (8)・(9) はとても重要な役割を担います．

3.3 部分モデルと高次スパイク相関
対数線形モデルの高次相互作項の θ-パラメータは高次のスパイク相関を表します．ただし θ-座標系によるスパイク相関は対応する η-座標系との関係において解

図2

高次スパイク相関の幾何的な説明．

釈されます．具体的には非ゼロの r 次以上の相互作用
項 {θI }（I ∈ {Ωr , . . . , ΩN }，r ≤ N ）は r 次以上の同期スパイクレート {ηI }（I ∈ {Ωr , . . . , ΩN }）が r 次未

用います）
．するとこの作業は p̃ (x) のエントロピーを満の同期スパイクレート {ηI }（i ∈ {Ω1 , . . . , Ωr−1 }）最大にすることに等しく，p (x) に最も近い Mr 上のモから偶発的に起こるとして期待される値（チャンスレベル）より多いかあるいは少なくなることを表していデルは最大エントロピーモデル p∗ (x) であることが示せます．次に下図の Er の中で q (x) に最も近いモデ

ます．甘利ら8〜10) は θ-座標系と η-座標系の双対性を巧みに利用した情報幾何の枠組みによりこれを簡潔にル p∗∗ (x) を探します．すると p∗∗ (x) は最大エントロ
ピーモデル p∗ (x) に一致することが（θ-座標系と η-座

説明しました．ここでは情報幾何の詳細には踏み込ま標系の双対性に基づくピタゴラスの定理によって）示ず，甘利9) に従ってそのエッセンスだけを紹介します．
N 個のニューロンのスパイクパターンの分布 q (x)されます．
q (x) をデータ分布（経験分布）とし，これに最も近を考え，その分布の r 次以下の同期スパイクレートをい Er 上の分布を探す作業を最尤法といいます．従っ{ηi }（i ∈ {Ω1 , . . . , Ωr }）とします．次に r 次以下の同期スパイクレートが q (x) と同じモデル p̃ (x) を考えます．この制約の下で p̃ (x) は r 次より大きい同期ス
パイクレート {η̃I } （i ∈ {Ωr+1 , . . . , ΩN }）を自由にとることができ，これによって p̃ (x) のスパイクの高次相関構造が変わります．図 2 上部ではこのモデル群を
Mr と表しています．この中で p̃ (x) のエントロピーを最大にするモデルがあるはずです．これを p∗ (x) として r 次以下の同期スパイクレートをデータ分布に一致させるという制約の下で最大エントロピーモデルを求ておきます．一方，対数線形モデルにおいて使用する相互作用項を 1 次から r 次までのみに制限した部分モ

めることは，データに対して r 次までの対数線形モデ
ルを最尤法によってフィッティングすることと同じです．多くの文献では最尤推定を行っており，そのため
対数線形モデルは最大エントロピーモデルと呼ばれることも多いです．r 次の最大エントロピーモデルでは

r + 1 次以上の同期スパイクレートに対して制約が設けてられていないため，これらは（r 次以下の同時ス
パイクレートの制約の下で）偶発的に生じます．逆にデルを考え Er と書くことにします（図 2 下）
． Er に非ゼロの r + 1 次以上の相互作用項 {θI } はこの偶然含まれるモデルをここでは p (x) とします．再び上図に戻り Mr 上で p (x) に最も近い分布を探します（こに起こる頻度からのずれを表します．これが対数線形モデルの相互作用項によって定義される高次スパイク

こで ‘近い’・‘遠い’ という概念にはカルバック距離を相関の意味になります．

198

日本神経回路学会誌 Vol. 18, No. 4（2011）1   tFI X .
nn

4. 定常モデルのデータへの適用

kI =4.1 最尤推定

(11)

t=1

次に r 次の定常対数線形モデルのデータへのフィッ
ティングに話を移しましょう．まず考えられるのが上で挙げた最尤法によるフィッティングです．
同時記録された複数のニューロンのスパイク時系列
を幅 Δ[s] のビンに区切ります．ビンの数は n 個あったとします（観測時間が nΔ）．t 番目のビンの観測をt
二値データで表し Xt = (X1t , X2t , . . . , XN
) とします．

ここで ‘1’ の値を含むビンはその中にひとつあるいはそれ以上のスパイクが存在したことを示します．またつまり kI は同時スパイクレート（単位は Δ）の観測値
であり，ηI の推定量です．式 (10) を見ると r 次までの相互作用項を入れたモデルの推定には r 次の同時ス
パイクレート k（k = [k1 , . . . , k12 , . . . , k1···r , . . .] ）があれば十分であることがわかります．例えば 2 次までの相互作用を考えるモデルの最尤推定には，発火レー
トと二つのニューロン間の同時スパイクレート（ただしすべての組について）があれば十分で，それ以上の‘0’ の値はそのビンの中にスパイクがひとつもなかった

高次の同時スパイクレートはパラメータの推定値に影響を与えません．このような性質を持つデータの縮約ことを示します．対数線形モデルの式にデータを投入
してパラメータの関数と見る尤度関数を考えます．

量を十分統計量といいます．
最尤推定量を θ MLE とすると形式的に



p X1 , . . . , Xn |θ
=

n



p Xt |θ

t=1
n

=



exp[

(12)

θ



で与えられますが，実際に最適解を求めるには次のようにします．目的関数を

I∈Ω1:rI∈Ω1:r= exp[n(θ MLE = arg max log p (k|θ) .

t=1= exp[n(I∈Ω1:rt

θI FI Xn1f (θ) = − log p (k|θ)− ψ(θ)]

t= −n k θ − ψ (θ) .

(13)

FI X

− ψ(θ t ))]として f (θ) を最小にする θ を求めることを考えます．

kI θI − ψ(θ))].

(10)f (θ) は θ に関して凸関数になることが保証されています．そこで f (θ) を簡単な二次近似の式に直してか

θI

n

t=1

始めの等式ではデータの定常性を仮定して，かつビン毎に n 回独立なサンプリングを行ったという仮定を用いています．実際にビン毎にサンプルが独立かというとかなり怪しい面があることは否めません．神経細胞には不応期が存在し，ある時刻におけるスパイクの発生率は過去のスパイク，なかでも直前のスパイクからの経過時間に大きく依存すると考えられています．また複数神経細胞を考えた場合，直接的なシナプス結合
及びネットワークを介した間接的な影響で過去の他の神経細胞のスパイクの活動履歴がある神経細胞のスパ
イク生成率に影響を与えることが考えられます．実際，

ら近似式の最小値を求めることを繰り返し行います．









θ new = min{f θ old + f θ old , θ − θ oldθ

1+θ − θ old }.
2

(14)

近似式の最小値はその一回微分が 0 になるところです．



f θ old +1θ − θ old = 0.
(15)

従って更新式θ new = θ old −   f θ old



(16)

スパイク履歴の影響を導入しその代わりに同時刻の条
件付独立を仮定するモデルが発展しています14) ．一方が得られます．この式は最急降下法（Gradient descent）の公式です．もう少し 2 次の項の近似を良くす

で本稿のように複数ニューロンの同時刻での同時分布

るように二階微分（ヘシアン）を入れます．

を考えたい場合，加えて観測履歴もうまく導入するためには少し腰をすえて考える必要がありそうです．二兎追うものは一兎も得ず… とは言いませんが，今回はビン毎の独立性を仮定して話を進めましょう．
最後の等式で出てくる kI（I ∈ {Ω1 , . . . , Ωr }）は次のように定義しました．









θ new = min{f θ old + f θ old , θ − θ old
θ

 

−1 1 +
θ −θ old f θ oldθ −θ old }.
2
(17)
この式から自然勾配法15)と呼ばれる式が得られます．

199

島崎秀昭：対数線形モデルによるマルチニューロンスパイクデータ解析

θ new = θ old −  f θ old

−1





 f θ old . (18)

著に現れます．なにが起きたのでしょうか？ 理由は「データにスパイク同期発火がなかったから」です．例

 = 1 のときは Newton-Raphson 法です．
最尤法の場合の具体的な更新式を求めるために f (θ)

えば 2 つのニューロンで考えた場合，この 2 ニューロンが観測期間中に一度も同時に発火しなかったときの勾配とヘシアンを求めます．式 (8)・(9) は部分モデルでも成り立つことに注意して，f (θ) の勾配は

θ12 → −∞ となります．つまり最尤法では，2 つのf (θ) = −n {k − ψ (θ)}

ニューロンは非常に強い負の相関関係にあって両者が同時に発火するのを拒んでいるという解釈になります．

(19)

これは与えられたデータに対しては邪念のない解釈でここで式 (8)（ψ (θ) = η ）を用いました．次に f (θ)

す．しかし私達が扱うスパイクデータの発火率は低く，たとえニューロンが独立（データが 0 に偏っている）= −n (k−η) .
のヘシアンはに発火していても同期発火確率は非常に小さくなります．従って短いデータ区間では同期発火がたまたま見f (θ) = n  ψ (θ)
≡ nG.

(20)

ここで式 (9) の要素を持つフィッシャー情報行列 G を定義しました．
従って最尤法の場合の更新式は最急降下法・自然勾
配法それぞれ






−1


old

θ new = θ old + n k−η old .
θ new = θ old + G

k−η

(21)
.

(22)

となります．自然勾配法ではフィッシャー情報行列もη old を用いて計算します．これらの最適化公式は誤差
項 k−η を最小化するように θ を更新します．極値である f (θ) = 0 の解が η = k であり，最尤法ではモデルの同時スパイクレート ηI が実測値 kI と一致す
るように選ばれます．
さて，どちらの式も θ old に対応する η old を計算する必要があります．周辺量である η を直接得るためには対数正規化項 ψ (θ) の計算が必要になります．10 個程度の少数ニューロンならばすべての確率関数を計算してしまうことも十分できます．しかしニューロン数が増えた場合，実はこの変換は容易ではありません．
対数線形モデルの難しさはほとんどすべてこの座標変換が容易でないことに起因するといっても過言ではありません．通常はモンテカルロ法を使って正規化項を計算することなく η を求めることが行われています．
ここでは η を求めることができるとして先に進みましょう．

られなかった可能性もあり，最尤推定量によるモデルは新しいデータを説明するのには適さないかもしれません．この問題は良く知られたモデルの過剰適合（オーバーフィッティング）の問題です．この問題を回避するためには，適切な次数の対数線形モデルを選ぶための交差検定や赤池情報量基準（AIC）によるモデル選
択，あるいはパラメータの正則化を行う必要があります．ここでは先を急ぎ事前分布を用いたベイズの定理に基づく正則化を考えます．

4.2 MAP 推定
ベイズの定理を用いる推定法ではパラメータ θ に分
布（事前分布）を想定します．ここでは事前分布として
平均 0，共分散行列 λ−1 I の正規分布を考えましょう．

p(θ|λ) =

λ
2π

d/2

1
exp(− λ θ).
2

データが与えられた下でのパラメータ θ の密度分布（事後確率密度関数）はp (θ|k, λ) =p (k|θ) p (θ|λ)
p (k|λ)(24)で与えられます．この事後分布の最頻値（Maximum

a posteriori，MAP 推定量）
θ MAP = arg max p (θ|k, λ)

(25)

θ

を求めたいのですが，これを探すのに事後分布の分母の計算は必要ありません．そこで目的関数を最尤法では最急降下法を用いるのが一般的です（そf (θ) = − log p (k|θ) − log p (θ|λ)の理由は後で説明します）．それではこの方法を使ってデータに対数線形モデルをフィッティングしてみま= −n k θ − ψ (θ) −しょう．するとどうでしょう，いくつかの θI について

(23)

d
λ
λ
log
+ θ
2
2π
2
(26)

その最尤推定量が負の非常に大きな値になってしまいます！ この現象はデータサイズ n が小さいとき，加として f (θ) を最小化する θ を求めます．最小化の手
順は最尤法のときと同じで，式 (16) もしくは式 (18)

えて高次スパイク相関をモデルに組み込んだときに顕に対して上記の f (θ) を用います．そこで f (θ) の勾200

日本神経回路学会誌 Vol. 18, No. 4（2011）

配及びヘシアンを求めておきます．

5. 動的スパイク相関の推定

 f (θ) = −n (k−η) + λθ,(27) f (θ) = nG + λI.

(28)

最尤法における式 (19)・(20) と比べると式 (27)・(28) には正則化項が第二項として加えられています．式 (27)を見るとこの項によって過剰適合が避けられることがわかります．一方 MAP 推定量によるモデルの同時ス
パイクレートは実測値からずれることがわかります．
また前節の最尤法では，モデルに含まれていながら観測されない同期スパイクがあるときに自然勾配法（式(18)）による更新はうまくいきません．フィッシャー
情報量行列の対応する行（列）の値がすべて 0 に近づき逆行列を求められなくなるからです．事前分布を用いた正則化を行うと式 (28) 第二項の存在によりこの問題点も回避することができます．
事後分布を用いた推定では分布の広がりをみること
で信用区間（信頼区間のベイズ版）を得ることができます．特に事後分布を正規分布で近似するとその平均
は MAP 推定量 θ MAP ，共分散行列 W は

W = −  log p (θ|k, λ)|θ=θ MAP

= f (θ)|θ=θ MAP

−1

−1

= [nG + λI]−1

ら記録されたスパイク時系列データの発火頻度や同期
発火頻度は変動します．本節ではこうした実データに対応できるよう状態空間モデルの枠組みを用いて対数線形モデルを時間依存の形に拡張します7) ．準備が長くなりましたがここまでくれば動的スパイク相関の推定に必要な手段は揃っています．
5.1 動的スパイク相関の状態空間モデル
状態空間モデルでは観測されない状態変数を考え，隠れ状態の時間的遷移を仮定します．ここでは対数線形モデルの相互作用パラメータ θ を隠れ状態とし，θ が時刻（ビン）t (= 1, . . . , T ) に依存するとします．時刻 t における相互作用パラメータを θ t =tt[θ1t , . . . , θ12, . . . , θ1···r, . . .] とし，変動の様子を記述す
る状態方程式を仮定します：

θ t = θ t−1 + ξ t ,(31)(t = 2, . . . , T )．ここで ξ t は平均 0・共分散 Q の多変量正規分布に従う乱数とします．初期値は θ1 ∼ N (μ, Σ)とします．これは観測期間にわたる隠れ状態 θ 1:T ={θ 1 , θ 2 , . . . , θ T } の同時分布として，事前分布(29)で与えられます．この近似はラプラス近似といい後述の動的スパイク相関の推定においても使用します．
尤度と事前分布のどちらを重視するかは超パラメー
タ λ によって決まります．通常 λ は周辺尤度のパラメータとしてここまで同時確率関数 p (x) は時間によらず一定であることを仮定してきました．しかし行動中の動物か

p (θ 1:T ) = p (θ 1 |μ, Σ)Tp (θ t |θ t−1 )

(32)

t=2

を想定していることになります．ここで p (θ t |θ t−1 )は状態遷移を表す確率密度で式 (31) から正規分布

N (θ t−1 , Q) になります．共分散行列 Q は相互作用(30)パラメータが時間的にどの程度速く変動するのかを決める平滑化パラメータの役割を果たします．事前分布を最大にするものが選ばれます（経験ベイズ法・TypeIIに含まれる超パラメータをまとめて w = [μ, Σ,Q] とします．

の最尤推定）．周辺尤度はラプラス近似を用いて近似式

状態空間モデルではさらに時刻（ビン）毎の条件付を得ることができます．次数の異なるモデルの比較も
周辺尤度比（ベイズ因子）を用いて行うことができま

独立を仮定します．すなわち時刻 t の θ t が決まった
らこの同時確率関数に従ってデータが生成されるとしす．上記の事前分布ではひとつの超パラメータを考えましたが，パラメータ毎に異なる複数の超パラメータと同じようにp (k|λ) =p (k|θ) p (θ|λ) dθ

を考えることもできます．このような方法は関連ベク
トルマシンと呼ばれスパース推定の文脈で研究されています．この場合の異なるモデル間の比較では超パラメータ数がパラメータ次元に依存するため，周辺尤度を用いた赤池ベイズ情報量基準（ABIC，単に AIC と呼ぶことも多いようです）や，さらに上位の階層ベイズによる周辺化を考える必要があります．

ます．従って時刻 t の条件付対数線形モデルは式 (10)p (x | θ t ) = exp[I∈Ω1:rθIt FI (x) − ψ(θ t )] (33)となります．このときモデルの同時発火レートは

ηIt = E [FI (x) |θ t ] ,(34)(I ∈ {Ω1 , . . . , Ωr })．これをベクトル η t で表します，tt
η t = [η1t , . . . , η12, . . . , η1···r, . . .] ．

島崎秀昭：対数線形モデルによるマルチニューロンスパイクデータ解析

電気生理実験では同じ試行（トライアル）を繰り返
し行い，その加算平均を取ることがしばしば行われます．そこでマルチニューロンの同時記録スパイク時系
列が n 回の同一試行により得られたとします．このスパイク時系列を適当な幅 Δ [s] のビンに区切りビン
中に 1 つ以上のスパイクがあれば ‘1’ とし，なければ ‘0’ とします．ビンの数を T 個とします．この二値データを Xt,l = (X1 , X2 , . . . , XN ) であらわしま
t,l

t,l

t,l

t,l

す．要素 Xi

は i(= 1, . . . , N ) 番目のニューロンのl(= 1, . . . , n) 試行目の t(= 1, . . . , T ) 番目のビンのス
パイクの有無を表します．ここで定常モデルのときとサンプル数 n の定義が変わっていることに注意してください．定常性を仮定できるのは同一試行を繰り返し行ったときの同一時刻のサンプルであるとするからです．そこで時刻 t の観測同期スパイクレートを

1   t,l 
kIt =
FI X
n
n

(35)

l=1

(I ∈ {Ω1 , . . . , Ωr }) とします．これを用いて全観測
データ k1:T = {k1 , k2 , . . . , kT } によるパラメータθ 1:T = {θ 1 , θ 2 , . . . , θ T } の尤度関数は

=

exp[



=


t=1



exp[n(

I∈Ω1:rl=1 t=1TI∈Ω1:rt,lθIt FI X− ψ(θ t )]p (k1:T , θ 1:T |w) dθ 1:T(38)を最大化するものを選びます．長くなるので割愛しま
すがこの式もまたラプラス近似により近似式を求める
ことができます．また異なる次数の相関を組み込んだモデルの適合度の比較にも周辺尤度を用いることができます．ただし超パラメータがパラメータ数に依存する場合，超パラメータの次元のペナルティーを入れた
赤池ベイズ情報量基準等でモデル選択をすることにな
ります．
実際の運用では式 (37) のような大掛かりな事後分
布の構築と超パラメータ最適化を同時に行うことはせず，EM アルゴリズムを用いて交互に行うことが普通です16〜18) ．この方法では超パラメータ w を決めて事後分布を構築し（E-step），この事後分布を用いて w を最適化する（M-step）ことを交互に行います．E-step
における事後分布の構築は時間前方へ再帰的に事後分布を構築するフィルタリングを行い，その後時間後方に向けて事後分布を構築する平滑化（スムージング）を行うことで完了します19) ．本稿では最も重要なフィ
ルタリングとスムージングについて説明します．

θIt kIt − ψ(θ t ))]exp[n(kt θ t − ψ(θ t ))]まず初めに一期先予測分布 p (θ t |k1:t−1 ) を計算します．これは時刻 t − 1 までの観測が与えられたもとで事前分布に基づいて時刻 t の状態を予測する分布です．
一期先予測分布は Chapman-Kolmogorov の式として
次のように書くことができます．



T

=


p (k1:T |w) =5.2 フィルタ分布・平滑化分布

p (k1:T |θ 1:T )
n
T



201

(36)

p (θ t |k1:t−1 ) =p (θ t |θ t−1 ) p (θ t−1 |k1:t−1 ) dθ t−1 .
(39)

t=1

と書くことができます．
本稿での目標はデータが与えられた下で相互作用パラメータの事後分布を求めることです．尤度関数（式ここで右辺の遷移確率 p (θ t |θ t−1 ) は事前分布である
状態モデル（式 (31)）で与えられる平均 θ t−1 ・共分(36)）と事前分布（式 (31)・(32)）を与えたので原理

散 Q の正規分布です．また右辺 2 つめの密度関数は
時刻 t − 1 におけるフィルタ密度関数 p (θ t−1 |k1:t−1 )的にはベイズの定理から事後分布はです．ここで時刻 t − 1 におけるフィルタ分布を仮に

p (θ 1:T |k1:T , w) =p (k1:T |θ 1:T ) p (θ 1:T |w)(37)
p (k1:T |w)

で与えられます．事後分布がわかればその MAP 推定量を得ることでデータに照らして最も適切なスパイク
相関のダイナミクスの経路がわかります．また分布を求めることで推定の確からしさを表す信用区間を得る

平均 θ t−1|t−1・共分散 Wt−1|t−1 の正規分布とすると一期先予測分布も正規分布となり，その平均 θ t|t−1 及び共分散行列 Wt|t−1 は次式で与えられます．

θ t|t−1 = θ t−1|t−1 ,(40)

Wt|t−1 = Wt−1|t−1 + Q.

(41)

こともできます．ただし事後分布は超パラメータ w の選択に依存します．超パラメータは定常モデルのとき次にフィルタ分布 p (θ t |k1:t ) を求めます．フィルタ分
布は時刻 t までのスパイク時系列を観測した下での，と同じように周辺尤度

時刻 t における θ t の事後分布です．フィルタ分布は

202

日本神経回路学会誌 Vol. 18, No. 4（2011）

ベイズの定理を使って尤度と一期先予測分布を用いて次のように書くことができます．

p (θ t |k1:t )=p (kt |θ t ) p (θ t |k1:t−1 )p (kt |k1:t−1 ) ∝ exp n kt θ t − ψ (θ t )− −1 1
θ t − θ t|t−1 Wt|t−1
θ t − θ t|t−1 .
2(42)

図 3 フィルタ分布構築の模式図．

ここで MAP 推定量を求めるため定常モデルのときとして第二項のフィッシャー情報行列は時刻 t の観測に同じように目的関数を

f (θ) = − log p (kt |θ) − log p (θ|k1:t−1 )(43)

よって得られる推定精度の向上を示します．以上を模式的に示すと図 3 のようになります．一般に時刻 t の

θ の一期先予測分布の裾野は時刻 t − 1 におけるフィルタ分布に比べて広がります（式 (41) 参照）．時刻 tとおきます．目的関数の勾配およびヘシアンは−1∇f (θ) = −n (kt − η) + Wt|t−1
θ − θ t|t−1−1∇∇f (θ) = nG + Wt|t−1においてデータを観測することにより一期先予測によ(44)

る分布の期待値は修正され（式 (46)），新たな事後分(45)

布は分散も小さくなります（式 (47)）．
最後に平滑化事後分布16, 19)で得られます．式 (16) もしくは (18) に対して上述の導関数を適用し MAP 推定量を得ることができます．

p (θ t |k1:T ) =p (θ t |k1:t )

さらにここでは事後分布を正規分布によって近似します（ラプラス近似）．これには MAP 推定量をフィルタ分布の期待値
  −1 θ t|t とし，共分散行列を Wt|t =∇∇f θ t|tとします．フィルタ分布を正規分布で近似することで一期先予測分布の導出で行った正規分布の仮定と矛盾が生じません．このようにして一期先予測の式 (40)・(41) とあわせて t = 1, . . . , T について再帰的にフィルタ分布の近似分布を求めて行きます．
上述のアルゴリズムは次のような非線形再帰フィル
タ公式を解いていることになります7)．

θ t|t = θ t|t−1 + nWt|t−1 (kt − η t|t ),

(46)

−1
−1
Wt|t
= Wt|t−1+ nGt|t .

(47)

ここで式 (46) は ∇f (θ) = 0 から得ました．式 (46)で θ t|t 及び η t|t は同じ分布の異なる座標系による表記

p (θ t+1 |k1:T ) p (θ t+1 |θ t )dθ t+1p (θ t+1 |k1:t )(48)の構築を後方への再帰的アルゴリズムを用いて行います．これにはフィルタリングにより得られた θ T |T 及
び WT |T から始めて平滑化分布の平均と分散をθ t|T = θ t|t + At θ t+1|T − θ t+1|t ,

Wt|T = Wt|t + At Wt+1|T − Wt+1|t



(49)
At
(50)

−1
（ただし At = Wt|t F Wt+1|t）に従って t = T −1, T −2, . . . 1 の順に求めて行きます．これにより全観測デー
タを用いて時刻 t のスパイク相関状態を推定することができます．

ですからこの式は非線形方程式であり，上述のように解かなければなりません．しかしながら式 (46)・(47)

6. お わ り にによる表記はこの非線形近似フィルタがなにを行って

本稿は対数線形モデルを用いたマルチニューロンス

いるかを端的に表しています（図 3 参照）．まず式 (46)の第一項は時刻 t のデータが与えられる前の事前分布

パイクデータ解析を紹介しました．このモデルによって明らかにできること，特に高次スパイク相関についに基づいた期待値の一期先の予測です．これに対し第て紹介しました．次に定常モデルのパラメータの具体

二項は時刻 t の観測が行われたことによる修正を表します．同様に式 (47) の第一項は事前分布によって得ら的な推定法（最尤推定・MAP 推定）を経て，動的なスパイク相関構造の推定を可能にする状態空間モデル

れる一期先予測の推定の正確さを表します．これに対による定式化までを順を追って紹介しました．それで島崎秀昭：対数線形モデルによるマルチニューロンスパイクデータ解析

は高次のスパイク相関は本当に動物の行動とともに動的に生成されるのでしょうか．それだけでなく動物の高次の認知機能に積極的に関与しているのでしょうか．
その答えは．．．近々日の目をみるであろう私達の論文にぜひ目を通してみて下さい．
謝辞

対数線形モデルを教えて頂いた甘利俊一先生・

中原裕之先生に感謝いたします．また豊泉太郎氏との議論により理解が深まりました．第 2 章で用いたパンチカードの例のヒントを与えてくれた Elad Schneidman
氏にも感謝致します．
参 考

文 献

1) Schneidman, E., Berry, M.J., Segev, R., andBialek, W. (2006): Weak pairwise correlations
imply strongly correlated network states in aneural population, Nature, Vol.440, No.7087,
pp.1007–1012
2) Santos, G.S., Gireesh, E.D., Plenz, D., andNakahara, H. (2010): Hierarchical interaction
structure of neural activities in cortical slice
cultures, J. Neurosci., Vol.30, No.26, p.8720
3) Riehle, A., Grün, S., Diesmann, M., andAertsen, A. (1997): Spike synchronization
and rate modulation diﬀerentially involved
in motor cortical function, Science, Vol.278,
No.5345, pp.1950–1953
4) Sakurai, Y. and Takahashi, S. (2006): Dynamic synchrony of ﬁring in the monkey prefrontal cortex during working-memory tasks,
J. Neurosci., Vol.26, Nol.40, pp.10141–10153
5) Fujisawa, S., Amarasingham, A., Harrison,
M.T., and Buzsáki, G. (2008): Behaviordependent short-term assembly dynamics inthe medial prefrontal cortex, Nat. Neurosci.,
Vol.11, No.7, pp.823–833
6) Ito, H., Maldonado, P.E., and Gray, C.M.
(2010): Dynamics of stimulus-evoked spiketiming correlations in the cat lateral geniculate nucleus, J. Neurophysiol., Vol.104, No.6,
pp.3276–3292
7) Shimazaki, H., Amari, S.-I., Brown, E.N.,
and Grün, S. (2009): State-space analysis on time-varying correlations in parallel
spike sequences, In Proc. IEEE ICASSP2009,203

pp.3501–3504
8) Amari, S.-I. and Nagaoka, H. (2000): Methods of Information Geometry (Translations ofMathematical Monographs). American Mathematical Society9) Amari, S.-I. (2001) Information geometry
on hierarchy of probability distributions, Information Theory, IEEE Transactions on,
Vol.47, No.5, pp.1701–1711
10) Nakahara, H. and Amari, S.-I. (2002):
Information-geometric measure for neuralspikes, Neural Comput., Vol.14, No.10,
pp.2269–2316
11) 赤穂昭太郎（2005）：情報幾何と機械学習, 計測と制御, Vol.44, No.5, pp.299–306
12) 岡田真人（2009）：脳の計算論（シリーズ脳科学1）, chapter スパイクの確率論, pp.223–251. 東京大学出版会

13) 中原裕之（2009）：脳の計算論（シリーズ脳科学 1）,
chapter 意思決定とその学習理論, pp.159–221.
東京大学出版会

14) Truccolo, W., Eden, U.T., Fellows, M.R.,
Donoghue, J.P., and Brown, E.N. (2005): A
point process framework for relating neuralspiking activity to spiking history, neural ensemble, and extrinsic covariate eﬀects, J. Neurophysiol., Vol.93, No.2, pp.1074–1089
15) Amari, S.-I. (1998): Natural gradient works
eﬃciently in learning, Neural Comput.,
Vol.10, pp.251–276
16) Kitagawa, G. and Gersch, W. (1996):
Smoothness Priors Analysis of Time Series(Lecture Notes in Statistics). Springer
17) Brown, E.N., Frank, L.M., Tang, D., Quirk,
M.C., and Wilson, M.A. (1998): A statistical paradigm for neural spike train decoding
applied to position prediction from ensemble
ﬁring patterns of rat hippocampal place cells,
J. Neurosci., Vol.18, No.18, pp.7411–7425
18) Smith, A.C. and Brown, E.N. (2003): Estimating a state-space model from point process
observations, Neural Comput., Vol.15, No.5,
pp.965–991
19) Koyama, S. and Paninski, L. (2010): Eﬃcientcomputation of the maximum a posteriori
path and parameter estimation in integrateand-ﬁre and more general state-space models,J. Comput. Neurosci., Vol.89, No.105, pp.1–2