I am trying to implement simple birth death process. Why my code does not work? Any suggestion. Thanks. Cross posted http://community.wolfram.com/groups/-/m/t/1205656
birthDeath[λ_, μ_, initialPopulation_, numOfReaction_] :=
NestList[(
Δt1 = RandomVariate[ExponentialDistribution[λ #]];
Δt2 = RandomVariate[ExponentialDistribution[μ #]];
If[Δt1 < Δt2, {# + 1}, {# - 1}]) &, initialPopulation, numOfReaction]
birthDeath[3, 1, 10, 2]
Answer
How can I store Δt = Min[Δt1, Δt2] so that I can plot time vs population
One way is
birthDeath2[λ_, μ_, initialPopulation_, numOfReaction_] :=
NestList[(Δt1 = RandomVariate[ExponentialDistribution[λ #[[2]]]];
Δt2 = RandomVariate[ExponentialDistribution[μ #[[2]]]];
{Min[Δt1, Δt2], If[Δt1 < Δt2, #[[2]] + 1, #[[2]] - 1]}) &,
{0, initialPopulation}, numOfReaction]
{Δt, pop} = Transpose[birthDeath2[3, 1, 10, 10]];
ListLinePlot[Transpose[{Accumulate@Δt, pop}], InterpolationOrder -> 0]
Comments
Post a Comment