HOME PAGE

Spike Density Estimation

- Histogram Method

- Kernel Method

Share it with your friend

Tweet

スパイク密度推定

ヒストグラム密度推定    (PSTH)

電気神経生理学の動物実験では感覚刺激・行動・注意等と神経細胞の発火頻度（レート）の相関関係がよく調べられる. 広く使われているレート推定の手法に, 同一刺激下で行われた複数回の試行のスパイク時系列を適当な時間幅をもつ区間に分割し, その中でのイベント生成率（発火率）を棒グラフとして表す Peristimulus Time Histogram (PSTH)がある. PSTHの形状は分割する区間の時間幅に依存するにもかかわらず, 区間幅は多くの場合研究者により恣意的に与えられている.

ここでは平均二乗誤差最小化の観点から最適区間幅を決定する簡便な公式を導出する. またレートが一般の定常確率過程の場合の最適区間幅の理論値のスケーリング則を導出し, 転移点近傍での振る舞いを調べる. この理論の応用としてコスト関数を外挿することで, データからヒストグラム作成に必要な最小試行数を求めたり, 背後のレート過程がなめらかな過程か否かを推定することができる.

スパイク時系列からのヒストグラムの作成    

長さ$T$の時間依存Poisson過程のレート（強度過程）を $\lambda_{t}$ ( $t\in\lbrack0,\,T]$)とする. 区間 $[0,\,\Delta]$の棒ヒストグラムの真の高さは

$$\theta=\frac{1}{\Delta}\int_{0}^{\Delta}{\lambda_{t}dt}.$$ (3.1)

で与えられる. この区間内の$n$回の試行数の総スパイク数$k$は次のPoisson分布で与えられる.

$$p(k\vert n\theta\Delta)=\frac{\left( n\theta\Delta\right) ^{k}}{k!}% e^{-n\theta\Delta}%$$ (3.2)

従って ${\theta}$の不偏推定量である ${{\hat{\theta}}}=k/(n\Delta)$がデータから求められるヒストグラムの高さである.

平均積分二乗誤差及びコスト関数の導入

スパイク時系列のレート $\lambda_{t}$とヒストグラム $\hat{\lambda}_{t}$の当てはまりの良さは平均積分二乗誤差(Mean Integrated Squared Error, MISE)で評価する. 十分長い定常なスパイク時系列が与えられた場合MISEは次式で与えられる[Shimazaki and Shinomoto, 2007].

$$\equiv\frac{1}{T}\int_{0}^{T}E\,(\hat{\lambda}_{t}-\lambda _{t})^{2}\,dt$$ MISE (3.3)

$$=\left\langle {E\left[ {\frac{1}{\Delta}\int_{0}^{\Delta}{(\,{\lambda _{t}-}\hat{\theta}\,)^{2}dt}}\right] }\right\rangle$$ (3.4)

ここで $\left\langle \,\cdot\,\right\rangle$はレート過程 $\lambda_{t}$の経路によるアンサンブル平均を意味し, $E_{\Lambda }\left[ \,\cdot\,\right]$はレート $\lambda_{t}%$の区間幅$\Delta$内での時間平均${\theta}$がである場合のスパイク数の条件付き確率分布(式3.2)による平均操作を表す.

MISEを区間幅$\Delta$内でのレートのゆらぎとスパイク生成のゆらぎに分割し, さらに区間幅$\Delta$の選択に依らない項を除いたコスト関数を導入することができる.

$$C(\Delta) \equiv\mathbf{MISE}-\left\langle {\left( {\lambda_{t}% -\langle\theta\rangle}\right) ^{2}}\right\rangle$$

$$=\left\langle E\left[ (\hat{\theta}-\theta)^{2}\right] \right\ran... ...left\langle {\left( \theta{-\langle\theta\rangle}\right) ^{2}% }\right\rangle%$$ (3.5)

式3.5は$\theta$の分散を第２項に含んでいるので, 観測量のみからなる式に書き直し次式を得る.

