統計数理（2012）
第 60 巻 第 1 号 173–188
c 2012 統計数理研究所
［原著論文］

状態空間法による神経スパイク発火率の推定
小山 慎介

†

（受付 2011 年 8 月 2 日；改訂 12 月 7 日；採択 12 月 14 日）

要旨

単一神経スパイク時系列から発火率を推定する問題を扱う．この問題を状態空間モデルで表現し，発火率推定を状態の平滑化として定式化する．平滑化のための効率的な近似アルゴリズ
ムを構成し，EM アルゴリズムによるモデルパラメータ推定，及び周辺尤度に基づくモデル選
択の方法を提案する．推定するときに仮定するモデル，および積分発火モデルを用いた数値実験により，本稿で提案する非ポアソン性を仮定した手法の有効性を示す．
キーワード： 神経スパイク時系列，発火率推定，状態空間モデル．
1.

はじめに

神経科学において，記録された神経スパイク列から発火率を推定することは，神経系の機能を調べるための第一歩である．外界からの刺激や行動とスパイク発火率との相関を見いだす研究は，Adrian らの先駆的な試みに始まり（Adrian and Zotterman, 1926），現在でも神経コーディ
ングという名で神経科学の一大分野となっている（Rieke et al., 1997; Dayan and Abbott, 2001）．
スパイク発火率とは，単位時間あたりに発生するスパイク数の期待値で定義される点過程の強度関数で与えられる．発火率 λ(t) は各時刻 t において定義され，微小時間 [t, t + dt) にひとつのスパイクが見出される確率は λ(t)dt + o(dt) である．与えられたスパイク発火時刻の
系列 {t0 , t1 , . . . , tn } から λ(t) を推定する問題を考えよう．最も素朴な方法は，スパイク間隔
を yi = ti − ti−1 (i = 1, 2, . . . , n) として，その逆数を時間順に並べたものである（Bessou et al.,1968）．図 1（a）からわかるように，この方法だと推定した発火率はノイジーである．大脳皮質神経細胞のスパイク発火時刻には規則性がなく，また試行毎の再現性もないことが知られて
おり（Softky and Koch, 1993），そのような確率的とみなせる単一試行のスパイク列データからよりよい発火率推定を得るためには平滑化が必要であることがよく知られていた（Kass et al.,
2003; Shinomoto, 2010）．またこのことから，発火率を推定する際にスパイク発火のポアソン性を仮定することが多い．例えば Shimazaki らはヒストグラム法およびカーネル法（Shimazaki，Koyama らは状態空間モデル（Koyama et al., 2007; Koyama et al.,and Shinomoto, 2007, 2010），また Kass はスプライン平滑化を用いた発火率推定を提案しているが（Kass, 2008），こ2010）
れらはいずれもポアソン性を仮定している．
一方で近年，一見ランダムだと思われてきた大脳皮質のスパイク発火の時系列パターンにも細胞ごとに固有な非ランダム性・非ポアソン性があり，大脳皮質の機能とも相関していることが発見された（Shinomoto et al., 2005; Shinomoto et al., 2009; Maimon and Assad, 2009）．このようなスパイク発火の非ポアソン性の適切な特徴付けが，発火率推定の精度を改善することが†

統計数理研究所：〒190–8562 東京都立川市緑町 10–3

統計数理

174

図 1.

第 60 巻

第 1 号 2012

スパイク列からの発火率推定の例．灰色太線は真の発火率，実線は推定発火率を表す．発火率は図中のラスタープロットで示されたスパイク時系列から推定した．スパイク時系
列は，真の発火率を用いて時間リスケーリングされたガンマ過程から生成した（Berman,1981）．（a）はスパイク間隔の逆数を取った素朴な発火率推定であり，（b）はポアソン過程を想定して（つまり，スパイク間隔分布が指数分布にしたがうと仮定して）本稿で提案する方法で平滑化を行った結果である．（c）はスパイク間隔分布がガンマ分布にしたがうと仮定して平滑化を行った結果である．パラメータの推定値は γ̂ = 0.42，η̂ = −1.09，φ̂ = 5.84 である．

報告されている．Omi らは，スパイク数データの非ポアソン性（スパイク数の平均と分散が一致しないこと）を考慮に入れることで，ヒストグラム法のビン幅が最適化されることを示した（Omi and Shinomoto, 2010）．また，Koyama and Shinomoto（2005）および Cunningham et al.
（2008）は非ポアソン発火をガンマ分布で表現し，発火率の平滑化にガウス過程を用いたベイズ
推定方法を提案した．
本稿では，状態空間モデルを用いた発火率推定法を提案する．状態空間モデルの利点は，ス
パイク発火の統計的モデリングから状態の平滑化，モデルパラメータの推定を統一的に扱えることにある．Koyama and Shinomoto（2005）および Cunningham et al.（2008）の提案した方法
も状態空間モデルに含まれるが，本稿で提案する枠組みにおいては，スパイク発火の非ポアソ
ン性を指数型分布族で表現することで，より幅広いモデリングを統一的に扱うことが可能となる．図 1 に，発火率推定において適切なスパイク発火のモデリングと平滑化が必要であることを模式的に示した．図 1（b）は，スパイク発火の非ポアソン性を考慮せずポアソン過程を仮定
して平滑化した推定発火率であるが，推定した発火率の変動はもとの発火率変動よりも小さく見積もられてしまっている．同じ発火率のもとでポアソン過程のほうが非ポアソンよりも変動を許すために本来の発火率変動の一部がポアソン性で説明されてしまうためである．図 1（c）
は，スパイク発火の非ポアソン性を適切にモデリングして平滑化したものである．（b）に比べてもとの発火率をよりよく推定していることが見てとれる．
2 節では提案する手法の詳細を解説し，3 節ではシミュレーションで生成したスパイク列デー
タに当てはめて提案手法の有効性を示す．
2.

