HOME PAGE

Spike Density Estimation

- Histogram Method

- Kernel Method

Share it with your friend

Tweet

非定常ポアソン過程

時間伸縮

Figure 2.1: 変動スパイクレートと時間伸縮

本節では, 定常ポアソン過程の拡張として, 図-aのように時間的に変動するレートから生成される非定常ポアソン過程を考えたい. レートの低いところではスパイクの発生率は低く, レートの高いところではスパイクの発生率は高くなるだろう. 変動レートを $\lambda_{t}$で表し, 変動レートを時刻$t$までを積分した

$$\Lambda\left( t\right) =\int_{0}^{t}\lambda\left( u\right) du %$$ (2.1)

を考える. 図-bに関数 $\Lambda\left( t\right)$を示す. $\Lambda$は無次元量で, いわば $\lambda_{t}$により規格化された時間である. 変動レートの小さなところでは, 実時間に対して$\Lambda$はゆっくりと進み, 変動レートの大きな所では$\Lambda$は早く進む. そのため, 2.1の変数変換を行うことを時間伸縮(Time-rescaling）をすると言う. 例えば図-bでは, リスケールされた時間軸$\Lambda$上にある２点 $\Lambda _{1},\Lambda_{2}$の間隔と２点 $\Lambda _{3},\Lambda_{4}$の間隔はほぼ同じぐらいだが, 実時間上ではレートの $t_{1},t_{2}$の間隔は $t_{3},t_{4}$の間隔よりずっと大きい.

リスケールされた時間軸$\Lambda$上で平均$1$の標準指数分布に従う定常ポアソン過程を考え, 実時間へ写像したスパイク時系列が非定常ポアソン過程である. $\Lambda$軸上のイベント系列を $\left\{ \Lambda_{1},\Lambda_{2},\cdots\right\}$とすると, イベント間隔 $\Lambda_{i}-\Lambda_{i-1}$が標準指数分布

$$g\left( z\right) =\exp\left[ -z\right]$$

に従う. 実時間との対応関係は $\Lambda_{i}=\Lambda\left( t_{i}\right)$だから

$$z\left( t_{i}\vert t_{i-1},\lambda_{t_{i-1}:t_{i}}\right) \equiv\... ...Lambda\left( t_{i-1}\right) =\int_{t_{i-1}}^{t_{i}}% \lambda\left( u\right) du$$

非定常ポアソン時系列を数値的に実現するためには, 指数分布に従う乱数を発生させ, これを満たす$t_{i}$を逐次求めていけばよい. この方法は時間伸縮(Time-rescaling）理論に基づく方法であり, 次のようにまとめられる[Brown et al., 2001].

アルゴリズム 8   時間伸縮理論を用いた非定常ポアソン過程の実現

1. $i=1,t_{1}=0.$

2. レート$1$の指数分布に従う乱数$\eta$を生成する.

3. $\eta=\int_{t_{i-1}}^{t_{i}}\lambda\left( u\right) du$を満たす$t_{i}$を求める.

4. $t_{i}\geq T$ならば終了.

   $t_{i}<T$ならば $i\leftarrow i+1.\ 2$へ.

アルゴリズムの第３ステップで指数分布に従う乱数を発生させるには次に説明する逆関数法を用いる.

条件付きスパイク密度分布

変動ポアソン過程において時刻$t_{i-1}$でスパイクが生じた下で, 次のスパイク時刻$t_{i}$の条件付スパイク密度分布（ISI分布）を求めてみよう：

$$p\left( t_{i}\vert t_{i-1},\lambda_{t_{i-1}:t_{i}}\right) =\left\vert \frac{dz}{dt_{i}}\right\vert g\left( z\vert\lambda_{t_{i-1}:t_{i}}\right)$$

$$=\lambda\left( t_{i}\right) \exp\left[ -\int_{t_{i-1}}^{t_{i}}% \lambda\left( u\right) du\right]%$$ (2.2)

条件つきスパイク密度分布2.2は, 定常ポアソン過程のスパイク密度分布である指数分布 $\lambda e^{-\lambda t}$の自然な拡張になっていることがわかる. 逆関数法を用いれば, この条件付スパイク密度分布に従う時系列を作れる.

アルゴリズム 9   逆関数法を用いた非定常ポアソン過程の実現

1. $i=1,t_{1}=0$.

2. $\left[ 0,1\right]$の区間の一様乱数$\xi$を生成する.

3. $\xi=\lambda\left( t_{i}\right) \exp\left[ -\int_{t_{i-1}}^{t_{i}% }\lambda\left( u\right) du\right]$を満たす$t_{i}%$を求める.

4. $t_{i}\geq T$ならば終了.

   $t_{i}<T$ならば $i\leftarrow i+1.\ 2$へ.

   
   

条件付きスパイク密度分布2.2は. $t_{i}=0,t_{i-1}=t$と見れば定常リニューアル過程におけるスパイク密度分布の定義と同義である. 定常リニューアル過程のスパイク密度分布とハザード関数との間には, 1.16のような関係があるから, 2.2とを見比べて

$$r\left( t\right) =\lambda\left( t\right)%$$ (2.3)

であることがわかる. ハザード関数も変動レートも, 時刻$t$における瞬間スパイク生成率を与え, 両者は同一視できる. ただしリニューアル過程の場合には, スパイク発生から経過した時間を変数とする特定のハザード関数を考え, スパイクの度に時刻が原点にリセットされた. これに対して, 非定常ポアソン過程の場合, 時刻$t$の瞬間発火率は過去のスパイク発火履歴に拠らない $\lambda\left( t\right)$のみできまる. 瞬間スパイク生成率がスパイク履歴に拠らないのは, 定常・非定常を問わずポアソン過程の特徴である.

直説法

次に瞬間スパイク生成率（ハザード関数） $\lambda\left( t\right)$を用いる方法が考えれる. 時間を十分小さな幅$\Delta$の区間に区切って, スパイクの発生をベルヌーイ過程で近似する. 各区間で一様乱数を生成し, スパイクの $\lambda\left( t\right) \Delta$確率で発生させる. スパイクが生じない確率は $1-\lambda\left( t\right) \Delta$である. 非定常ポアソン過程は履歴を考慮せず微小区間毎に独立に計算できる. 十分小さな$\Delta$を用いればポアソン過程・リニューアル過程を近似することができるはずである.

アルゴリズム 10   強度関数を用いた直接的な非定常ポアソン過程の実現

1. 長さ$T$の時系列を十分に小さな幅$\Delta$の $N=T/\Delta$個の区間に分ける.

2. $i=1$, $j=1$.

3. $j\leftarrow j+1$, $\quad j=N$で終了.

4. $\left[ 0,1\right]$の区間の一様乱数$\xi$を生成する.

   $\qquad\lambda\left( j\Delta\right) \Delta\geq\xi$    ならば     $t_{i}=j\Delta,i\leftarrow i+1$.

   $\qquad\lambda\left( j\Delta\right) \Delta<\xi$    ならば    $3$へ.

この直接的な方法はもっとも実装が簡潔だが, 乱数を$\ N$回発生させる必要があり計算効率は著しく悪い. 毎回乱数を生成させるよりは, 指数分布に従う乱数を一回発生させ, その値になるまで $\lambda\left( j\Delta\right)$を数値的に足し込んでいく（積分する）アルゴリズムの方が効率がよい. ただし, この方法も微小区間に区切った関数の数値積分が必要になってくる. それでは次に, 簡潔でかつ効率のよい数値計算法としてよく用いられる希薄化による手法を紹介しよう.

希薄化による数値計算法

希薄化(Thinning[Ogata, 1981,Daley and Vere-Jones, 1988,Heyman and Sobel, 1990])と呼ばれる方法は積分操作を必要としないため, 計算時間が短くて済む場合がある.

まず$[0,T]$の範囲にレートが$M$の定常ポアソン過程を生成する. ただし, $M$はどの時刻の変動レートよりも大きい値とする（つまり $t\in\lbrack0,T]$に対して $\lambda\left( t\right) \leq M$）. 生成されたスパイク時系列を $\left\{ t_{1,}t_{2}\right\}$. 時刻$t_{j}$におけるスパイクは, 規格化した変動レート $p\left( t_{j}\right) =\lambda \left( t_{j}\right) /M$の確率で残す. 残ったスパイク時系列 $\left\{ t_{1}^{\ast},t_{2}^{\ast},\cdots\right\}$は変動レート $\lambda\left( t\right)$のポアソン過程に従う.

アルゴリズム 11   Shedler-Lewis Thinning Algorithmによる非定常ポアソン過程の実現

1. $\left[ 0,T\right]$の区間の変動レート $\lambda\left( t\right)$を用意する

2. $M=\max\lambda\left( t\right)$をもとめる

3. $i=1$

4. レート$M$の指数分布に従う乱数$\xi$を生成する

5. $s=s+X$, $\quad s>T$ならば終了

6. $\left[ 0,1\right]$の区間の一様乱数$\xi$を生成する

   $\qquad\lambda\left( s\right) /M\geq\xi$    ならば $\quad t_{i}=s,$     $i\leftarrow i+1$

7. $4$に戻る