$$C\left( \Delta\right) =2\left\langle E\left[ (\hat{\theta}-\theta... ...ngle E\left[ (\hat{\theta}% {-\langle\theta\rangle})^{2}\right] \right\rangle%$$ (3.6)

ヒストグラムの最適区間幅決定の手順を以下にまとめる[Shimazaki and Shinomoto, 2007].

(i)

$n$回の試行により得られたスパイク時系列について， その観測期間 $T$ を幅$\Delta$の $\ N$ 個のビンに区切る．$i$番目のビンに入るスパイクの数を数え， $k_{i}$ とする．

(ii)

スパイク数$\{k_{i}\}$ の平均と分散を計算する

$$\bar{k}\equiv\frac{1}{N}\sum_{i=1}^{N}k_{i}$$, and $$\displaystyle v\equiv\frac{1}{N}\sum_{i=1}^{N}(k_{i}-\bar{k})^{2}.$$

(iii)

コスト関数を計算する,

$$C_{n}(\Delta)=\frac{2\bar{k}-v}{(n\Delta)^{2}}.$$

(iv)

異なるビン幅$\Delta$に対してi から iii を繰り返し，コスト関数 $C_{n}% (\Delta)$の最小値を与える $\Delta^{\ast}$を探す．

最適区間幅の理論値のスケーリング則と発散

式3.5の第一項にCramér-Raoの不等式を適用することで, コスト関数の下限が$\Lambda$の統計量で与えられる.

$$C\left( \Delta\right) \geq\frac{\left\langle \theta\right\rangle ... ...left\langle {\left( \theta{-\langle\theta\rangle}\right) ^{2}% }\right\rangle%$$ (3.7)

右辺の極値を考えることで, レートが平均$\mu$, 相関関数$\phi(t)$なる定常確率過程について最適幅の解析解が求まる.

$n$が十分大きい場合, $\phi(t)$の原点付近での展開式を用いて式3.7右辺の極値を与える$\Delta$を求める. $\phi(t)$が原点でCusp型となるときは漸近値 $\phi^{\prime}(0+)$を用いて最適幅は $\Delta^{*} \sim \sqrt{{{ - 3\mu} \mathord{\left/ {\vphantom {{ - 3\mu } {\phi ... ...right)n}}} \right. \kern-\nulldelimiterspace} {\phi '\left( {0 + } \right)n}}}$で与えられる. $\phi(t)$が原点でなめらかなときは対称性から $\phi^{\prime}% (0)=0$であり, 最適幅は $\Delta^{*} \sim\left( {{{ - 6\mu} \mathord{\left/ {\vphantom {{ - 6\mu } {\phi... ... ^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}$となる.

$n$が小さい転移点付近では $\langle(\theta-\langle\theta\rangle)^{2}% \rangle\simeq\mu/n_{c}(1/\Delta)-u(1/\Delta)^{2}+O\left( (1/\Delta )^{3}\right)$と展開する($n_{c}$, $u$は定数). このとき臨界点は$n_{c}$でありランダウの２次相転移の理論が適用できる. $n>n_{c}$では最適幅の振る舞いは $\Delta^{\ast}\sim nn_{c}/(n-n_{c})$で表される.

ヒストグラム作成に必要な最小試行回数の推定    

$n$回の試行数のスパイク統計$\bar{k}$, $v$を用いて, 試行数が$m$回の場合の最適幅を, 外挿したコスト関数

$$C_{m}\left( \Delta\right) =\left( \frac{1}{m}+\frac{1}{n}\right) \frac{\bar{k}}{n\Delta^{2}}-\frac{v}{(n\Delta)^{2}}%$$ (3.8)

から推定することができる. これによりデータから次の２つの数値を推定できる. (i) スケーリング指数: 最適幅の試行数に関するスケーリング指数を調べることができる. 指数が$-1/2$のときは背後のレートは微分不可能であり, $-1/3$のときには微分の存在するなめらかな確率過程であると推定される. (ii) ヒストグラム作成に必要な最小試行回数: 転移点で発散やとびを示す指標を用いることで, 転移点を与える$n_{c}$を推定することができる. これにより実験者は少ない試行数からヒストグラム作成に最低限必要な試行数を予測することができる.

カーネル密度推定

神経科学の古典的な実験では，記録されたスパイク時系列は刺激・行動などの開始時刻でそろえてならべる．これをラスタープロットという．これらスパイク時系列を重ね合わせたスパイク時系列データ ${t}_{i}$ ( $i=1,2,\cdots,N$) を次のようにデルタ関数を用いて表す．

$$x_{t}=\frac{1}{n}\sum_{i=1}^{n}{{\delta\left( {t-t}_{i}\right) }}, %$$ (3.9)

ここで $n$ は繰り返し試行の数である．カーネル密度推定はこのデータ時系列$x_{t}$とカーネル関数$k(s)$ の畳み込み積分で与えられる．

$$\hat{\lambda}_{t}=\int x_{t-s}k{\left( {s}\right) }\,ds. %$$ (3.10)

以下積分 $\int$ は $\int_{-\infty}^{\infty}%$を表す. カーネル関数は密度の条件, $\int k(s)\,ds=1$, を満たし，中心がゼロで, $\int sk(s)\,ds=0$, かつ有限のバンド幅を有する, $w^{2}=\int s^{2}k(s)\,ds<\infty$.

一般に良く用いられるのはガウスカーネルであり，これは正規密度分布のこと．

$$k_{w}(s)=\frac{1}{\sqrt{2\pi}w}\exp{\left( -\frac{s^{2}}{2w^{2}}\right) }, %$$ (3.11)

ここで $w$ がバンド幅である．以下このバンド幅を最適化する公式を導く．

カーネルバンド幅最適化

最適カーネルバンド幅は以下のコスト関数を最小化するバンド幅により得ることができる(Shimazaki & Shinomoto 2010[Shimazaki and Shinomoto, 2010] )．

$$\hat{C}_{n}\left( w\right) =\frac{1}{n^{2}}\sum_{i,j}\int{k}_{w}{... ...}\right) }\,dt-\frac{2}{n^{2}}% \sum_{i\neq j}k_{w}\left( t_{i}-t_{j}\right) ,$$

多次元のカーネルの場合も同様に

$$C\left( \mathbf{w}\right) =\sum_{i,j}\int k_{\mathbf{w}}\left( \m... ...bf{x-2}\sum_{i\neq j}k_{\mathbf{w}}\left( \mathbf{x}_{i}-\mathbf{x}_{j}\right)$$

ここでは $\mathbf{x=}\left( x_{1},x_{2},\ldots ,x_{p}\right) .$ $\mathbf{x=}\left( x,y\right) \mathbf{x}_{i}% \mathbf{=}\left( x_{i},y_{i}\right)$ $i=1,\ldots,N$    . ただし, $\left\Vert \mathbf{x}\right\Vert \mathbf{=}\sum_{i}% ^{p}x_{i}^{2}$とする. 二次元の場合は $\left\Vert \mathbf{x}\right\Vert \mathbf{=}x^{2}+y^{2}$．

多次元カーネル幅最適化公式の導出概略

デルタ関数を用いて，データを次のように表すことができる．

$$s\left( \mathbf{x}\right) =\sum_{i=1}^{n}{{\delta\left( \vert\vert{x-x}% _{i}\vert\vert+\vert\vert{y-y}_{i}\vert\vert\right) }}%$$ (3.12)

カーネル密度推定は

$$\hat{\lambda}\left( x,y\right) =\int s\left( \mathbf{x}\right) k{\left( \mathbf{x}-\mathbf{u}\right) d}\mathbf{u}.$$ (3.13)

$$=\sum_{i=1}^{n}k\left( \mathbf{x}-\mathbf{x}_{i}\right)$$ (3.14)

MISE

   MISE$$=\int E[(\lambda\left( \mathbf{x}\right) -\hat{\lambda}\left( \mathbf{x}\right) )^{2}]d\mathbf{x},$$

からコスト関数を導入する

$$C_{n}\left( w\right) = -\int\lambda\left( \mathbf{x}\right) ^{2}\,d\mathbf{x}$$

$$=\int E\hat{\lambda}\left( \mathbf{x}\right) ^{2}\,d\mathbf{x}% -... ...a\left( \mathbf{x}\right) E\hat{\lambda}\left( \mathbf{x}% \right) d\mathbf{x}.$$ (3.15)

第二項は

$$Es\left( \mathbf{x}\right) E\hat{\lambda}\left( \mathbf{x}\right)... ...t{\lambda }\left( \mathbf{x}\right) -E\hat{\lambda}\left( \mathbf{x}\right) )],$$ (3.16)

さらに，上式の第二項は

$$E[(s\left( \mathbf{x}\right) -Es\left( \mathbf{x}\right) )(\hat {\lambda}\left( \mathbf{x}\right) -E\hat{\lambda}\left( \mathbf{x}\right) )]$$

$$\hskip2em=\int\int k{\left( \mathbf{x}-\mathbf{u}\right) }E\left[... ...t( s\left( \mathbf{x}\right) -Es\left( \mathbf{x}\right) \right) \right] \,dudv$$

$$\hskip2em=\int\int k{\left( \mathbf{x}-\mathbf{u}\right) }\left[ ... ...hbf{x}-\mathbf{u}\right) \frac{{1}}{n}Es\left( \mathbf{x}\right) \right] \,dudv$$

$$\hskip2em=\frac{{1}}{n}{k}_{w}(\mathbf{0})Es\left( \mathbf{x}\right) .$$ (3.17)

これらをコスト関数に代入して

$$C_{n}\left( w\right) =\int E\hat{\lambda}\left( \mathbf{x}\right) ^{2}\,d\mathbf{x}-2\... ...bda\left( \mathbf{x}\right) E\hat{\lambda}\left( \mathbf{x}\right) d\mathbf{x}.$$

$$=\int E\hat{\lambda}\left( \mathbf{x}\right) ^{2}\,d\mathbf{x}% -... ...ac{{1}}{n}{k}_{w}(\mathbf{0})Es\left( \mathbf{x}% \right) \right] d\mathbf{x}%$$