手法

与えられたスパイク時系列 {t0 , t1 , . . . , tn } から発火率を推定する問題を考える．推定のため，スパイク間隔 yi = ti − ti−1 (i = 1, 2, . . . , n) の分布 p(yi |µi , η) を導入する．ここで，µi = E(yi ) は
i 番目のスパイク間隔 yi の平均パラメータであり，η は分布の形をコントロールするパラメー
タである．i 番目のスパイク間隔における発火率は λi = 1/µi で与えられる．推定発火率の平滑化を行うため，平均パラメータ µi = g(xi ) を状態 xi の関数とし，状態空間上で平滑化事前分布

状態空間法による神経スパイク発火率の推定175p({xi }ni=1 |γ) を導入する．ここで γ は滑らかさの度合いを表す正則化パラメータである．推定問題は状態の事後分布を求めることに還元され，状態の事後推定から発火率 {λ1 , λ2 , . . . , λn } の推定値が求まる．以上が概略であるが，以下で順を追って詳しく見ていく．
2.1 状態空間モデル
スパイク間隔の統計モデルとして，指数型分布族(2.1)

p(yi |θi , φ) = exp[φ{θi yi − b(θi )} + c(yi , φ)]を考える（Jørgensen, 1987; McCullagh and Nelder, 1989）．ここで θi は自然パラメータで，以下で見るように平均パラメータ µi の関数である．1/φ は拡散パラメータ，b(θ) はキュムラント
関数と呼ばれる．b(θ) は凸関数であると仮定する．yi の平均 µi と分散は(2.2)

µi = E(yi ) =∂b(θi ),
∂θi

Var(yi ) =1 ∂ 2 b(θi ) 1≡ V (µi )φ ∂θi2φとなる．ここで V (µ) は平均 µ と分散の関係を表す因子で，分散関数と呼ばれる．yi の分散がφ でスケールされていることからわかるように，φ はスパイク間隔の乱雑さの度合いを表すパラメータである．b(θ) の凸性から µi は θi の単調関数になるので，θi と µi は 1 対 1 写像によって変換される．このことを明記するときは θ(µi ) と書くことにする．
（2.1）
は φ についても指数型であり，観測データからの拡散パラメータ φ の推定のため，分布φ がその自然パラメータであると仮定する．すなわち，c(yi , φ) = φa(yi ) − ϕ(φ) + d(yi ) という形を持ち，分布は(2.3)

p(yi |θi , φ) = exp[φT (yi , θi ) − ϕ(φ) + d(yi )]という形に書き直すことができる．ここで(2.4)

T (yi , θi ) = θyi − b(θi ) + a(yi )は φ の十分統計量であり，ϕ(φ) は φ に関してのポテンシャルである．T (yi , θ) の期待値はポテ
ンシャルの 1 回微分(2.5)

η = E[T (yi , θ)] =

∂ϕ(φ)
∂φ

で与えられる． ϕ(φ) の凸性を仮定すれば φ と η は 1 対 1 対応であり，このことを明記すると
きは φ(η) と書く．
以上で見たように，分布は自然パラメータ (θi , φ) もしくは期待値パラメータ (µi , η) で表すことができるが，文脈から明らかな場合は特に断らずに期待値パラメータを用いた場合もp(yi |µi , η) = p(yi |θi , φ),T (yi , µi ) = T (yi , θi )などと表示することにする．ここで導入した指数分布族は，スパイク間隔分布としてよく用いられるガンマ分布，逆ガウス分布および対数正規分布を含む．これらの具体的な分布については 2.5 節にまとめる．
状態空間上での平滑化を行うために状態変数 x1:n = {x1 , x2 , . . . , xn } を導入し，yi の期待値
µi = g(xi ) を xi の単調関数とする．
（一般化線型モデルでは通常，g の逆関数 g −1 (µ) をリンク関数と呼ぶ
（McCullagh and Nelder, 1989）．）すると， i 番目のスパイク間隔 yi の観測モデルは
g を通して xi の条件付き確率分布(2.6)

p(yi |xi , η) = p(yi |g(xi ), η)

統計数理

176

第 60 巻として与えられることになる．
平滑化のための状態モデル

xi − xi−1 =第 1 号 2012

γ(yi + yi−1 )i2を導入する．i は平均 0，分散 1 の各時刻で互いに独立な正規分布にしたがうとする．ここで時間軸が {y1 , y2 , . . . , yn } で不等間隔に離散化されているため，時間軸上で一様に平滑化を行う
ために状態モデルに yi を含めた．確率分布で書き表すと，(2.7)

p(x1:n |γ) = p(x1 )n

p(xi |xi−1 , γ)i=2となる．ここで p(x1 ) は初期分布であり，
p(xi |xi−1 , γ) = (xi − xi−1 )21
exp −γ(yi + yi−1 )πγ(yi + yi−1 )は状態モデルの遷移確率である．
以上，観測モデル（2.6）と状態モデル（2.7）によって状態空間モデルが定義された．推定問題
は，スパイク間隔列 y1:n = {y1 , y2 , . . . , yn } が与えられたときのモデルパラメータ (η̂, γ̂) の推定，
および状態 x1:n の事後推定に還元される．モデルパラメータは対数周辺尤度関数


