Message Boards Message Boards

0
|
4867 Views
|
4 Replies
|
1 Total Likes
View groups...
Share
Share this post:

[?] Implement simple birth death process?

Posted 6 years ago

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
4 Replies

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
Posted 6 years ago

Thanks

POSTED BY: Okkes Dulgerci
Posted 6 years 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
Posted 6 years ago

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

POSTED BY: Okkes Dulgerci
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract