情報論的学習理論テクニカルレポート 2009
Technical Report on Information-Based Induction Sciences 2009 (IBIS2009)

非定常２値時系列データから隠れた構造を読み取る- 神経科学データへの応用 Hidden Structures Detection in Nonstationary Binary Time Series Data
- Application to Neuroscientiﬁc Data 瀧山健∗
Ken Takiyama

岡田真人†
Masato Okada

Abstract: We propose the algorithm which can estimate event rate, timimngs of change
point and the number of states in a binary time series data. Event rate can be estimatedwith high accuracy even if the rate has nonstationary temporal correlation and mean value.
Our algorithm consists of nongaussian switching state space model, variational bayes, andlocal variational. We demonstrate the algorithm applying to neuroscientiﬁc data. Ouralgorithm can estimate nonstationary event rate, change point and the number of neuralstates using only one observation data, which is a neuroscientiﬁc requirement in somecases. Synthetic data analysis reveals our algorithm can discriminate neural states basedon both mean event rates and temporal correlations though many previous algorithm candetect change points based only on mean event rate. Our algotrithm can detect change
points in area MT neural data and estimate the number of neural states based on temporalcorrelation. These results indicate our algotihm is applicable to wide range of nonstationary
binary time series data, to neural data speciﬁcally.
Keywords: nongaussian switching state space model, variational bayes, local variational,
pruning, binary time series data, ﬁring rate estimation

1

まえがき
時系列データから変化点を検出することは機械学習,

データマイニングや統計学における重要な課題の１つである [1, 2]. 変化点とはその前後で時系列データの時間相関や平均値が急激に変化する点を意味し, 非定常性の原因となる.
時系列データ解析にはしばしば状態空間モデル (StateSpace Model(SSM)) が用いられる. SSM は観測された時系列データの生成過程, 観測過程をモデル化し, デー
タの予測や観測データのノイズ除去を行うことができる. しかし, SSM は 1 つの確率モデルを用いるため, 変∗ 東京大学大学院新領域創成科学, 277-8561 千葉県柏市柏の葉 51-5 tel.04-7136-4128, e-mail takiyama@mns.k.u-tokyo.ac.jp
Graduate School of Frontier Sciences, The University of Tokyo
Kashiwanoha 5-1-5, Kashiwa-shi,Chiba, 277–8561 Japan† 東京大学大学院新領域創成科学, 理化学研究所脳科学総合研究所
e-mail okada@mns.u-tokyo.ac.jp,
Graduate School of Frontier Sciences, The University of Tokyo,
Brain Science Institute, RIKEN

化点前後での時系列データの急激な性質の変化を捉えることはできない. SSM の拡張であり, 変化点を含む時系
列データを記述する確率モデルとして, スイッチング状
態空間モデル (Switching State Space Model (SSSM) )[3] が提唱されており, 臨床データ解析, 経済時系列デー
タ解析, 画像処理など幅広く応用されている [3, 4, 5].

SSSM では複数の生成過程を定義し, ある時刻の観測は
1 つの生成過程から行われるとモデル化する. 用いる生成過程がスイッチする時刻を変化点として検出し, 非定常なデータを定常なデータに分割する. そして各定常区
間毎に異なる生成過程のフィッティングを行い, 時系列
データの非定常性を記述する.
観測する時系列データが 0, 1 の 2 値のみをとる場面は数多く存在する. 単一分子間の癒着の発生 [6], 地震の発生 [7], ニューロンのスパイクの発生 [8] などが 2 値時
系列データの例であり, 各々において発生確率の推定手法が提唱されている. Brown らは地震データの解析に用入力

0...スパイクなし
1...スパイクあり
(a)

x

N
1

x

N
2

x

N
3

x

x

x

x

1
3

x

ラベル
番号

z

z

z

z

スパイク
データ

η

η

η

N
M

発火率
1
1

N試行目
発火率の推定値

ニューロン
スパイクを出力
010 … 110

1試行目
2試行目

10本

0

20本 27本

1000
2000
時間(msec)
(b)

1
2

1

2

3

η

1

1
M

2

M

3

M

(c)

図 1: (a): 入力が加えられた時, ニューロンはスパイクと呼ばれる 0, 1 の系列を出力として掃きだす. (b): 従来の発火率の推定手法. トライアル間で発火率は一定であると仮定し, 複数トライアルで得られたスパイク列をトライアル
平均することで発火率を推定する. (c): スイッチング状態空間モデルのグラフィカルモデル.

いられていた Kolmogorov-Smilnov テスト [7] をニューが望ましい [24]. しかしながら, 変化点を含む非定常な

ロンデータの解析へ応用できることを示した [9]. つま

発火率をワンショットから推定し, 更に神経状態数も同り, 2 値時系列データの解析手法は様々な分野に共通し時に推定する手法は未だ存在しない. 本手法をスパイクて適用可能であると考えられる.

データへ適用し, これらの問題点を同時に解決可能であ

本研究は時間変動し, 変化点を含む非定常なイベント