サンプルからの推定は以下の式で与えられる．

$$\hat{C}_{n}\left( w\right) =\int\hat{\lambda}\left( \mathbf{x}\right) ^{2}\,d\mathbf{x}-2\in... ... -\frac{{1}}{n}{k}_{w}(\mathbf{0})s\left( \mathbf{x}\right) \right] d\mathbf{x}$$

$$=\int\hat{\lambda}\left( \mathbf{x}\right) ^{2}\,d\mathbf{x-2}\le... ...lambda}\left( \mathbf{x}_{i}\right) -\frac{{1}}{n}{k}% _{w}(\mathbf{0})N\right]$$

$$=\sum_{i,j}\int k\left( \mathbf{x}-\mathbf{x}_{i}\right) k\left( ... ...\right) d\mathbf{x-2}\sum_{i\neq j}k\left( \mathbf{x}_{i}-\mathbf{x}_{j}\right)$$

対称なカーネルのコスト関数はさらに簡略化できる．

$$\hat{C}_{n}\left( w\right) =\sum_{i,j}\psi\left( \mathbf{x}% _{i},\mathbf{x}_{j}\right) \mathbf{-2}\sum_{i\neq j}k\left( \mathbf{x}% _{i}-\mathbf{x}_{j}\right)$$

$$=\sum_{i}\psi\left( \mathbf{x}_{i},\mathbf{x}_{i}\right) +2\sum_{... ...,\mathbf{x}_{j}\right) -4\sum_{i<j}k\left( \mathbf{x}_{i}-\mathbf{x}_{j}\right)$$

