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