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

廿TT

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

強度関数(intensity function)が区分的に定数(piecewise constant)なポアソン過程のシミュレーション

確率過程 R

※コードが思いっきりまちがっていたので訂正しました.
lambda[i]* (t[i+1] - t[i]) としなければいけないところを lambda[i] としていました.
lambda は1単位長さ辺りに生成する乱数の個数に相当するので, (t[i+1] - t[i]) を掛けなきゃいけなかったんでした.すみません.(6月21日)


強度関数(intensity function)が区分的に定数(piecewise constant)なときの,すなわち,
f:id:abrahamcow:20140614032652p:plain
 t_0 < t_1 < \cdots < t_n
ポアソン過程.

R で乱数生成.

引数 t: t_i,lambda:  \lambda_i

#関数を定義
PoissonProcessStep <- function(t, lambda) {
  stopifnot( length(t) == length(lambda)+1)
  X <- c()
  N <- c()
  for (i in 1:length(lambda)) {
   N[i] <- rpois(1, lambda[i]* (t[i+1] - t[i]) ) #ポアソン乱数
   X <- c(X, runif(N[i], t[i], t[i+1])) #一様乱数
  }
  N <- c(N,0)
  list(rv=X,count=cbind(t,N))
}
#
r1 =PoissonProcessStep(c(1:5),c(20,100,50,70))
#グラフ
plot(r1$count[,1],r1$count[,2], type="s",
     ylab="frequency",xlab="t")
tmp <- cbind(r1$rv,0)
points(tmp, col=densCols(tmp),pch=16)

f:id:abrahamcow:20140614033850p:plain
これは,
\lambda=(20,100,50,70)',
 t_0=1,t_1=2,\ldots,t_5=5n=5とした例.


関連エントリ:
一様乱数の点と点の間隔は指数分布に従う。 - 廿TT


参考文献:

An Introduction to Statistical Computing: A Simulation-based Approach (Wiley Series in Computational Statistics)

An Introduction to Statistical Computing: A Simulation-based Approach (Wiley Series in Computational Statistics)