3

Gillespie SSA を使用して、感染 (寄生虫) の確率モデルを作成しました。モデルは「GillespieSSA」パッケージ ( https://cran.r-project.org/web/packages/GillespieSSA/index.html ) を使用しました。要するに、コードは個別のコンパートメントの集団をモデル化します。コンパートメント間の移動は、ユーザー定義のレート方程式に依存します。SSA アルゴリズムは、特定のタイム ステップ (タウ) の各レート方程式によって生成されるイベントの数を計算し、それに応じて母集団を更新します。プロセスは特定の時点まで繰り返されます。問題は、イベントの数がポアソン分布 (Poisson(rate[i]*tau)) であると想定されるため、人口数が負になる場合を含め、率が負の場合にエラーが発生することです。

# Parameter Values 
sir.parms <- c(deltaHinfinity=0.00299, CHi=0.00586, deltaH0=0.0854, aH=0.5,
               muH=0.02, SigmaW=0.1, SigmaM =0.8, SigmaL=104, phi=1.15, f = 0.6674,
               deltaVo=0.0166, CVo=0.0205, alphaVo=0.5968, beta=52, mbeta=7300 ,muV=52, g=0.0096, N=100)
# Inital Population Values
sir.x0 <- c(W=20,M=10,L=0.02)
# Rate Equations
sir.a <- c("((deltaH0+deltaHinfinity*CHi*mbeta*L)/(1+CHi*mbeta*L))*mbeta*L*N"
           ,"SigmaW*W*N", "muH*W*N", "((1/2)*phi*f)*W*N", "SigmaM*M*N", "muH*M*N",
           "(deltaVo/(1+CVo*M))*beta*M*N", "SigmaL*L*N", "muV*L*N", "alphaVo*M*L*N", "(aH/g)*L*N")
# Population change for even
sir.nu <- matrix(c(+0.01,0,0,
                   -0.01,0,0,
                   -0.01,0,0,
                   0,+0.01,0,
                   0,-0.01,0,
                   0,-0.01,0,
                   0,0,+0.01/230,
                   0,0,-0.01/230,
                   0,0,-0.01/230,
                   0,0,-0.01/230,
                   0,0,-0.01/32),nrow=3,ncol=11,byrow=FALSE)
runs <- 10
set.seed(1)

# Data Frame of output
sir.out <- data.frame(time=numeric(),W=numeric(),M=numeric(),L=numeric())
# Multiple runs and combining data and SSA methods 
for(i in 1:runs){
  sim <- ssa(sir.x0,sir.a,sir.nu,sir.parms, method="ETL", tau=1/12, tf=140, simName="SIR")
  sim.out <- data.frame(time=sim$data[,1],W=sim$data[,2],M=sim$data[,3],L=sim$data[,4])

  sim.out$run <- i
  sir.out <- rbind(sir.out,sim.out)
}

したがって、レートが計算され、モデルは各時間ステップの人口値を更新し、データ フレーム内のデータ ストアを使用して、以前の実行と一緒にアタッチします。ただし、集団のレベルが非常に低くなると、集団を減少させるイベントの数がコンパートメント内の数よりも多くなるようなイベントが発生する可能性があります。1 つの方法は、タイム ステップを非常に小さくすることですが、これによりシミュレーションの長さが非常に長くなります。

私の質問は、各時間ステップでデータが作成/計算されるときに、負の人口数の値が 0 に変換されるようにコードを拡張する方法はありますか?

私はこの問題に取り組んでみましたが、シミュレーションが完了すると値を変更する方法しか思いつかないようで、負の値は依然として実行自体に問題を引き起こしています。例
(sir.out$L < 0) sir.out$L == 0 の場合

どんな助けでもいただければ幸いです

4

1 に答える 1

2

問題は、ssa 関数で設定したメソッド (「ETL」) にあると思います。ETL メソッドは、最終的に負の数を生成します。タウ跳躍シミュレーション方法の効率的なステップ サイズの選択に基づいて、「OTL」方法を試すことができます。微調整できるパラメーターがさらにいくつかありますが、基本的なコマンドは次のとおりです。

ssa(sir.x0,sir.a,sir.nu,sir.parms, method="OTL", tf=140, simName="SIR")

または、負の数をまったく生成しない直接的な方法：

ssa(sir.x0,sir.a,sir.nu,sir.parms, method="D", tf=140, simName="SIR")
于 2016-09-13T15:53:31.697 に答える