C++を使用した出産・死亡プロセスによるCTMCシミュレーション

Question

数学は私の強いスキルではありません。しかし、私は誕生のための連続時間標識鎖（CTMC）の移行時間をシミュレートする必要があります&死のプロセスは、C + +を使用しています。

これは、すべてのラムダの行合計が1になる通常のCTMCをシミュレートするgithub projectです。しかし、出生死のプロセス（M/M/c/K）の場合、ゼロになります。だから私はそれを私の目的のために正確に使うことはできません。

M/M/c/K CTMCをシミュレートするアルゴリズムはどこにありますか？見つけたらアルゴリズムをコーディングすることができます。しかし、私は自分でアルゴリズムを構築することはできません、数学は私の頭の上に行く。

ポアソン分布を使用してM/M/c/Kキューにイベントを送信するには、このシミュレーションが必要です。そうすれば、最大到着率（ラムダ）の下でサーバ（c）の要件を把握することができます。

Answer 1

あなたの心が望むものはまだ100％明確ではありませんが、話題はかなり面白いです。

最初に、あなたが言及しているgithubプロジェクトは、正確に何が入力として期待されていると言えるほど十分に文書化されておらず、マルコフプロセスについての知識が不十分であると言うと、あまり意味がありません。

2番目：遷移率行列の行の合計は、出生死のプロセスだけでなく、すべてのマルコフプロセスで0です。

私は、実行をシミュレートする方法を次のとおりです。、（P'=PQ経由で定義された）開始状態（おそらく0）S、遷移率行列Q遷移関心のn数：

考える

1. を。

2. 出力：times - 遷移が発生した時刻 - 訪問された一連のステート。states - 一連の訪問ステート。

は、ここでは、C++と擬似コードの野生のミックスで行く：

std::default_random_engine rnd;//ini 

double current_time=0.0; 
int current_state=S; 
vector<doubles> times={current_time}; 
vector<int> states={current_state}; 

for (int i=0;i<n;i++){ 
//Part 1: simulate the time staying in current state: 
    double decay_rate=-Q[current_state][current_state]; 
    if(decay_rate==0.0){ 
    //that means we are not going anywhere anymore and staying for ever in this state: 
     return; 
    } 
    //we don't do error checking and expect rate to be positive, because diagonal elements of Q must be 0.0 or negative 
    //The current state will decay with the decay_rate, so the life time in this state until the decay is exponentially distributed with parameter decay_rate 
//simulate the life time: 
std::exponential_distribution<> life_time_generator(decay_rate); 
double life_time=life_time_generator(rnd); 
times.push_back(times.back()+life_time); 

//Part2: to which state have we actually decayed? 
// The probability to decay to the state new_state is proportional to transition rate Q[current_state][new_state] 
//thus generate an uniformly distributed random variable [0, decay_rate] (decay_rate - sum of transition rates of all possible new states) and map it on the states: 
double target_value=std::generate_canonical<double,10>(rnd)*decay_rate; 
double sum=0.0; 
for (int new_state=0;new_state<Q.size();new_state++){ 
    if (new_state==current_state)//don't forget to skip the state itself 
     continue; 
    sum+=Q[current_state][new_state]; 
    if (sum>target_value){//we found our next state! 
     current_state=new_state; 
     states.push_back(current_state); 
     break; 
    } 
    //there are still a some precision issues, if the sum of all transition rates is slightly under 1.0 
    //the issues should be handled somehow but not in this pseudo-code. 
} 

} 
// do whatever you want with times/states 

私はそれはあなたが考えていたものです願っています。

編集：

として簡単な説明：最初の部分のため

1. - 遷移するまでの時間。このコードは、既知の事実に基づいています。到着時間がレートλでポアソン分布の場合、待機時間はパラメータλで指数分布します。第二の部分のための例here

2. を参照 - それだけでconditional probability：時間dtの非常に短い期間のために遷移確率が-Q[current_state][current_state]dtあります。この状態は、遷移が起こったことを知っているときに満たされます。new_stateの状態になる確率はQ[current_state][new_state]*dtですが、遷移が発生したという条件の下ではQ[current_state][new_state]/-Q[current_state][current_state] - となり、2番目の部分で計算されます。