(2.8)
L(η, γ; y1:n ) = log · · · p(x1:n , y1:n |η, γ)dx1 · · · dxnを最大化するものに選ぶ．ここでp(x1:n , y1:n |η, γ) =

 n



p(yi |xi , η) p(x1:n |γ)i=1は観測と状態の同時分布であり，
（2.6）
式と
（2.7）
式から求まる．このように，対数周辺尤度を最大化することによってパラメータ選択をする処方箋は経験ベイズ法と呼ばれる（石黒 他，2004）．
モデルパラメータの値 (η̂, γ̂) が決まれば，状態の推定は事後分布(2.9)

p(x1:n |y1:n , η̂, γ̂) =p(y1:n |x1:n , η̂)p(x1:n |γ̂)p(y1:n |η̂, γ̂)から求まる．
以下で，具体的な状態の事後分布の計算とパラメータ推定のアルゴリズムを構築する．
2.2 事後分布のラプラス近似
観測値 y1:n が与えられたときの状態 x1:n の事後分布（2.9）を，そのモードを中心とする正規分布で近似する．このような近似はラプラス近似と呼ばれる（Kass, 1997）．ここでは，ひとまずモデルパラメータ (η, γ) の値は既知であるとし，これらの推定については次節で考えることにする．
まず，事後分布のモードを求める．状態を列ベクトル x = (x1 , x2 , . . . , xn )T で表すことにする．ここで ()T はベクトルの転置を表す．事後分布のモードは，事後分布の対数
n(2.10)

nlog p(xi |xi−1 , γ) +l(x) = log p(x1 ) +i=2

log p(yi |xi , η) + const.
i=1を最大化する x̂ で与えられる．ここで，観測モデル（2.6）および状態モデル（2.7）は {xi } について対数凹関数であるから l(x) も対数凹関数であり，x は凸集合に含まれるから唯一の大域

状態空間法による神経スパイク発火率の推定

177

解 x̂ の存在が保証される（Paninski, 2004; Koyama and Paninski, 2010）．ここでは，この最大化をフィッシャーのスコア法でおこなう（Fahrmeir and Kaufmann, 1991）．まず，状態の初期分布 p(x1 ) には拡散事前分布（任意の平均で分散が十分に大きい正規分布）を用いることで，その寄与を取り除くことができる（Durbin and Koopman, 2001）．すると，l(x) の勾配ベクトル∇l(x) =

∂l(x) ∂l(x), ∂x2 , . . . , ∂l(x)
∂x1
∂xn

T

は

⎧ 2(x − x )φ(η){y1 − g(x1 )} dg(x1 )
2
1
⎪
,
+
⎪
⎪
⎪
γ(y
+
y
)
V (µ)
dx1
2
1
⎪
⎪
⎪
2(xi+1 − xi )φ(η){yi − g(xi )} dg(xi )∂l(x) ⎨ 2(xi − xi−1 )
+
+
,
= −
⎪
γ(yi + yi−1 )γ(yi+1 + yi )
V (µ)
dxi
∂xi
⎪
⎪
⎪
⎪
⎪
⎪
⎩ − 2(xn − xn−1 ) + φ(η){yn − g(xn )} dg(xn ) ,γ(yn + yn−1 )
V (µ)
dxn

i=1
i = 2, 3, . . . , n − 1i=nとなる．また，H obs = ∇∇ log p(x|γ) を状態モデル（2.7）の対数のヘッセ行列，J state (x) = −E[∇∇ log p(y1:n |x, η)]


 2
∂
log p(yi |xi , η)= diag − E∂x2iを観測モデル（2.6）のフィッシャー情報行列とすると， J(x) = −H obs + J state (x) は⎞⎛D1 B1,2
0
···
0
⎟
⎜
..
..
⎟
⎜
.
⎟
⎜ B1,2 D2
B2,3
.
⎟
⎜
⎟
⎜
.
.
.
..
..
..
J(x) = ⎜ 0
⎟
0
⎟
⎜
⎟
⎜ .
..
⎟
⎜ .
.
Bn−2,n−1 Dn−1 Bn−1,n ⎠
⎝ .
0 ···
0
Bn−1,n Dnとなる．ここで
⎧
2

⎪
⎪
2
φ(η) dg(x1 )
⎪
⎪
+
,
i=1
⎪
⎪
V (µ)
dx1
⎪
⎪ γ(y2 + y1 )
⎪
⎪
⎪
2

⎪
⎨
2
φ(η) dg(xi )
2
Di =
+
+
, i = 2, 3, . . . , n − 1⎪γ(yi + yi−1 )γ(yi+1 + yi )
V (µ)
dxi
⎪
⎪
⎪
⎪
⎪
2