ることを示す.

発生確率, そして変化点の時刻を 2 値時系列データから本稿は以下の内容から構成されている. ２節では確率推定する手法を構築する. 本手法は非ガウス SSSM からモデルの定義, ３，４節では各々変分ベイズ法, EM アル

構成され, 推定, 学習法として変分ベイズ法 [10, 11], 局

ゴリズム [25] の構築を行う. ５節で本アルゴリズムを人所変分法 [12, 13] を用いる. 変分ベイズ法による関連度

工データへ適用し, 従来の手法との性能を比較する. そ自動決定 [14] により, 同時に最適な生成過程の数を推定して実データへ適用する. ６節でまとめと議論を行う.

を行う. そして 2 値時系列データの 1 例であるニューロンのスパイクデータ解析へ適用する.

Abeles らは遅延反応課題時のサルの前頭葉ニューロ
ンデータを解析した結果, ニューロンは異なる複数の神

2

確率モデル

2.1

スイッチング状態空間モデル (SSSM)

状態は異なる入力を反映し, 変化点前後で起こる神経状

SSSM を構成する生成過程の数を N とし, 各々の生成過程に事前分布を定義する. 以下では事前分布の数 N態遷移は, 入力の切り替わりやオンオフと対応しているをラベル数と呼ぶ. 時刻 m の観測が n 番目の事前分布

