はじめに

こちらの本を読んでいます。

点過程の時系列解析 (統計学One Point)

作者:崇宏, 近江,俊一, 野村
発売日: 2019/06/08
メディア: 単行本

毎度のこと内容を簡単にまとめます。

はじめに
点過程について
ポアソン過程
- 定常ポアソン過程
  - Rでポアソン過程
- 非定常ポアソン過程
更新過程
生存時間分析との関係
参考

点過程について

点過程とは、空間上・時系列上のランダムに分布する点（イベント）の集合に関する確率過程である。

購買間隔やイベントの発生タイミングについて検討することになる。

例

商品の購買の次回列
Twitter上ので呟きの時系列
脳神経細胞のスパイク発火
余震を含む地震の発生

各要素の定義

ある観測期間 $[0,T$ ]にn回のイベントが起きた時、各イベントの発生時刻を $t_i$ とし、その集合を $t_n= \left \{ t_1, t_2, \dots, t_n \right \}$ と表現する。

各イベントの発生時刻には順序がある。

ある区間 $[0,T$ ]に発生するイベント数 $N(x,y)$ で表す。

地震の発生の時系列のように、イベントの発生に何らかの情報(規模, 発生場所, etc.)が付随する場合がある。

このような時系列は、マーク付き点過程(marked point process)と呼ばれている。(今回は扱わない)

確率密度関数

強度関数

瞬間のイベントの発生のしやすさと強度と呼び、強度が時間に応じて変化する場合に各時刻の強度の変化を強度関数で定義する。

また、過去のイベントの発生に強度が依存する場合は、条件付き強度関数と呼ぶ。

イベントの配列 $t_n$ の発生は試行ごとに異なるため、あるイベント配列 $t_n$ の起こりやすさを確率密度関数 $P_{ \left[ 0, T \right] } ( t_n )$ で表現する。

イベント数の分布

点過程では、ある観測期間に得られるイベント数 $n$ も確率に得られる変数である。

そのため、複数の観測データが与えられている場合は、イベント数の分布はデータから調べることができる。

イベント間隔の分布

あるイベントの後の次にイベントが発生する時間をイベント間間隔 $\tau_i = t_{i+1}-t_i$ と呼ぶ。

イベント間隔は観測データから調べることができるため、点過程のデータを特徴付けるのに重要である。

また、イベントの発生に関係なく任意の時刻から次のイベントが発生するまでの時間を待ち時間と呼ぶ。

ポアソン過程

ポアソン過程とは、観測期間 $[ 0,T$ ]に発生するイベント数 $n$ がポアソン分布に従う過程のことを指す。

強度関数 $\lambda(t)$ のポアソン過程では、微小な幅 $\Delta$ の買う間にイベントが発生する場合と発生しない場合の確率は次のようになる。

$P[N(t,t+\Delta)=1 ] = \lambda(t) \Delta$

$P[N(t,t+\Delta)=0 ] = 1-\lambda(t) \Delta$

ある時刻でのイベントの発生確率は他のイベントには依存しない。

また、微小な $\Delta$ において、次のような近似が成り立つ。

$1-\lambda(t) \Delta \approx exp[-\lambda(t) \Delta$ ]

そしてポアソン過程には、定常ポアソン過程と非定常ポアソン過程が存在する。

定常ポアソン過程

定常ポアソン過程では強度が時間によって変化せず $\lambda(t)=\lambda$ となり、発生確率の関数は次のようになる。

$P[N(t,t+\Delta)=1 ] = \lambda \Delta$

$P[N(t,t+\Delta)=0 ] = 1-\lambda \Delta \approx 1- exp[-\lambda \Delta$ ]

確率密度関数は、次のようになる。

$\displaystyle p_{ [ 0,T ] } (t_n) = \lambda^n exp[-\lambda T$ ]

期間 $[ 0,T$ ]イベント数の分布は、平均が $\lambda(t)$ 、分散 $\lambda (t) T$ がポアソン分布に従うことになる。

そして、イベント間間隔 $\tau$ の確率密度関数 $p(\tau)$ は、期待値が $1/\lambda$ 、分散が $1 / \lambda ^2$ の指数分布にしたがう。

Rでポアソン過程

Rでポアソン過程を扱うNHPoisson パッケージがあります。

一日1440分で $\lambda=1/60$ とした時のポアソン過程のデータを生成した。

# 定常ポアソン過程の生成
lambda = 1/60 #強度
time_span = 60*24 #1日:1440分

ps_sim <- simNHP.fun(rep(lambda,time_span))

