非同次/確率的な強度λを持つ指数分布に従う確率変数τ(死亡/消滅時刻)の生成

非同次な指数分布

非同次なポアソン分布ってのがある。日本語で定義がちゃんと載ってるページがなかったので、英語版のwikipediaを参照すると

Poisson point process - Wikipedia

が該当する。要するにこれは、今まで、一定だと思っていたポアソン分布の強度(普通のポアソン分布だと平均だと思っていたパラメータ) $\lambda$ を時間や空間に依存させてみたりするってもんだ。あと、この拡張として、この $\lambda$ 自身を確率過程だと思っても以下の話は大体同じになる。正確な測度論の記述じゃないけど。

そして、

指数分布とポアソン分布のいけない関係 from tera monagi

の資料に説明があるように、その非同次なポアソン分布の裏側には、非同次な指数分布が生まれてるんだろうなぁと考えることが出来るわけで、ここではそれを考えたい。具体的には、非同次な指数分布に従う乱数を生成するにはどうしたらよいのかという話だ。日本語の教科書にこの辺の記載があるのかないのかもよくわからないが、金融だと、

の「3.5 Process with Jumps」や「9.4 Credit Risk」を見ておくと良い。日本語だとたぶん信頼性工学系の書籍を漁ることになるんだとおもう、多分。

非同次な指数分布の数式展開

さて、(非)同次な指数分布に従う確率変数を $\tau$ と書くことにすると、まず、同次な指数分布の分布関数(時刻 $T$ までに死亡/消滅していない確率)は
$P(\tau < T) = 1-\exp\left\{-\lambda T \right\}$
と書くことができるが、今、非同次なケースを考えようとすると、強度 $\lambda$ が変数(時間・空間）依存性を持つので、以下のように積分形式で書くことができる。
$P(\tau < T) = 1-\exp\left\{-\int_0^T \lambda(t)dt \right\}$
これは $\lambda$ が一定値だと思うと、ちょうど同次のケースに戻るので、これで定義としては良さそうだぞと、逆に、これを持って非同次な指数分布の分布関数とするぞとそういうことです。

ここで、ここから先、逆変換法を使っていくが、逆変換法については

Rで学ぶ逆変換（逆関数）法 from tera monagi

あたりを見てください*1。

さて、非同次な指数分布に対して逆変換法を適用するために区間[0, 1]の一様乱数 $u$ を導入し、
$u = 1-\exp\left\{-\int_0^T \lambda(t)dt \right\}$
と書く。あとは、これをTについて解いてやればよいので、ちょっと変形して
$\log(1-u) = \int_0^T \lambda(t)dt$
とする。ここから先は解析的に解くのがキツいので、こいつはこのままで放置するとして、この数式を満たすような最小の $T$ が非同次な指数分布に従うというわけだ。もう少しちゃんと書くと
$\tau = \inf\left\{ T>0; \int_0^T \lambda(t)dt = \log(1-u) \right\}$
という形になる。更に、ここで、
$\log(1-u)$
は逆変換法の視点から解釈するに、 $\lambda=1$ とした指数分布になっているので、それを $\xi$ と書くことにすると、
$\tau = \inf\left\{ T>0; \int_0^T \lambda(t)dt = \xi \right\}$
となる。なので、シミュレーションの手順としては

1: $\lambda=1$ (平均1)とした指数分布から一個乱数 $\xi$ を取得する

2: $\int_0^T \lambda(t)dt$ の計算を数値積分なり解析解で頑張って計算する

3: $\int_0^T \lambda(t)dt = \xi$ を満たすような $T$ を求めて、こいつを非同次な指数分布からの乱数だと思う
という手順になる。

Rでやってみる

というわけで、上で書いたアルゴリズムを実際に書いてみた。ここでは、
$\lambda(t) = 15t^2 + 1$
と設定している。深い意味はなく、解析的に計算したかったのでこうしている。この関数の積分は解析的に計算できるので、全体として確率密度関数 $f(T)$ は
$f(T)dT = P(T < \tau < T+dt) = P(\tau < T+dt) - P(\tau < T) \sim \lambda(T) \exp\left\{ -\int_0^T \lambda(t)dt \right\} = (15T^2 + 1)\exp\left\{ -(5T^3+T) \right\} dT$
となるんで、これとサンプリング結果を比較する。結果をPLOTすると、以下のようにほぼ一致しているので、これで良さそうだ。

ソースは以下。コメントにあるように平均1の指数分布と強度 $\lambda$ の積分が一致する点はめんどくさいのでソルバー(uniroot)関数に投げちゃってる。あと、乱数生成は汎用的にしたかったので数値積分で処理しているので遅いぞと。

#int_0^t \lambda(u)du - \xi
#な関数。無理やり根を探す形で書くのでこうした
objective <- function(t, lambda, xi){
  integrate(lambda, 0, t)$value - xi
}
#適当な強度の関数形
lambda <- function(t){5*3*t^2+1}
#描画範囲
range_plot <- c(0, 2)
#サンプリング
set.seed(100)
x <- sapply(
  1:10^4, 
  function(i){
    uniroot(objective, range_plot, tol=10^(-3), lambda=lambda, xi=rexp(1))$root
  }
)
#解析解との比較
#上のラムダを積分（手でがんばる）
lambda_integrated <- function(t){5*t^3+t}
#解析解（確率密度関数）
analytical_solution <- function(t){
  lambda(t)*exp(-lambda_integrated(t))
}
#比較PLOT
library(ggplot2)
qplot(x, geom="blank") + 
  geom_histogram(aes(y=..density..), fill="aquamarine",color="black") + 
  stat_function(fun=analytical_solution, color="red", size=1) + 
  scale_x_continuous(limits=range_plot)