経状態間を遷移していることを示した [15]. 異なる神経[15, 16, 18, 19]. つまり, 変化点の検出, 神経状態数の推定は未だ明らかでないニューロンの情報表現の解明に重要な役割を果たす.
ニューロンはスパイクと呼ばれる０，１の２値系列をnnから出力されるなら zm= 1, そうでないなら zm=0を
1
N, ..., zM} を定とるラベル番号 z = {z11 , ..., z1N , ..., zM

義する. M は観測の全時刻である. n 番目のラベルの発n
火率を xn = {xn
1 , ..., xM } (n = 1, ..., N ) , 観測された出力として掃きだす (図 1(a)). また, 同じ実験を行ってスパイクの時刻を s とおき, 同時確率分布を以下の式でも出力するスパイク列は毎回異なり, 非常に不規則な振定義する.

る舞いを示す. 神経活動は不規則だが, スパイクの頻度p(s, x, z) = p(s|x, z)p(z)p(x1 )...p(xN )(1)である発火率によりニューロンは情報を表現していると考えられている [24]. 通常, 発火率は数多くのスパイク列

SSSM のグラフィカルモデルを図 1(c) に示す.

を平均して計算される (図 1(b)). しかし, 意志決定や運動計画などの内的な過程では, 実験毎にニューロンの神
経状態遷移時刻が異なり, １本のスパイク列 (ワンショット) から発火率と変化点を推定する必要がある [20, 21].
また, 神経活動を基に義手などのデバイスを操作するブ
レインマシンインターフェース [22, 23] では, 1 回の試
行毎に実験環境が変わりやすく, ワンショットでの推定

2.2

事前分布

ラベル番号, 各ラベルの発火率の事前分布を各々p(z 1 |π) =∏N

n=1 (πp(z m+1 |z m , a) =

n z1n

)

∏N
n=1

(2)
∏N

n k
zm+1
nk zm
k=1 (a )(3)√p(xn |β n , µn ) =となる. ただし, η̂m = 2|β n Λ|
βn
n
n T
n
n
exp(−
(x
−µ
)
Λ(x
−µ
))
2
(2π)N

(4)
と定義する. π n , ank は各々, 初期時刻でラベル n を選択する確率, ラベルが n から k へ切り替わる確率であ∑N∑N
る ( n π n = 1, k ank = 1). β n , µn は各々n 番目のラベルの発火率がもつ時間変動の滑らかさ, 平均である. Λ の構造は発火率の時間相関を表わし, ここでは∏np(xn ) ∝ m exp(− β2 (xm − xm−1 − µm )2 ) を満たす Λ3

隠れ変数 x, z の推定とパラメータ a, π の学習を行う. SSSM における周辺事後分布 p(・|s) の計算は困難であり [3], この計算困難を回避するため変分ベイズ法を用いる. 変分ベイズ法とは, 隠れ変数を w, パラメータを θ とすると

F[q]からラベル k への遷移のしやすさを制御するハイパーパ
ラメータである.

2.3

dwdθq(w)q(θ) log

q(w)q(θ)
p(s, w, θ)(9)と定義される変分自由エネルギーを最小にする近似事後
分布 q(・) を求める手法である. ただし,∫∫
U[q] = −dwdθq(w)q(θ) log p(s, w, θ),∫∫H[q] = −dwdθq(w)q(θ) log q(w)q(θ) であり, 各々
統計力学における内部エネルギー, エントロピーと対応する量である. この手法では尤度関数と事前分布が互いに共役分布でなくてはならない. しかし, (4), (8) 式は発火率に関し共役性がない. そこで我々は局所変分法[12, 13] を用い, 非ガウス SSSM に対する変分ベイズ法の実行を可能とした.
nを用いて (8) 式を局所変分法は変分パラメータ ξm

尤度関数

観測時間 T は K 個の長さ ∆ のビンから構成され (∆ ¿1), 各々のビンには高々1 本のスパイクが含まれるとする. k 番目のビンにスパイクが観測されたなら +1, そうでないなら −1 となる変数 ηk を定義する (k = 1, ...,K). k 番目のビンでの発火率を λk とし,
1+ηk
1−ηk
∏K
p(s|λ) = k (λk ∆) 2 (1 − λk ∆) 2

∫ ∫== U[q] − S[q]ディリクレ分布

刻でラベル n の選ばれやすさ, 時刻間におけるラベル n

u=1 η(m−1)r+u − C である.

非ガウス SSSM の推定と学習を用いる. 更に π, a の事前分布を, (2), (3) 式に共役な∏Nnp(π|γ n ) = C(γ n ) n=1 (π n )γ −1
(5)
∏
∏
nk
N
N
p(a|γ nk ) = n=1 [C(γ nk ) k=1 (ank )γ −1 ] (6)
∑N n
Γ(
γ )
nk
n=1
) =
で定義する. ただし C(γ n ) = Γ(γ 1 )...Γ(γN ) , C(γ
∑N nk
Γ(
γ )
k=1
であり, C(γ n ), C(γ nk ) は各々p(π|γ n ),Γ(γ n1 )...Γ(γ nN )p(a|γ nk ) の規格化定数である. γ n , γ nk は各々初期時∑rpξ (s|x, z) ∝

∏M,N

n
Ln
n
n 2 zmmm,n [exp(− 2 (xm − η̃m ) )](10)として (4) 式の共役分布に変形する手法である. ただし,nn). (10) 式を用いた変η̃m= η̂m /Lnm , Lnm = ξCn tanh(ξmm
分自由エネルギー Fξ [q] は Fξ [q] ≥ F[q] を満たすが, EM

アルゴリズム [25] により ξ を推定することで q[・] に最(7)として尤度関数を定義する [26, 27]. 発火率とは時間ビンの長さをかけたものがニューロンのスパイクイベントの発生確率に相当し, 点過程においては条件付き集中度, 生
存解析においてはハザード関数と呼ばれている. ∆ ¿ 1
より λk ∆ ∈ [0, 1) である. 発火率が非負であることを厳λk ∆密に扱うため, exp(2xk ) = 1−λ(xk ∈ (−∞, ∞)) とk∆

してロジット変換を行う. 以下では xk を発火率と呼ぶ.

も近い分布を計算できる [28].
変分ベイズ法は多数の局所解が存在するため, 変分ベ
イズ法におけるアニーリング法 [29] により最適解を探索する. 温度パラメータ T を用いた変分自由エネルギーはF[q; T ] = U[q] − T S[q](11)で定義され, 温度 T を高い値から徐々に 1 に近づけていくことで, 大域解を得ることができる. 隠れ変数 xn , z∆ ¿ 1 より, K は大きい値をもち, 計算量やメモリーの推定値は近似分布√|W n |
1nq(x ) =
exp(−(xn − µ̂n )T W n (xn − µ̂n ))
(2π)N
2T

の蓄積に困難が伴う. そのため, 観測時間の粗視化を行(12)

本稿ではワンショット推定の定式化を行うが, 複数試行への拡張は容易である.

う. 観測時間 T は M 個の離散化された幅 C のビンから構成され, r∆ = C とする. この粗視化されたビン内で発火率は一定の値を持つとし, 多数のスパイクを含む.
最終的にラベル番号を含めた観測モデルは∏N,Mnp(s|x, z) = n,m [exp(η̂m xnm −C log 2 cosh xnm )]zm (8)

q(z)

∝

N
∏
n=1

exp(

N
M
b̂n n
π̂ n z1n ∏ ∏)exp( m )zm
TT
n=1 m=1

M
−1 ∏
N ∏
N
∏
m=1 n=1 k=1ânk n kexp( m )zm zm+1T(13)に従う. ただし,

W = CL + hβ iΛ
n

n

n

n −1

n

µ̂ = (W )(14)(w + hβ iΛµ )nnnπ̂ = hlog π innb̂nm = η̂m hxnm i−

(15)
(16)

C n
n 2
n
l(ξ )(h(xnm )2 i−(ξm) )−C log 2 cosh ξm
2 m
(17)

ânk = hlog ank inであり, Ln は hzm
i

(18)

n
tanh ξmn
ξm

を (m, m) 成分にもつ対角行niη̂m をもつベクトルである.
列, wn は (1, m) 成分に hzm

h·i は近似分布による平均を意味する.
パラメータ π, a は近似分布q(π) ∝

N∏(π n )

hz n i+γ 1 −1+T
1
−1
T

n=1

q(a) ∝N ∏N∏(ank )(19)(20)

アルゴリズム [30] から計算される.
温度の更新側は Ti+1 = 12 Ti + 12 (T1 = 100) を用いる

[3]. i とは更新回数である.

[8] (κ = 2.4) からスパイクを生成し, 生成したスパイク
データを用いて本アルゴリズムにより発火率を推定した(図 2(a)). スパイク発生時刻の時間間隔はガンマ分布に従うことが数多く報告されており [31], 実データ解析を想定しガンマ過程からスパイクを生成した. 入力が加わると発火率に不連続な変化が生じると考えられている[32, 33]. 過渡的な入力, 持続的な入力がニューロンに加えられる状況を想定し, 発火率を推定した (図 2(a)). 入

可能であるか否かを検証する.

2.5(n 6= k) を用いた. 推定結果を図 2 に示す. 図 2(a) 上
図より, ワンショットから高精度で非定常な発火率を推定していることがわかる. 図 2(a) 下図は, 各時刻 m にn
おけるラベル番号 n の近似事後分布による平均 hzmiを
ni は 1 に近い値を持つ.
表しており, 白色が濃いほど hzm

発火率の推定値を示している. ラベル 1, 2, 3 は各々想定された自発活動, 過渡的な入力, 持続的な入力が加わ変分パラメータ ξ, 発火率の事前分布に含まれる µ, β
を EM アルゴリズムにより更新する. ただし, Q 関数

(21)における平均操作 h·i は, パラメータ θ(t) を用いた近似
分布により行う. ただし, θ = {ξ, µ, β} であり, t は更新回数を意味する. 本来の EM アルゴリズムでは平均操作
は事後分布により行われるが, SSSM では事後分布による平均が計算困難であるため [3] 近似分布により平均操作を行う. (21) 式で与えられる Q 関数は (9) 式に含まれる U[q; θ, θ t ] の符号反転したものに等しく, (21) 式の最大化は変分自由エネルギー (9) 式の最小化に等しい.
∂Q(θ|θ (t) )= 0 を解くことによりパラメータの更新側∂θnξm= h(xnm )2 i

(22)µn = hxn i

(23)

M
β =
Tr[Λ((Wn )−1 + (hxn i − µn )(hxn i − µn )T )](24)nを導出することができる. EM アルゴリズムと変分ベ
イズ法の反復計算により, 変分自由エネルギーを最小化する.

任意に設定した発火率から人工的に非斉次ガンマ過程

図 2(d) は削除されたラベルを除き, 各ラベルが推定した

EM アルゴリズム

Q(θkθ ) = hlog p(η, x, z|θ)i

人工データ

以下では各パラメータとして T = 4.0, ∆ = 0.001, C =0.04, γ n = 1(n = 1, ..., 5), γ nk = 100(n = k), γ nk =に従う. z の１次統計量, ２次統計量は Forward-Backwardt

5.1

タから発火率の推定, 変化点の検出, ラベル数の推定がn=1 k=14

結果

力の時刻が未知であるとき, 本手法はワンショットデー

∑M −1

n zk
hzm
i+γ nk −1+T
m+1
m=1
−1
T

5

る状態の発火率を推定し, これらを組み合わせ非定常な発火率を推定している. 変分ベイズにおける関連度自動決定により不要なラベルであるラベル 4, 5 は削除され,この人工データにおいてラベル数は 3 と推定された.
関連度自動決定によるニューロンのラベル数推定の精度を確認するため, 変分自由エネルギーを用いたモデル
選択 [11] を行った. 図 2(e) は各ラベル数に対する変分
自由エネルギーの値を示している. ラベル数 3 が最も低い値を示し, モデル選択されていることがわかる. 図2(b) はラベル数 7 とラベル数 3 の場合の発火率の推定
nn値, hzmi の比較である. ラベル数が多いほど hzm
i の曖昧さは減っているものの, 推定発火率, 変化点, 削除されず有効的に用いられるラベルの数は, ほぼ一致している
ことがわかる.

5.2

他の手法との比較

従来の発火率推定との比較を行う. 比較する手法はカーネル平滑化 (Kernel Smoothing(KS))[34], カーネル
幅最適化 (Kernel Band Optimization(KBO))[33], 適応
カーネル平滑化 (Adaptive Kernel Smoothing(KSA))[35],ベイジアン適応的回帰スプライン (Bayesian Adaptive

MT 野のデータ
300

120

ラベル数3での推定値
ラベル数7での推定値
真の発火率

200

0
0

1000

2000

時刻(msec)

3000

4000

0

80

発火率 (Hz)040

変分自由エネルギー

ラベル 1 の推定値
ラベル 2 の推定値
ラベル 3 の推定値
真の発火率

1000

2000

時刻(msec)

3000

0
1000

2000

時刻(msec)

3000

0

4000

0

1000

4000

487
486
485
484
483
482
481
4801

3

(d)

4

ラベル数

4000

3000

各ラベルの推定値

変分自由エネルギー

2

2000

時刻(msec)

(c)

(b)

(a)
各ラベルの推定値

00< z >1

1
2
3
4
5

ラベル番号

0

ラベル番号

1
2
3
4
5

K-Sプロット

100

40

40

<z>
1

３ラベル ７ラベル< z >1

推定した発火率

発火率 (Hz)80

発火率 (Hz)120

推定した発火率
真の発火率80

発火率 (Hz)

人工データ

5

6

ラベル 2 の推定値
ラベル 4 の推定値

7

0

500

1000
1500
時刻(msec)

2000

(f)

(e)

図 2: 本アルゴリズムを人工データ, 実データに適用した結果. (a): 上図は人工データから推定した発火率と真の発ni の値. ただし, h·i
火率. 横軸は時間 (msec), 縦軸は発火率 (Hz) である. 下図は各時刻, 各ラベル毎に計算した hzm

は近似分布による平均を表す. (b): ラベル数 3, ラベル数 7 での発火率の推定値. (c): MT 野のデータから推定した発火率. 枠内は Kolmogorov-Smilnov テスト (K-S テスト) の結果であり, 灰色の潜が信頼区間, 黒色の線が推定結果を表す. (d): (a) の結果における各ラベルの推定値. 刈り取られたラベルの推定値は表示していない. (e): (a) に示したスパイク列を用いた推定における, 各ラベル数に対する変分自由エネルギーの値. (f): (c) の結果における各ラベ
ルの推定値. 刈り取られたラベルの推定値は表示していない.

Regression Spilne(BARS))[36], ベイジアンビンニング(Bayesian Binning)[37] である. KS にはガウシアンカー
ネル
f (x, y, σ) = √11
exp(− 2 (x − y)2 )
2σ
2πσ

(25)

を用い, カーネルの幅を決める σ には σ = 30(msec)(KS30),

Space Model(VB-SSSM)) と表記する. 各比較手法の利点, 欠点は Cunningham et al.[24] に述べられている.
平均絶対誤差 (Mean Absolute Error (MAE) ) を
1∑|λt − λ̂t |
T t=1
T