得られたイベントタイミングを確認してみました。

ppplot(ps_sim, time_span) # 自作関数

午前中の方が少しイベント発生件数が多いように見えます。

定常ポアソン過程の強度を推定するには、fitPP.fun関数を利用します。

# 強度の推定
# b0=0 が推定する強度の初期値
pp_res <- fitPP.fun(posE = ps_sim$posNH, n = time_span, start = list(b0=0))

推定された強度の確認をします。λ=exp(β)であるため自然対数をとる必要があります。

> pp_res@coef %>% exp
        b0 
0.01388889

概ね設定通りの強度を推定できていると考えられます。

非定常ポアソン過程

非定常ポアソン過程では、強度関数 $\lambda(t)$ に従い時間によって変化する。

次のようなものが、非定常ポアソン過程を考えたものです。

確率密度関数 $P_{ \left[ 0, T \right] } ( t_n )$ は次のようになる。

$\displaystyle p_{[0,T]} (t_n) = \prod_{i=1}^n \lambda(t_i) exp \left[- \int_0^T \lambda (t) dt \right]$

観測期間に発生するイベント数は、以下の時間変換定理を適用した場合、定常ポアソン過程の場合と等しくなる。そんため、イベント数は期待値 $\Lambda (T) = \int^T_0 \lambda (t) dt$ のポアソン分布に従うことになる。

イベント間間隔 $\tau_i=t_{i+1}-t_i$ の確率密度関数 $p(\tau|t_i)$ は、時間変換定理を用いると、次のように互いに独立に期待値１の指数分布に従うことになる。

$\displaystyle \tau_i'=\Lambda (t_{t+i})-\Lambda (t_i) = \int_{t_i}^{t_i+1}\lambda (s)ds$

また、時間変換後のイベント間間隔はものとイベント間間隔 $\tau_i$ を用いると次にように表される。

$\displaystyle \tau_i'= \int_{t_i}^{t_i+\tau_i}\lambda (s)ds$

そして、時間変換定理の逆変換を考えれば、イベント間間隔 $\tau_i$ の確率密度関数[p(\tau_i|t_i)]は次のように得られる。

$\displaystyle p(\tau_i|t_i)= \lambda (t_i+\tau_i)exp \left[ -\int_{t_i}^{t_i+\tau_i} \lambda (s)ds \right]$

時間変換定理

非定常ポアソン過程は、イベントの発生間隔に合わせて時間変換を行うと、定常ポアソン過程になる。

次の時間変換定理が適用できる時に定常ポアソン過程に変換することができる。

観測期間 $[ 0,T$ ] においてイベントが強度関数 $\lambda(t)$ に従っているとする時、

次の変換式によって得られるイベントは、観測期間 $[ 0,\Lambda (T)$ ] において強度1 の定常ポアソン過程に従う。

$\displaystyle t' = \Lambda (t) = \int_0^t \lambda (s) ds$

このことから、ある過程を時間変換することにより定常ポアソン過程に従う場合、その過程は非定常ポアソン過程であると言える。

Rで非定常ポアソン過程

時変の強度関数を作ります。

all.seconds　<-　seq(1, time_span, length.out=time_span)
lambdas <- 0.05 * exp(-0.005 * all.seconds) + 0.005

可視化すると次のようになります。深夜から明け方にかけてイベント発生強度が強いようなイメージです。

f:id:saltcooky:20201227022239p:plain

シミュレーションすると次のような結果が得られました。

ps_sim <- simNHP.fun(lambdas)

ppplot(ps_sim)

f:id:saltcooky:20201227020734p:plain

強度関数を推定します。

# 推定
npp_res <- fitPP.fun(tind = TRUE,
                 covariates = cbind(all.seconds),
                 posE = ps_sim$posNH,
                 start = list(b0=0, b1=0),
　　　　  modSim=TRUE)

推定された強度関数は次のようになりました。このグラフはfitPP.fun関数を実行すると自動的に出力されます。

概ね設定通りの強度関数を推定できました。

f:id:saltcooky:20201227021816p:plain

条件付き強度関数

条件付き強度関数とは、過去のイベント履歴 $H_t=t_i$ が与えられた時のイベント発生の強度の時系列変化を表現する関数である。

条件付き強度関数を用いた場合の点過程での、微小区間でのイベント発生確率は次のように定義できる。

$P[N(t,t+\Delta)=1 ] = \lambda (t|H_t) \Delta$

