CTMC Simulation mit Geburt und Tod Prozess mit C++

Question

Math ist nicht meine starke Fähigkeit. Allerdings muss ich die Übergangszeit Markov Kette (CTMC) Übergangszeiten für die Geburt & Tod Prozess mit C++ simulieren.

Ich stieß auf diese github project, die regelmäßige CTMC simuliert, wo die Zeilensumme aller Lambda wird 1. Aber im Fall der Geburt-Tod-Prozess (M/M/c/K) wird es Null sein. Also kann ich es nicht genau für meinen Zweck verwenden.

Wo finde ich den Algorithmus zur Simulation von M/M/c/K CTMC? Ich kann den Algorithmus programmieren, wenn ich ihn finde. Aber ich bin nicht in der Lage, den Algorithmus selbst aufzubauen, die Mathematik geht über meinen Kopf.

Ich brauche diese Simulation, um Ereignisse an eine M/M/c/K-Warteschlange zu senden, mit Poisson-Verteilung. Auf diese Weise kann ich die Server (c) -Anforderung unter verschiedenen Ankunftsraten (Lambda) ermitteln und gleichzeitig eine maximale Serverauslastung sicherstellen.

Answer 1

Es ist immer noch nicht 100% klar, was dein Herz begehrt, aber das Thema ist ziemlich interessant.

Erstens: Das Github-Projekt, auf das Sie sich beziehen, ist nicht gut genug dokumentiert, um zu sagen, was genau es als Input erwartet und ich habe nicht genug Wissen über Markov-Prozesse, um zu sagen, dass der kontinuierliche Teil falsch ist, aber für mich macht nicht viel Sinn.

Zweitens: Die Summe einer Zeile einer Übergangsrate Matrix ist 0 für alle Markov-Prozesse, nicht nur für Geburt-Tod-Prozesse.

Das ist, wie ich einen Lauf simulieren würde:

1. Gegeben: der Startzustand (wahrscheinlich 0) S, Übergangsrate Matrix Q (definiert über P'=PQ), die Anzahl der Übergänge n von Interesse.

2. Ausgabe: times - Zeiten, zu denen die Übergänge aufgetreten sind, states - die Reihe der besuchten Staaten.

Hier gehen wir mit einer wilden Mischung aus C++ und Pseudo-Code:

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 

Ich hoffe, dass es etwas, das Sie im Sinn hatte.

Edit:

Als kurze Erklärung:

1. Für den ersten Teil - die Zeit bis zum Übergang. Der Code basiert auf der bekannten Tatsache, dass, wenn die Ankünfte Poisson mit der Rate Lambda verteilt sind, die Wartezeiten exponentiell mit dem Parameter Lambda verteilt sind. Siehe zum Beispiel here

2. Für den zweiten Teil - es ist nur die conditional probability: Die Übergangswahrscheinlichkeit für einen sehr kurzen Zeitraum dt ist -Q[current_state][current_state]dt. Diese Bedingung ist erfüllt, da wir wissen, dass der Übergang stattgefunden hat.Die Wahrscheinlichkeit, in den Zustand new_state zu gelangen, ist Q[current_state][new_state]*dt, aber unter der Bedingung, dass der Übergang passiert ist, ist es Q[current_state][new_state]/-Q[current_state][current_state] - und das ist, was im zweiten Teil berechnet wird.