MAE =

(26)

σ = 50(msec)(KS50), σ = 100(msec)(KS100) を用いた.

として定義する. ただし, λt とは真の発火率, λ̂t とは推カーネル幅はヒューリスティックに決めなくてはならない

定発火率を表す. 図 2(a) に示した真の発火率から 7 本

ため, このように３種類のカーネル幅を用いて発火率を推のスパイク列を生成し, 各々のスパイク列毎に発火率推定した. また, KSA においてもヒューリスティックにカー定を行った. 7 本のスパイク列に対する各手法の MAE

ネル幅を決める必要があるが, ここでは KBO により推定の平均値, 標準偏差を図 3(a) に示す. 図 3(a) より, 非定されたカーネル幅を KSA に用いた. BARS, BB は各々

常な発火率のワンショット推定において本研究が高精度

ウェブ上 http://lib.stat.cmu.edu/ kass/bars/bars.html,

であることがわかる.

http://mloss.org/software/view/67/から入手可能なコー

図 3(b), (c) は従来スパイクデータから神経状態遷移

ドを利用した. 図中では, 本手法を変分ベイズスイッチン

時刻を検出するために用いていた, 多変数ポワソン分布

グ状態空間モデル (Variational Bayes Switching Stateを出力分布としてもつ隠れマルコフモデル [15, 16, 17,

0

120

10

80

(Hz)

40

(Hz) 30< z >1 000KS30 KS50 KS100 LOC KSA BARS BB VB-SSSM

mPHMM

2

mPHMM

4

0

60

本手法

6

真の発火率
変化点

90

本手法

平均絶対誤差< z >1

8

発火率

発火率

12

真の発火率

1000

(a)

2000
3000
時刻(msec)

4000

0

1000

(b)

2000
3000
時刻(msec)(c)

図 3: (a): 本手法と他の手法を図 2(a) で示されている発火率の推定に適用したときの推定誤差. 縦軸は推定発火率と真の発火率の平均絶対誤差を表わす. グラフは 7 本のスパイク列ごとの推定誤差の平均値 ± 標準偏差. 平均絶対誤差の定義は文中を参照. (b): 平均発火率が異なる複数の状態を含む発火率を, 本手法と多変数ポワソン分布を用いた隠れマルコフモデル (multivariate Poisson - Hidden Markov Model (mPHMM)) で区分した結果. (c): 時間変動の滑らかさが異なる複数の状態を含む発火率を, 本手法と mPHMM で区分した結果.

18, 19](以下, mPHMM と表記) による変化点検出性能

ラベル数

1

2

3

計

と, 本手法との変化点検出性能を比較した結果である.

データ数

5

20

3

28

本手法はワンショット, mPHMM は同一発火率から生成した 10 本のスパイク列を用いて変化点検出を行った. 図

表 1: 28 本のスパイクデータを解析し, 本手法により推3(b) より, 平均発火率の異なるラベルの変化点は両手法で検出可能である. 時間変動の滑らかさが異なるラベル

定されたラベル数.

の変化点検出を両手法で行った結果が図 3(c) である. 本k), γ nk = 2.5(n 6= k) を用いた.
Neural Signal Archive から入手できる MT 野データの内, nsa2004.1, j001 に含まれるコヒーレンス 6.4%の手法では変化点検出, ラベル数の推定が可能であること
がわかる.