$\displaystyle P[N(t,t+\Delta)=0 ] = 1-\lambda (t|H_t) \Delta \approx 1- exp[-\lambda (t|H_t) \Delta$ ]

また、確率密度間数は次のように定義できる。

$\displaystyle p_{[0,T]} (t_n) = \prod_{i=1}^n \lambda(t_i|H_{t_i}) exp \left[- \int_0^T \lambda (s|H_s) ds \right]$

NHPoisson packageでは、ポアソン過程の強度の変化に説明変数を回帰させるモデルに拡張することも可能である。

更新過程

定常ポアソン過程ではイベント間間隔がi.i.dの指数分布に従う性質を持っていたが、更新過程ではこれを一般化したものです。

具体的には、イベント間間隔 $\tau_i=t_{i+1}-t_1>0$ が独立に確率密度関数 $f$ と累積分布関数 $F$ を持つ同一の連続確率分布に従う点過程のことを指している。

この時用いる確率分布としては、次のような分布が代表的である。

指数分布
ガンマ分布
ワイブル分布
逆ガウス分布
対数正規分布

定常更新過程

イベント発生時刻の集合に対する同時分布の確率密度関数 $p(t_n)$ は次のようになる。

$\displaystyle p(t_n) = f(t_1) \prod_{i=1}^{n-1} f(t_{i+1}-t_i)$

つまり、イベント発生時刻の同時分布の確率密度関数は、イベント間間隔に対する確率密度間数の積で表現される。

次にある観測期間 $[ 0,T$ ]における確率密度間数は、次のように最後のイベント発生時刻 $t=t_n$ から観測終了時刻 $t=T$ までイベントが起こらない確率を掛け合わせたもになる。

$\displaystyle p_{[ 0,T ] } (t_n) = f(t_1) \prod_{i=1}^{n-1} f(t_{i+1}-t_i) \{ 1-F(T-t_n)\}$

また、i番目のイベント発生時刻における確率密度関数 $p(t_i)$ は次のように表現できる。

$\displaystyle f^{k*}(t_k) = \int^\infty_0 f^{(k-1)*}(t_{k-1}) f(t_k - t_{k-1})dt_{k-1}$

これは各イベントの発生時刻に関する確率密度関数を再帰的に得ることができることを意味する。

また、周辺分布の確率密度関数から、k番目のイベント発生時刻 $t_k$ の累積分布関数は次のように得られる。

$\displaystyle F^{k*}(t_k) = \int^{t_k}_0 f^{k*}(t) dt$

観測期間内のイベント数 $N(0,T)$ がk個以下になることと、 $t+1$ 回目のイベント発生時刻 $t_{k+1}$ が観測終了時刻 $t=T$ より後に来ることは同値な事象であることから、次の関係式が得られる。

$\displaystyle p[N(0,T) \leq k ] = p(t_{k+1}>T)=1-F^{(k+1)*} (T)$

これを用いると特定のイベント回数 $k$ となる時の確率関数は次のようになる。

$\displaystyle p(k) = p [N(0,T) \leq k ] - p [N(0,T) \leq k-1 ]=F^{k*} (T) -F^{(k+1)*} (T)$

期間内に発生するイベントの期待値は、次にようになります。

$\displaystyle E[N(0,T) ]=\sum_{k=1}^\infty F^{k*} (T)$

そして、イベント間間隔の期待値と次の関係が成り立つ。

$\displaystyle \lim_{T \to \infty } \frac{E[N(0,T)]}{T} = \frac{1}{\mu}$

非定常更新過程

非定常更新過程では、ポアソン過程における非定常過程と同様に定常更新過程から時間変換によって構築できる。

そのため、ある過程が以下の時間変換を行うことで、平均強度関数 $1$ に従っている時にイベント $t_n$ は平均強度関数 $v(t)$ の非定常更新過程であると言える。

$\displaystyle t'=\Lambda (t) = \int_0^t v(s)ds$

ここで、変換後のイベントは $t'_n = \left \{ \Lambda (t_1), \Lambda (t_2), \dots , \Lambda (t_n) \right \}$ である。

期間 $[0,T$ ]の平均強度関数 $v(t)$ 非定常更新過程の確率密度関数は、次のように定義できる。

$\displaystyle p_{[ 0,T ] } (t_n) = \prod_{i=1}^n v(t_i) × \{1-F(\Lambda (t_1))\} \prod_{i=1}^{n-1} f(\Lambda (t_{i+1}- \Lambda (t_i))) \{ 1- F( \Lambda (T)-\Lambda (t_n) )\}$