読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ワイブル過程のシミュレーション

R 確率過程

逆関数法による乱数の生成

ワイブル過程(Weibull process)とは,累積ハザード関数  H(t) =\left( \frac{t}{\alpha} \right)^\beta を持つ非斉次的(non homogeneous)なポアソン過程の一種.

この乱数列  t_1, t_2, \ldots, t_n は分布関数 F(t)=\exp(-H(t)) に従う.

 z_1, z_2, \ldots, z_n を一様乱数として,逆関数法により, t_1=\alpha(- \log (z_1))^{1/\beta}

 t_k は条件付き確率分布
\displaystyle P(T>t_k|T>t_{k-1})= \frac{\exp\left( -( \frac{t_k}{\alpha} )^\beta \right) } {\exp \left(- ( \frac{t_{k-1}}{\alpha})^\beta \right)}
に従うから,逆関数法により,
\displaystyle  t_k=\alpha \left(- \log (z_k) + \left(\frac{t_{k-1}}{\alpha}\right)^\beta  \right)^{1/\beta}

R のコード:

Weibull_process = function(n,alpha,beta){
  z=runif(n)
  w = alpha*(-log(z[1]))^(1/beta)
  if(n >= 2){
    for(i in 2:n){
      w[i] =alpha*( (-log(z[i]))^(1/beta) + (w[i-1]/alpha)^beta)^(1/beta)
    }    
  }
  w
}

シミュレーション

 \alpha >1 のときは事象の起こりやすさ(強度)は時間とともに増加する,

w = Weibull_process(50, 2, 3)
stripchart(w)

f:id:abrahamcow:20140615051356p:plain

 \alpha =1 のときは事象の起こりやすさは時間に依存しない. 点の集合は定常ポアソン過程になる. よって点と点の間隔は指数分布に従う.

w2 = Weibull_process(1000, 2, 1)
dx2 <- diff(w2)
hist(dx2,freq=FALSE)
curve(dexp(x,1/2),add=TRUE)

f:id:abrahamcow:20150320163125p:plain

qqplot(qexp(ppoints(length(dx2))), dx2)
qqline(dx2, distribution=function(p){qexp(p)})

f:id:abrahamcow:20150320163146p:plain

 \alpha <1 のときは事象の起こりやすさは時間とともに減少する.

w3 = Weibull_process(50, 2, 0.3)
stripchart(w3)

f:id:abrahamcow:20150320165423p:plain

ワイブル過程はワイブル分布とは別ものである. しかし, 最初の1つ目の事象のみを取り出すと, ワイブル過程から生成される乱数はワイブル分布に従う.

w = sapply(1:10000,function(i){Weibull_process(1, 2, 3)})
hist(w,breaks="scott",freq=FALSE)
curve(dweibull(x,3,2),add=TRUE)

f:id:abrahamcow:20150421232951p:plain

ポアソン分布の確率分布関数にワイブル過程の累積ハザード関数を代入して,

\displaystyle P(T>t) = e^{-H(t)}  \frac{H(t)^0}{0!} = \exp\left( -\left( \frac{t}{\alpha} \right)^\beta \right)

これはワイブル分布.