rで化合物ポアソン過程をシミュレートする

Question

私はrで複合ポアソン過程をシミュレートしようとしています。プロセスは$ \ sum_ {j = 1}^{N_t} Y_j $で定義されます。ここで、$ Y_n $はi.i.dシーケンスに依存しません。$ N（0,1）$ valuesと$ N_t $は$ 1 $パラメータのポアソンプロセスです。私は運のないrでこれをシミュレートしようとしています。 Simutale 0からTまでのCPP：

が開始：私は次のようにこれを計算するためのアルゴリズムを持っている$ K = 0 $

繰り返しをしながら$ \ sum_ {i = 1}^K T_I < Tの$

セット$ K = K + 1 $

シミュレー$ T_K \シムのEXP（\ラムダ）（= 1 $私の場合の$ \ラムダで）$

シミュレー$ Y_K \ simのN（0 、1）$（これは単なる特殊なケースですが、これを任意のディストリビューションに変更できます）軌道は次のように与えられる。ここで、軌道は次のように与えられる。ここで、軌道は次のように与えられる。ここで、$ N（t）= sup（k：\ sum_ {i = 1}^k T_i \ leq t） $

私はこのプロセスをプロットできるように、誰かが私をrでシミュレートする手助けをすることができますか？私は試しましたが、それをやり遂げることはできません。

Answer 1

N_tとX_tを決定する累積合計にはcumsumを使用してください。この例示的なコードは、nをシミュレートする回数を指定し、時間をn.tにシミュレートし、値をxにして（それが行ったことを表示するために）軌跡をプロットします。

n <- 1e2 
n.t <- cumsum(rexp(n)) 
x <- c(0,cumsum(rnorm(n))) 
plot(stepfun(n.t, x), xlab="t", ylab="X") 

このアルゴリズム、それは低レベルの最適化機能に依存するためには、高速である：私はそれをテストした6歳のシステムは、百万人以上（時間、値）のペアを生成します。毎秒。シミュレーションのために十分に通常は良いのですが、それは非常に私たちは、上記のコードを活用することができ、時間Tに出シミュレーションを生成するように求める問題は、満足していませんが、解決策は少しトリッキーです

---

。それは、時間Tの前にポアソン過程において何回起こるかについて合理的な上限を計算する。それは到着間時間を生成する。これは、時間Tに実際に到達していない（まれな）イベントで手順を繰り返すループでラップされます。

追加の複雑さは、漸近計算時間を変更しません。

T <- 1e2   # Specify the end time 
T.max <- 0   # Last time encountered 
n.t <- numeric(0) # Inter-arrival times 
while (T.max < T) { 
    # 
    # Estimate how many random values to generate before exceeding T. 
    # 
    T.remaining <- T - T.max 
    n <- ceiling(T.remaining + 3*sqrt(T.remaining)) 
    # 
    # Continue the Poisson process. 
    # 
    n.new <- rexp(n) 
    n.t <- c(n.t, n.new) 
    T.max <- T.max + sum(n.new) 
} 
# 
# Sum the inter-arrival times and cut them off after time T. 
# 
n.t <- cumsum(n.t) 
n.t <- n.t[n.t <= T] 
# 
# Generate the iid random values and accumulate their sums. 
# 
x <- c(0,cumsum(rnorm(length(n.t)))) 
# 
# Display the result. 
# 
plot(stepfun(n.t, x), xlab="t", ylab="X", sub=paste("n =", length(n.t)))