⎪
⎪
⎪
φ(η)
2
)
dg(x
⎪
n
⎪
,
i=n
⎪
⎩ γ(yn + yn−1 ) + V (µ)
dxn
2γ(yi+1 + yi )である．∇l(x) と J(x) を用いると，Fisher のスコア法を構成できる．今，m 反復での状態の推定値を x(m) とする．このとき，m + 1 反復での推定値をBi,i+1 = −x(m+1) = x(m) + J(x(m) )−1 ∇l(x(m) )で更新する．これを推定値が収束するまで繰り返す．
このようにして求まった事後分布のモード x̂ を平均ベクトルとし， J(x̂)−1 を共分散行列とする多変量正規分布で状態の事後分布を近似する．また，発火率 λi = 1/g(xi ) の MAP 推定値

178

統計数理

第 60 巻

第 1 号 2012と事後共分散行列は，状態 x のそれらから直接に λ̂i = 1/g(x̂i ) およびdg(x̂i )−1dg(x̂i )−1−1
diag
J(x̂) diag
dxi
dxi
と求まる．
以上は，
（2.1）
式で与えられる指数型分布族一般に対して事後分布の正規分布近似を導入したが，この分布族に含まれる特定の分布に適用する際に指定すべきモデル要素は，{g(x), V (µ), φ(η)}
である．ガンマ分布，逆ガウス分布，対数正規分布に対するそれぞれの具体形は 2.5 節で与える．
2.3 EM アルゴリズム
スパイク間隔列 y1:n からのパラメータ (η, γ) 推定は，先に述べたように（2.8）式で与えられる対数周辺尤度を最大化することで得られる．ここでは，対数周辺尤度関数の極大値を探索するために EM アルゴリズムを用いる（Dempster et al., 1977; 汪 他, 2003）．今 m 反復でのパラメータ値を (η(m) , γ(m) ) とする．まず，完全データ対数尤度関数の条件付き期待値
Q(η, γ|η(m) , γ(m) ) = E log p(x1:n , y1:n |η, γ)|y1:n , η(m) , γ(m)を計算し，これを (η, γ) について最大化したものをパラメータの更新値 (η(m+1) , γ(m+1) ) とす
る．このステップを値が収束するまで繰り返す．
われわれの状態空間モデルを代入して計算すると，まず η に対しては
n

 

∂φ ∂
∂
φE T (yi , µ(xi ))|y1:n , η(m) , γ(m) − ϕ(φ) + d(yi )Q(η, γ|η(m) , γ(m) ) =∂η∂η ∂φ i=1

 n


∂φ
E T (yi , µ(xi ))|y1:n , η(m) , γ(m) − nη=∂η i=1

ここで最後の式を得るために（2.5）
式を用いた．これを 0 とおくと，m + 1 反復でのパラメータ推定値は(2.11)

η(m+1) =n
1
E T (yi , µ(xi ))|y1:n , η(m) , γ(m)n i=1と求まり，η の十分統計量の条件付き期待値で与えられる．このように，パラメータ更新が十分統計量の条件付き期待値の計算に帰着されたのは，指数型分布族の期待値パラメータを扱っているためである．
同様に，γ に対しては∂Q(η, γ|η(m) , γ(m) ) = 0∂γを γ について解くことにより
n
E (xi − xi−1 )2 |y1:n , η(m) , γ(m)2(2.12)
γ(m+1) =n − 1 i=2
yi + yi−1と求まる．
（2.11），（2.12）式の中に現れる条件付き期待値の計算は，前節で導入した正規分布で近似した事後分布についておこなう．また，以上の式は指数型分布族一般に対するものであるが，特定の分布に当てはめる際に必要なモデル要素は T (y, µ) のみである．ガンマ分布，逆ガウス分
布，対数正規分布に対するそれぞれの具体形は 2.5 節で与える．

状態空間法による神経スパイク発火率の推定

179

2.4 モデル選択
2.3 節では，ひとつのモデルのなかでデータにもっとも適合するパラメータ値 (η̂, γ̂) を EM ア
ルゴリズムにより探索したが，2.1 節で述べたように，指数分布族（2.3）は特別な場合としてス
パイク間隔分布によく用いられるガンマ分布，逆ガウス分布および対数正規分布を含む．これ
らの分布のうちどの分布がもっともデータに適合するかを決めるため，各々の分布に対する対数周辺尤度（2.8）の値を比較し，もっとも大きい値を持つ分布を選ぶことにする（Akaike, 1980;
Mackay, 1992）（．Akaike, 1980 は負の対数周辺尤度にパラメータの次元数を加えた ABIC 最小化を提案したが，今の場合はパラメータの次元数は同じだから，周辺尤度最大化と同じである．）
対数周辺尤度値を計算するため，（2.8）式をL(η̂, γ̂; y1:n ) = log p(y1:n |η̂, γ̂)nlog p(yi |y1:i−1 , η̂, γ̂)

=
i=1
n

=


log

p(yi |xi , η̂)p(xi |y1:i−1 , η̂, γ̂)dxi

i=1と書き換える．この式から，対数周辺尤度の計算は i = 1, 2, . . . , n に対して
(2.13)p(yi |xi , η̂)p(xi |y1:i−1 , η̂, γ̂)dxi
を数値的に計算することに帰着されることがわかる．
（2.13）式は粒子フィルタを用いて計算することができる（Kitagawa, 1987; Doucet et al., 2001; 伊庭 他, 2005）．すなわち，粒子フィル
タにより N 個のサンプル x̃(j)（2.13）式は
i , j = 1, 2, . . . , N を p(xi |y1:i−1 , η̂, γ̂) から生成すれば，
Np(yi |xi , η̂)p(xi |y1:i−1 , η̂, γ̂)dxi ≈1(j)p(yi |x̃i , η̂)N j=1と近似できる．粒子数 N が十分に大きければこれはよい近似を与える．拡散事前分布を用いた場合，i = 1 の因子は 0 になるからこの寄与を取り除いて，i = 2, 3, . . . , n の和


n
N
1
(j)
L(η̂, γ̂; y1:n ) ≈log
p(yi |x̃i , η̂)
N j=1
i=2
で対数周辺尤度の値を計算すればよい．
(j){x̃i }N
j=1 (i = 1, 2, . . . , n) を生成する粒子フィルタのアルゴリズムを以下にまとめる．
(j)1） N 個の粒子 x̃1 （j = 1, 2, . . . , N ）を初期分布 p(x1 ) から生成し，i ← 1 とする．ここで，p(x1 ) は任意の平均で分散が十分に大きい正規分布にとる（拡散事前分布）．
(j)2） 重み wi （j = 1, 2, . . . , N ）
を
(j)
(j)
wi = p(yi |x̃i , η̂)

で求める．
(j)
(j)
3） 重みの和が 1 になるように規格化する． j = 1, 2, . . . , N に対して wi の確率で x̃i を選(j)び，全部で N 個のサンプルを選び出す．それを x̂i （j = 1, 2, . . . , N ）とする．
(j)(j)4） 各々の x̂i （j = 1, 2, . . . , N ）に対して x̃(j)i+1 を確率分布 p(xi+1 |x̂i , γ̂) から生成する．
5） i ← i + 1 として 2）に戻る．
，3）の重み付けはベイズの公式をシミュレーションするための手続きであ
ここで，ステップ 2）
る（Brockwell et al., 2004）．

統計数理

180

第 60 巻

第 1 号 2012

2.5 具体的なモデル例
ここでは以下の数値実験に用いる具体的なモデルを挙げる．アルゴリズムに必要なモデルパ
ラメータを表 1 にまとめる．

ガンマ分布．ガンマ分布は(2.14)

p(y|µ, κ) =κκ y κ−1 − κye µµκ Γ(κ)で与えられる．ガンマ分布をスパイク間隔分布として用いた解析は，例えば（Baker and Lemon,
2000）
を見よ．指数型分布族（2.3）式の形にするためには，
θ = −1/µ,b(θ) = − log(−θ),φ = κ,

ϕ(φ) = log Γ(φ) − φ log φT (y, θ) = θy + log(−θ) + log y = log

y
y
−
µ
µ

とすればよい．これからただちに分散関数はV (µ) =

∂ 2 b(θ)= µ2∂θ2と求まる．ポテンシャル関数 ϕ(φ) を微分することにより
η=

∂ϕ(φ)
= ψ(φ) − log φ − 1∂φdと，φ とその期待値パラメータ η との関係が得られる．ここで ψ(φ) = dφlog Γ(φ) はディガンマ関数である．φ(η) は η の関数として陽には求まらないので，上の式を φ について数値的に解く必要がある．µ を非負にするため，µ = g(x) = e−xとする．このとき，スパイク発火率は λ = 1/µ = ex で与えられる．
逆ガウス分布．逆ガウス分布はp(y|µ, ξ) =ξ(y − µ)2
ξ
exp −
2πy 3
2µ2 y

で与えられる．スパイク間隔分布として逆ガウス分布を用いた先駆的な研究として（Gersteinを挙げる．ここでand Mandelbrot, 1964）
1
ϕ(φ) = − log φ
2


1 y
µ
1
1
T (y, θ) = θy + (−2θ)1/2 −
=−
+
+
2y
2µ µ
y
µ

θ=−

1
,
2µ2

φ = −ξ,b(θ) = −(−2θ)1/2 ,表 1. アルゴリズムに必要なモデルパラメータ．

状態空間法による神経スパイク発火率の推定181とすれば，指数型分布族（2.3）式の形になる．ガンマ分布の場合と同様にして，分散関数はV (µ) =

∂ 2 b(θ)= µ3∂θ2と求まり，φ と η の関係は，
η=

∂ϕ(φ)
1
=−
∂φ
2φ

あるいは

1
2η
と求まる．µ を非負にするため，ガンマ分布の場合と同様に µ = g(x) = e−x に選ぶ．
φ(η) = −対数正規分布．対数正規分布は


1(log y − µ)2p(y|µ, σ ) = √
exp −
2σ 2
y 2πσ 22

(2.15)

で与えられる．対数正規分布をスパイク間隔分布として用いた解析は，例えば（Levine, 1991）
を見よ．指数型分布族（2.3）式の形にするためには，y を log y に置き換えて，µ = E(log y) = θ,

φ=

1
,
σ2

T (y, θ) = θ log y −b(θ) =

θ2,
2

ϕ(φ) =

1
2π
log
2
φ

(log y)2(y − µ)2
θ2
−
=−
2
2
2

とすればよい．分散関数は∂ 2 b(θ)=1∂θ2と求まる．ポテンシャル関数 ϕ(φ) を微分することにより
V (µ) =

η=

∂ϕ(φ)
1
=−
∂φ
2φ

または φ を η の関数として

12ηと求まる．µ = E(log y) は (−∞, ∞) に値を取るから非負に制限する必要はないので，単純にφ(η) = −µ = g(x) = xに選ぶことにする．このときスパイク発火率は状態 x の関数として

1
λ = 1/E(y) = exp − x −2φで与えられる．
3.

数値実験

3.1 順モデルで生成したデータからの推定
本節では，順モデル，すなわち推定するときに仮定するモデルで生成したスパイク時系列か
ら発火率を推定し，推定精度の評価を行う．

統計数理

182

図 2.

第 60 巻

第 1 号 2012

発火率推定の精度（MISE）を発火率周期 T の関数としてプロットした図．各点は，500
個のスパイクを含むデータから発火率を推定して MISE を計算し，さらに 100 回の試
行平均をとった値である．実線，破線，点線はそれぞれガンマ分布，逆ガウス分布および対数正規分布の結果を示す．

スパイク時系列 {ti } は，まず発火率が 1 に規格化されたスパイク間隔分布から n 個のスパ
イク間隔 {ui }ni=1 を用意し，変動発火率 λ(t) を用いた時間リスケール変換
 ti
ui =
λ(t)dtti−1により生成した．数値実験では2πtTを用いた．ガンマ分布，逆ガウス分布および対正規分布を用いて生成したスパイク列 {ti } に対し，同じモデルを仮定した推定方法を適用した．発火率推定の精度は，推定値と真の発火率との平均積分二乗誤差（Mean Integrated Squared Error, MISE）で評価した．
図 2 に数値実験の結果を示す．どのモデルの場合も変動周期 T が大きいほど推定精度がよいことがわかる．特に変動周期 T が小さいときは，速く変動する発火率を推定することが困難になり，MISE は急激に大きくなる．
λ(t) = 1 + 0.6 sin

3.2 積分発火モデルで生成したデータからの推定
本節では，積分発火モデル（Burkitt, 2006）
によって生成したスパイク時系列 {ti } から，スパ
イク間隔分布としてガンマ分布，逆ガウス分布および対数正規分布を用いて発火率を推定し，その結果を比較する．積分発火モデルは閾値下の細胞膜電位 V (t) の時間発展を記述するモデ
ルで，確率微分方程式


V (t)
+ µ + I(t) dt + σdW (t)(3.1)
dV (t) = −τ

にしたがい，V (t) が閾値 vth に到達するとスパイクを発生し vr にリセットされる．ここで W (t)は標準ウィナー過程である．τ は膜電位の減衰時定数であり，µ, I(t) はそれぞれ定常入力電流および変動入力電流，そして σ はノイズ強度である．それぞれのパラメータは物理的単位を持
r
つが，膜電位と時間をそれぞれ X = vVth−v，t ← τt と変換して無次元化することにより，閾値
−vr

状態空間法による神経スパイク発火率の推定

183とリセット膜電位はそれぞれ xth = 1，xr = 0 と無次元化される．入力電流とノイズ強度もそれぞれ√τ µ − vr
τ I(t) − vrσ τµ←, I(t) ←, σ←vth − vr
vth − vr
vth − vrと無次元化され，方程式は(3.2)

dX(t) = [−X(t) + µ + I(t)]dt + σdW (t)とスケーリングされる．これによって，モデルの振る舞いを決める本質的なパラメータは µ, I(t), σの三つに減らされる．
定常入力のみを加えた積分発火モデル（I(t) = 0）の定性的な振る舞いは，入力の値 µ によって大きく三つに分類できる（Wan and Tuckwell, 1982）．µ = 1 のとき（3.2）式の定常解の平均値は閾値 xth = 1 に一致するので「閾値（threshold）」モードと呼び，µ < 1 のときは「閾値下（subthreshold）」モード，µ > 1 のときは「閾値上（suprathreshold）」モードと呼ぶことにする．
これら三つの動作モードにおいて正弦波入力
2πt
I(t) = a sin
Tを加えて方程式（3.2）
を数値的に解いてスパイク列 {ti } を生成し，それに対しガンマ分布，逆ガ
ウス分布および対数正規分布を用いた推定法を適用した．発火率推定の精度は，推定値と真の発火率との MISE で評価した．積分発火モデル（3.2）の真の発火率は Koyama and Kass（2008）にしたがい数値的に計算した．
閾値下モードの結果を図 3（a）と表 2 左列に示す．周辺尤度値はガンマ分布と対数正規分布が逆ガウス分布よりもわずかに大きく，より適合するモデルとして選ばれる．発火率推定の
MISE の観点からは，対数正規分布を用いたときが他の分布を用いたときよりもわずかによい．
閾値上モードの場合（図 3（b）および表 2 中列）は，周辺尤度及び発火率推定の MISE の観点か
ら逆ガウス分布がもっともよい．概して閾値下と閾値上モードにおいては，三つの分布を用いた発火率推定にはそれほど大きな違いが見られない．
より興味深い結果は，閾値モードで見られる（図 3（c）および表 2 右列）．周辺尤度値，また発火率推定の MISE の観点からも逆ガウス分布がもっともよい．三つの分布を用いたそれぞれの推定発火率を見てみると，ガンマ分布の場合は変動がノイジーで平滑化が十分ではなく（図 3（c）左），対数正規分布を用いた場合は逆に平滑化されすぎてもとの発火率変動を追従できていない（図 3（c）右）．ガンマ分布を仮定した場合の発火率推定の過剰な変動は，他のふたつの分布に比べて分布の裾が薄く大きな変動を許さないために，データに含まれるスパイク間隔の変動の一部が発火率変動として説明されるためである．対数正規分布を仮定した場合は，逆に分布の裾が厚くスパイク間隔の大きな変動を許すため，発火率変動の寄与が小さく見積もられるためと考えられる．モデル分布の適切な選択が発火率推定の平滑化に大きく影響を与えているのが見てとれる恰好の例である．
4.

まとめ

本稿では，神経スパイク時系列から発火率を推定するための方法を提案し，シミュレーショ
ンにより生成したデータに当てはめて推定のデモンストレーションを行った．
ノイズを含む非定常スパイク時系列データから発火率を推定するためには平滑化が必要である．ここでは，発火率変動の滑らかさについての正則化を状態空間モデルの枠組みで取り入れることにより，推定問題は観測データが与えられたもとでの状態の事後分布を求めることに帰着された．そしてラプラス近似を用いることにより，状態の事後分布を近似的に計算する効率

統計数理

184

第 60 巻

第 1 号 2012

表 2.

数値実験の結果．“likelihood” は粒子フィルタで計算した周辺尤度の値． 105 個の粒子
を用いた．

図 3.

発火率の推定例．灰色太線が真の発火率，実線が推定発火率，鎖線が 95％ベイズ信用区
間である．ラスタープロットは推定に用いたスパイク列を表す．（a）は閾値下，
（b）は閾
値上，
（c）は閾値モードの結果である．ガンマ分布（左），逆ガウス分布（中），対数正規分
布（右）を用いたときの推定である．図中にパラメータ γ ，η および φ の推定値を載せた．

状態空間法による神経スパイク発火率の推定

185

的なアルゴリズムを構築した．
神経スパイク発火の非ポアソン性を特徴付けるために，ふたつのパラメータ（平均パラメータと拡散パラメータ）を持つ指数型分布族をスパイク間隔分布として用いた．この分布族は，スパイク間隔分布を記述するためによく用いられるガンマ分布，逆ガウス分布，対数正規分布を含む．したがって，われわれの提案する枠組みは，従来では個別に扱われてきたこれらの分布を状態空間モデルに系統的に取り入れることができる．この枠組みは West らにより提案さ，れた動的一般化線型モデルに含まれるが（dynamic generalized linear model, West et al., 1985）
ここでの新しい点は，スパイク間隔分布の拡散パラメータに関しても指数型を仮定することでEM アルゴリズムを簡潔に構成できたことにある．
最後に，我々の提案する手法を，推定で仮定するモデルおよび積分発火モデルから生成した
スパイク時系列に当てはめた．前者の数値実験では，発火率変動の周期が小さいと推定精度が急に劣化することをみた．スパイク時系列は時間軸上のスパースなデータであるから，発火率変動が速いほど推定するのは困難になる．後者の数値実験では，積分発火モデルの動作モードによりもっとも適合するスパイク間隔分布が異なることをみた．特に，積分発火モデルが閾値
モードにあるとき，逆ガウス分布がもっともよく適合した．他方でガンマ分布と対数正規分布を用いた場合には，推定した発火率はもとの発火率から大きくずれていた．スパイク間隔分布の選択が発火率推定に大きな影響をおよぼすことを顕著に示す例である．
したがって，与えられたスパイク時系列から発火率を推定するためにはスパイク間隔分布の
適切な統計モデリングと平滑化が必要である．我々の提案する枠組みは，スパイク間隔分布によく用いられるモデルを統一的に扱い，効率的な平滑化およびパラメータ推定アルゴリズムとモデル選択法を提供する．

参 考 文 献Adrian, E. D. and Zotterman, Y.（1926）. The impulses produced by sensory nerve endings: Part II:
The response of a single end organ, Journal of Physiology, 61, 151–171.
. Likelihood and Bayes procedure, Bayesian Statistics (eds. J. M. Bernard et al.),
Akaike, H.（1980）
1–13, University of Valencia Press, Valencia, Spain.
. Precise spatiotemporal repeating patterns in monkey primary
Baker, S. N. and Lemon, R. N.（2000）
and supplementary motor areas occur at chance levels, Journal of Neurophysiology, 84, 1770–
1780.
Berman, M.（1981）. Inhomogeneous and modulated gamma processes, Biometrika, 68, 143–152.
Bessou, P., Laporte, Y. and Pages, B.（1968）. A method of analysing the response of spindle primaryendings to fusimotor stimulation, Journal of Physiology, 196, 37–45.
. Recursive Bayesian decoding of motor cortical
Brockwell, A. E., Rojas, A. L. and Kass, R. E.（2004）
signals by particle ﬁltering, Journal of Neurophysiology, 91, 1899–1907.
. A review of the integrate-and-ﬁre neuron model: II. Inhomogeneous synaptic
Burkitt, A. N.（2006）
input and network properties, Biological Cybernetics, 95, 97–112.
. Inferring neural ﬁring rates from
Cunningham, J. P., Yu, B. M., Shenoy, K. V. and Sahani, M.（2008）
spike trains using Gaussian processes, Advances in Neural Information Processing Systems（eds. J. Platt, D. Koller, Y. Singer and S. Roweis）, 20, MIT Press, Cambridge.
. Theoretical Neuroscience, MIT Press, Cambridge.
Dayan, P. and Abbott, L. F.（2001）
Dempster, A. P., Laird, N. M. and Rubin, D. B.（1977）. Maximum likelihood from incomplete data
via the EM algorithm, Journal of the Royal Statistical Society:B, 79, 1–38.
）2001）. Sequential Monte Carlo Methods in Practice,
Doucet, A., de Freitas, N. and Gordon, N.（eds. （

統計数理

186

第 60 巻

第 1 号 2012

Springer-Verlag, Berlin.
Durbin, J. and Koopman, S. J.（2001）. Time Series Analysis by State Space Methods, Oxford University Press, Oxford.
. On Kalman ﬁltering, posterior mode estimation and Fisher
Fahrmeir, L. and Kaufmann, H.（1991）
scoring in dynamic exponential family regression, Metrika, 38, 37–60.
. Random walk models for the spike activity of a single
Gerstein, G. L. and Mandelbrot, B.（1964）
neuron, Biophysical Journal, 4, 41–68.
伊庭幸人，種村正美，大森裕浩，和合肇，佐藤整尚，高橋明彦（2005）. 『計算統計 II：マルコフ連鎖モンテカルロ法とその周辺』，統計科学のフロンティア，岩波書店，東京．
石黒真木夫，松本隆，乾敏郎，田邉國士（2004）. 『階層ベイズモデルとその周辺』，統計科学のフロン
ティア，岩波書店，東京．
Jørgensen, B.（1987）. Exponential dispersion models, Journal of the Royal Statistical Society:B, 49,
127–162.
. Laplace’s method, The Encyclopedia of Statistical Sciences, Update Volume 1,
Kass, R. E.（1997）(eds. S. Kotz, C. Read and D. Banks), 347–354, Wiley, New York.
. Adaptive spline smoothing of neural data, Neural Signal Processing: Quantitative
Kass, R. E.（2008）
Analysis of Neural Activity（ed. P. Mitra）, 35–42, Society for Neuroscience, Washington, D. C.
Kass, R. E., Ventura, V. and Cai, C.（2003）. Statistical smoothing of neuronal data, Network : Computation in Neural Systems, 14, 5–15.
. Non-Gaussian state-space modeling of nonstationary time series, Journal of
Kitagawa, G.（1987）
the American Statistical Association, 82, 1032–1063.
. Spike train probability models for stimulus-driven leaky integrateKoyama, S. and Kass, R. E.（2008）
and-ﬁre neurons, Neural Computation, 20, 1776–1795.
. Eﬃcient computation of the maximum a posteriori path andKoyama, S. and Paninski, L.（2010）
parameter estimation in integrate-and-ﬁre and more general state-space models, Journal ofComputational Neuroscience, 29, 89–105.
Koyama, S. and Shinomoto, S.（2005）. Empirical Bayes interpretations of random point events, Journal of Physics A: Mathematical and General, 38, L531–L537.
. Phase transitions in the estimation of eventKoyama, S., Shimokawa, T. and Shinomoto, S.（2007）
rate: A path integral analysis, Journal of Physics A: Mathematical and General, 40, E383–
E390.
Koyama, S., Eden, U. T., Brown, E. N. and Kass, R. E.（2010）. Bayesian decoding of neural spiketrains, Annals of the Institute of Statistical Mathematics, 62, 37–59.
Levine, M. W.（1991）. The distribution of the intervals between neural impulses in the maintaineddischarges of retinal ganglion cells, Biological Cybernetics, 65, 459–467.
. Beyond Poisson: Increased spike-time regularity across primate
Maimon, G. and Assad, A.（2009）
parietal cortex, Neuron, 62, 426–440.
. Bayesian interpolation, Neural Computation, 4, 415–447.
MacKay, D. J. C.（1992）
McCullagh, P. and Nelder, J. P.（1989）. Generalized Linear Models, 2nd ed., Chapman and Hall,New York.
. Optimizing time histograms for non-Poissonian spike trains,
Omi, T. and Shinomoto, S.（2011）
Neural Computation, 23, 3125–3144.
. Maximum likelihood estimation of cascade point-process neural encoding models,
Paninski, L.（2004）
Network : Computation in Neural Systems, 15, 243–262.
. Spikes: Exploring the
Rieke, F., Warland, D., de Ruyter van Steveninck, R. and Bialek, W.（1997）
Neural Code, MIT Press, Cambridge.
Shimazaki, H. and Shinomoto, S.（2007）. A method for selecting the bin size of a time histogram,状態空間法による神経スパイク発火率の推定

187

Neural Computation, 19, 1503–1527.
Shimazaki, H. and Shinomoto, S.（2010）. Kernel bandwidth optimization in spike rate estimation,
Journal of Computational Neuroscience, 29, 171–182.
. Estimating the ﬁring rate, Analysis of Parallel Spike Train Data（eds. S. Gruen
Shinomoto, S.（2010）, 21–35, Springer, New York.
and S. Rotter）
Shinomoto, S., Miyazaki, Y., Tamura, H. and Fujita, I.（2005）. Regional and laminar diﬀerences inin vivo ﬁring patterns of primate cortical neurons, Journal of Neurophysiology, 94, 567–575.
Shinomoto, S., Kim, H., Shimokawa, T., Matsuno, N., Funahashi, S., Shima, K., Fujita, I., Tamura, H.,
Doi, T., Kawano, K., Inaba, N., Fukushima, K., Kurkin, S., Kurata, K., Taira, M., Tsutsui,K., Komatsu, H., Ogawa, T., Koida, K., Tanji, J. and Toyama, K.（2009）. Relating neuronalﬁring patterns to functional diﬀerentiation of cerebral cortex, PLoS Computational Biology,
5, e1000433.
. The highly irregular ﬁring of cortical cells is inconsistent with
Softky, W. R. and Koch, C.（1993）
temporal integration of random EPSPs, Journal of Neuroscience, 3, 334–350.
Wan, F. Y. M. and Tuckwell, H. C.（1982）. Neuronal ﬁring and input variability, Journal of Theoretical Neurobiology, 1, 197–218.

汪金芳，田栗正章，手塚集，樺島祥介，上田修功（2003）. 『計算統計 I: 確率計算の新しい手法』，統計科学のフロンティア，岩波書店，東京．
West, M., Harrison, P. J. and Migon, H. S.（1985）. Dynamic generalized linear models and Bayesian
forecasting, Journal of the American Statistical Association, 80, 73–83.

188

Proceedings of the Institute of Statistical Mathematics Vol. 60, No. 1, 173–188(2012)

State-space Method for Estimating Neuronal Firing Rates
Shinsuke Koyama
The Institute of Statistical Mathematics

It is an important task in the field of neural coding to estimate firing rates fromsingle spike trains. Here, we propose a state-space method for this purpose. Specifically,
we develop an approximate method for computing posterior estimates of the state, an EM
algorithm for estimating the model parameters, and a model selection technique basedon marginal likelihood. We demonstrate the performance of the proposed method byapplying it to the simulated spike trains to show the validity of the method.

Key words: Neuronal spike trains, ﬁring rate estimation, state-space models.