5.3

実データ

内側側頭葉 (Medial Temporal Area(MT 野)) ニュー
ロンは, 視覚刺激の運動方向に優先的に反応する [39]. 本
節ではランダムドット刺激を提示時のアカゲザル MT 野
ニューロンデータの解析を行う. このデータはウェブ上の Neural Signal Archive から入手可能であり [40], 詳細
な実験環境は Britten ら (1992) に記述されている [41].
図 2(c) はランダムドット (コヒーレンス 6.4%) 提示
時の MT 野データ (nsa2004.1 j001 T2) の解析結果である. 図 2(c) 上図は推定発火率とコルモゴロフ-スミル
ノフプロット (Kolmogorov-Smirnov plot (K-S plot))[9]を表しており, このデータでは振動発火率が推定されている. K-S プロット中の灰色の線が 95%の信頼区間を示しており, プロットした直線がほぼ信頼区間に収まっ

ランダムドット提示時のニューロンデータ 28 本の解析を行った. 28 本のスパイクデータを解析し, 推定されたニューロンのラベル数を表１に示す.
解析した 28 本のデータのうち, 2 ラベルと推定された20 データにおける, 発火率の時間相関 β n および発火率の平均 µn を推定した. 神経状態 1 から神経状態 2 への変化点時刻, そして神経状態 1, 神経状態 2 の推定発火率に含まれる β, µ を図 4 に示す. (c), (d) は各々µt の∑Tn n
1
時間平均, 標準偏差を表しており, hµnt i = Tn
t=1 µt ,
∑Tn n
1
n
n 2
hdµ i = Tn t=1 (µt − hµt i) と定義される. ただし,

