# [✓] Implement simple birth death process?

Posted 1 year ago
806 Views
|
4 Replies
|
1 Total Likes
|
 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] 
4 Replies
Sort By:
Posted 1 year 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 1 year ago
 How can I store Δt=Min[Δt1,Δt2] so that I can plot time vs population
 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]