$$=\sum_{i}\psi\left( \mathbf{x}_{i},\mathbf{x}_{i}\right) +2\sum _... ...i},\mathbf{x}_{j}\right) -2k\left( \mathbf{x}_{i}-\mathbf{x}_{j}\right) \right]$$

            

ガウスカーネルの場合，これを用いて以下の公式を得ることができる．２次元のガウスカーネルは

$$k_{w}(\mathbf{x})=\frac{1}{2\pi w^{2}}\exp{\left( -\frac{\vert\vert\mathbf{x\vert\vert}% }{2w^{2}}\right) },$$ (3.18)

ゆえに，

$$k_{w}\left( \mathbf{x}_{i}-\mathbf{x}_{j}\right) =\frac{1}{2\pi w... ...t[ -\frac{\left\Vert \mathbf{x}_{i}-\mathbf{x}_{j}\right\Vert }{2w^{2}}\right]$$

及び

$$\int k_{w}\left( \mathbf{x}-\mathbf{x}_{i}\right) k_{w}\left( \ma... ...t[ -\frac{\left\Vert \mathbf{x}_{i}-\mathbf{x}_{j}\right\Vert }{4w^{2}}\right]$$

を用いて以下のコスト関数が導出できる．

$$\hat{C}_{n}\left( w\right) =\frac{N}{4\pi w^{2}}+\frac{1}{2\pi w^... ...c{\left\Vert \mathbf{x}_{i}-\mathbf{x}_{j}\right\Vert }{2w^{2}}\right] \right]$$