Tn とはラベル n から他のラベルへの変化点時刻, もし
くは全観測時間を意味し, 時刻 t は変化点を 0 と見なす.

6

議論ていることから, 推定発火率は信頼できると言える. 図

本研究は 2 値時系列データから, 平均, 時間相関が非2(c) 下図, 図 2(f) より, ラベル 2, 4 が各々入力直後の状定常であるイベント発生確率, 変化点, ラベル数を推定

態, 時間が経った後の状態を推定しており, その他の不するアルゴリズムを提唱した. 非定常なイベント発生

要なラベルは削除されていることがわかる. 以下では入確率, 変化点は SSSM により推定され, ラベル数は変分

力直後のラベル, 時間が経った後のラベルを各々神経状

ベイズにおける関連度自動決定により推定可能であった

態 1, 神経状態 2 と呼ぶ. このとき T = 2.0, ∆ = 0.001(図 2(a), 2(c)).
MT 野の解析で抽出された神経状態 1, 神経状態 2 はであり, C = 0.02, γ n = 1(n = 1, ..., 5), γ nk = 100(n =

4000

5
10
15
トライアル番号

(a)

20

−1

2

−2

1.5

<dμ>

0.4
0.2
0
0

β

x 10 5
4
3.5
3
2.5
2
1.5
1
0.5
0
0
5
10
15
トライアル番号−3<μ>遷移時刻 (sec)

1.4
1.2
1
0.8
0.6

−4

0.5

−5
20

(b)

−6
0

1

5
10
15
トライアル番号

(c)

20

0
0

5
10
15
トライアル番号20(d)

図 4: 2 ラベルと推定された 20 本のデータの推定量. 横軸はトライアル番号を番号が小さい順に並べたものである.

(a): 神経状態 1 から神経状態 2 への変化点の時刻. (b): 神経状態 1(青), 神経状態 2(赤) の β の値. (c): 神経状態
1(青), 神経状態 2(赤) の µt の値を時間平均したもの. (d): 神経状態 1(青), 神経状態 2(赤) の µt の標準偏差.

各々過渡状態, 持続状態と呼ばれている状態 [42, 43, 44]

ライアル間でばらつきが大きい. そのため, 解析に用いに相当しているものと考えられる. 両状態での発火率のたデータはトライアル平均を行うことで変化点時刻の情強さは物体の速度に依存し [43], 過渡状態での発火率の報をあいまいにしてしまう. 解析に用いた MT 野データ

強さは物体の加速度に依存している [44]. 図 2 に示した

はワンショットデータから発火率, 変化点の時刻を検出

解析結果においては, 図 2(f) よりラベル 2 が過渡状態することが望ましいことを示唆する結果となった.

を, ラベル 4 が持続状態の発火率を推定している. 以下

図 3(b), (c) に示されるように, 本手法は従来の手法とでは神経状態 1 を過渡状態, 神経状態 2 を持続状態と呼同様に平均発火率に基づいた変化点検出に加え, 時間変ぶ. また, 各ラベルが各神経状態と対応していることか動の滑らかさに基づいた変化点検出を行えることを示し

ら, 以下ではラベルを神経状態と呼ぶ.

た. 図 4(c), (d) より平均発火率はトライアル間で変動

図 2(c) に示した MT 野の推定発火率は, 過渡状態, 持が大きいものの, 図 4(b) より時間変動の滑らかさはト

続状態の両神経状態において振動している. その理由

ライアル間で共通し, 更に過渡状態と持続状態とで大きは, バースト発火を示すニューロンを解析したためであな差がみられる. つまりランダムドット提示時の MT 野

る [42]. バースト発火が起きている時刻で高発火率, そ

データに含まれる状態を推定するには, 時間変動の滑らうでない時刻で低発火率と推定されている. 従来は刺激

かさに基づき変化点を検出できる本手法が適しているとを加えた後, 一定時間以内は過渡状態と見なす [42]. も考えられる. 実際, 多変量ポワソン分布に基づく隠れマ

しくはトライアル平均により求めた平均発火率からの差

ルコフモデルでは MT 野データの変化点は検出できなにより, 過渡状態から持続状態へ切り替わる変化点を検いことも確認している (非公開データ).

出していた [43]. これらはヒューリスティックな手法でニューロンは入力を加えられたときのみではなく, 自あり, 本手法はワンショットから統計的に 2 つの神経状

発活動時も神経状態遷移を示す [45, 46]. 決められたコー

態を区別することができる点が MT 野データの解析に用スの上を移動するトラッキング課題後の睡眠 (自発活動)

いた場合の特長の一つである. 統計的な手法である本手

中, 徐波睡眠状態でニューロンは低発火率の DOWN 状

法は他の領野のニューロンデータの解析にも用いること

態, 高発火率の UP 状態の 2 つの神経状態間を遷移する.

ができる. 注意の神経機構を解明するため, 近年 MT 野と頭頂間溝外側部との同時記録が盛んに行われている.

UP 状態は記憶の形成に深く関わると示唆されており,
自発活動中のスパイクデータから UP 状態と DOWN 状

本手法を各々の領野のスパイクデータに適用し, 変化点

態を区別することは記憶形成メカニズムの解明に重要

時刻, 神経状態数の推定結果から 2 つの領野間の相互作な役割を果たす. 近年 Chen et al(2009) は一般化線形用をワンショットデータから検証できる可能性がある.

モデル, 隠れセミマルコフモデルを用いて UP 状態と20 本のスパイクデータで推定された 2 状態が必ずしも全て過渡状態, 持続状態と対応しているとは限らない.

DOWN 状態を区別するアルゴリズムを提唱した [46].
しかし Chen らのアルゴリズムはサンプリング手法の一

しかし, Jones ら (2007)[17] が味覚野で示した結果に基つであるリバーシブルジャンプマルコフ連鎖モンテカル

づくならば, トライアル間で神経状態間の遷移時刻は異ロ法 [47] を用い, 計算時間がかかってしまう. さらに, 解なるが, 神経状態の系列はトライアル間で一致する. つ析する全てのニューロンの変化点時刻が同じであり, 神

まり, 状態 1 が過渡状態, 過渡状態から遷移する状態 2

経状態の違いは平均発火率の違いのみであるという仮定

は持続状態であるとみなすことができる. 図 4(a) で示を置いている. 従来のアルゴリズムも同様の仮定を用いした通り, 過渡状態から持続状態への変化点の時刻はトており [15, 16, 19, 17, 18], 本稿で用いた MT 野データ

のように発火率の時間変動の滑らかさの違いにより変化点を検出することは難しいであろう. 図 3(a), 図 3(b),図 3(c) より本手法は高い発火率推定精度, 変化点検出能力を示し, 更に図 2(b), 図 2(e) より神経状態数の推定も同時に行えることが, 先行研究と比較した際の本手法の特長である.

Abeles ら [15], Gat ら [16] は遅延到達運動時における前頭葉内の articulate sulcus に含まれるニューロンデー
タから, mPHMM を用いて神経状態遷移時刻の検出を行った. 結果としてデータは 6 個の神経状態に区分され, 各々の神経状態内でクロスコレログラムを計算し,ニューロン間の相互作用が異なることを発見した. この解析ではまず平均発火率に基づきデータを区分している. つまり本手法のように時間変動の滑らかさ, 時間
相関に基づきデータを区分することで彼らの結果とは異なる変化点時刻, 神経状態数が推定される可能性がある. よって, 本手法の多ニューロン解析への拡張が今後の課題の一つとして挙げられる. また, 近年の計測技術の向上により, 数１００オーダーのニューロン活動を同時に計測可能となった. ネットワーク全体においてニューロンの状態遷移が担う役割を解明するためにも, 本手法
を多ニューロン解析へ拡張する必要がある. 解析する
全てのニューロンの変化点時刻が同じであるという仮
定 [15, 16, 19, 17, 18, 46] を用いるならば, 本手法の多ニューロン解析への拡張は容易い. より詳細な多ニューロン解析を行うため, 発火率の事前分布にニューロン間相互作用を取り入れる, または観測モデルにニューロン間の相関を取り入れることが必要であると考えられる.
更に今後の展望の 1 つとして離散化したビン幅の最適化が挙げられる. Shimazaki ら (2007)[48] は刺激前後時
間ヒストグラムを用いた推定時に, 平均積分二乗誤差が最小になる離散化したビン幅を計算する手法を提案ていしている. その計算されたビン幅を用いて本手法を適用するということが考えられる. これは計算時間もかからず高い性能を示す可能性はある. ベイズ最適化の枠組み
で離散化したビン幅を推定する枠組みとして, Watanabe
ら [38] は周辺尤度最大化による最適なビン幅の推定方法を提案している. 本研究も同様に最適なビン幅を推定できる可能性がある.

参考文献[1] V. Guralnik and J. Srivastava, ACM-SIGKDD ,(1999) 33.
[2] J. Takeuchi and K. Yamanishi, IEEE Trans. Know.
Data Eng. , 18(4) (2006) 482.
[3] Z. Gharamani and G. E. Hinton. Neural Compt.
12, (2000) 831.
[4] Hamilton, J. D. Econometrica 57, (1989) 357.

[5] V. Pavlovic et al. Advances in NIPS 13, (2001)
981.
[6] Y. Hung et al. J. American Stat. Association
103(483), (2008) 1248.
[7] Y. Ogata. J. American Stat. Association 83(401),(1988) 9.
[8] R. Barbieri. J. Neurosci. Methods 105, (2001) 25.
[9] E. N. Brown et al. Neural Compt. 14 (2001) 325.
[10] H. Attias Proc. 15th Conf. on UAI (1999)[11] Beal M. J. Ph.D thesis University College London(2003)[12] T. S. Jaakola and M. I. Jordan. Stat. and Compt.
10, (2000) 25.
[13] K. Watanabe, and M. Okada, Proc. of ICONIP(2008).
[14] A. Corduneanu and C. M. Bishop. (2001) Morgan
Kaufmann[15] M. Abeles et al. PNAS. 92, (1995) 8616.
[16] I. Gat et al. Network. 8, (1997) 297.
[17] L. M. Jones et al. PNAS. 104, (2007) 18772.
[18] C. Kemere et al. J. Neurophysiol.. 100, (2008)2441.
[19] N. Achtman et al. J. Neural Eng. 4 (2007) 336.
[20] M. Churchland et al. J. Neurosci. 26(14) (2006)
3697.
[21] B. Yu et al. J. Neurophysiol. 102, (2009) 614.
[22] J. P. Donoghue. Nat. Neurosci. Supple. 5, (2002)1085.
[23] L. R. Hochberg et al. Nature 442, (2006) 164.
[24] J. Cunningham et al. Neural Netw. 19, (2009) inpress.
[25] A. P. Dempster et al. J. Roy. Statist. Soc. B 39
(1977) 1.
[26] E. N. Brown et al. In Computational Neuroscience: Comprehensive Approach (2003)[27] W. Truccolo et al. J Neurophysiol. 93 (2005) 1074.
[28] C. M. Bishop and M. Svensen Proc. on UAI(2003) 57.
[29] K. Katahira et al. J. Phys. 95 (2008) 012015.
[30] L. E. Baum et al. Annals of Math. Stat. 41 (1970)
164.
[31] S. N. Baker and R. N. Lemon. J. Neurophysiol.
84(4) (2000) 1770.
[32] A. C. Smith and E. N. Brown, Neural Computation 15 (2003) 965.
[33] H. Shimazaki and S. Shinomoto, Neural Coding
2007 Montevideo, Uruguay.
[34] M. Nawrot et al. J. Neurosci. Methods 94 (1999)
81.
[35] B. Richmond et al. J. Neurophysiol. 64(2) (1990).
[36] I. DiMatteo et al. Biometrika 88 (2001) 1055.
[37] Endres et al. Advances in NIPS 20 (2008)[38] K. Watanabe et al. IEICE trans. Info. and Syst.,(in press).
[39] J. H. Maunsell and D. C. Essen, J. Neurophysiol
49 (1983) 1127.
[40] K. H. Britten et al. The Neural Signal Archive
nsa2004.1[41] K. H. Britten et al. J. Neurosci. 12(12) (1992)
4745.
[42] W. Bair and C. Koch, Neural Compt. 8 (1996)
1185.
[43] S. G. Lisberger and J. A. Movshon j. Neurosci.
19(6) (1999) 2224.
[44] L. C. Obsorne et al. j. Neurosci. 24(13) (2004)
3210.
[45] D. Ji and W. A. Wilson, Nat. Neurosci. 10(1)(2007) 100.
[46] Z. Chen et al. Neural Compt. 21 (2009) 1797.
[47] P. J. Green, Biometrika 82(4) (1995) 711.
[48] H. Shimazaki and S. Shinomoto, Neural Comput.
19 (2007) 1503.