Message Boards Message Boards

[✓] Implement simple birth death process?

GROUPS:

Hi All, I am trying to implement simple birth death process. Why my code does not work? Any suggestion. Thanks.

birthDeath[λ_, μ_, initialPopulation_, numOfReaction_] :=     
 NestList[(
    Δt1 =  RandomVariate[ExponentialDistribution[λ  #]];
    Δt2 = RandomVariate[ExponentialDistribution[μ #]]; 
    If[Δt1 < Δt2, {# + 1}, {# - 1}]) &,  initialPopulation, numOfReaction]

birthDeath[3, 1, 10, 2]
POSTED BY: Okkes Dulgerci
Answer
1 month ago

I fixed it.

exact[\[Lambda]_, \[Mu]_, initialPopulation_] := 
 Module[{}, 
  DSolveValue[{p'[t] == (\[Lambda] - \[Mu]) p[t], 
    p[0] == initialPopulation}, p[t], t]]; 
birthDeath[\[Lambda]_, \[Mu]_, initialPopulation_, numOfReaction_] := 
 NestList[(\[CapitalDelta]t1 = 
     RandomVariate[ExponentialDistribution[\[Lambda] #[[2]]]];
    \[CapitalDelta]t2 = 
     RandomVariate[ExponentialDistribution[\[Mu] #[[2]]]];
    \[CapitalDelta]t = Min[\[CapitalDelta]t1, \[CapitalDelta]t2];
    {#[[1]] + \[CapitalDelta]t, 
     If[\[CapitalDelta]t1 < \[CapitalDelta]t2, #[[2]] + 1, #[[2]] - 
       1]}) &, {0, initialPopulation}, numOfReaction]


With[{\[Lambda] = 3, \[Mu] = 1, initialPopulation = 10, 
  numOfReaction = 1000, numOfsim = 10},
 sim = Table[
   birthDeath[\[Lambda], \[Mu], initialPopulation, numOfReaction], 
   numOfsim];
 Show[ListStepPlot[sim, PlotLegends -> {"Simulation"}, 
   PlotStyle -> Directive[AbsoluteThickness[0.2]], Frame -> True, 
   PlotTheme -> "Detailed", FrameLabel -> {"Time", "Population"}, 
   ImageSize -> Large], 
  Plot[Evaluate@exact[\[Lambda], \[Mu], initialPopulation], {t, 0, 
    Max@sim[[All, -1, 1]]}, PlotStyle -> {Black, Thick}, 
   PlotLegends -> {"ODE"}], 
  ListLinePlot[Mean@sim, PlotLegends -> {"Avg of Simulation"}, 
   PlotStyle -> {Red, Dashed}]]]
POSTED BY: Okkes Dulgerci
Answer
1 month ago

How can I store Δt=Min[Δt1,Δt2] so that I can plot time vs population

POSTED BY: Okkes Dulgerci
Answer
1 month ago

This may be one way:

birthDeath2[\[Lambda]_, \[Mu]_, initialPopulation_, numOfReaction_] :=
  NestList[(\[CapitalDelta]t1 = 
     RandomVariate[ExponentialDistribution[\[Lambda] #[[3]]]];
    \[CapitalDelta]t2 = 
     RandomVariate[ExponentialDistribution[\[Mu] #[[3]]]];
    {\[CapitalDelta]t1, \[CapitalDelta]t2, 
     If[\[CapitalDelta]t1 < \[CapitalDelta]t2, #[[3]] + 1, #[[3]] - 
       1]}) &, {0, 0, initialPopulation}, numOfReaction]

birthDeath2[3, 1, 10, 4]
POSTED BY: Gianluca Gorni
Answer
1 month ago

Thanks

POSTED BY: Okkes Dulgerci
Answer
1 month ago

Group Abstract Group Abstract