0
$\begingroup$

I am simulating a birth-death population process. My code is 

RESULTS = {};
For[iterator = 1, iterator < 21, iterator++,
 n = 25; t = 0;
 results = {{t, n}};
 mu = 1; lam = 1;
 While[t < 200,
  death = lam*n;
  birth = mu*n;
  rate = death + birth;
  deltaT[r_] := -1/r*Log[RandomReal[]];
  t1 = deltaT[rate];
  t = t + t1;
  rand = RandomReal[]*rate;
  Which[
   rand <= death, n = n - 1,
   True, n = n + 1
   ];
  results = Append[results, {t, n}];
  ];
 AppendTo[RESULTS, results];
 ]

I am getting an error message that says 'Infinite expression 1/0 encountered. >>'.

I'm not quite sure why this is. I think it might be because the population dies out so the rate would be zero, hence deltaT[r_] would be trying to work out -1/0.

Is this correct?

How do I get around this?

Improve this question
$\endgroup$
5
  • $\begingroup$ I understand that's the reason why. I therefore want the simulation to run up until the point to which the rate becomes 0, then stop. How do I adapt my code so that it does this? @HenrikSchumacher $\endgroup$
    – Mlo27
    Commented Apr 14, 2018 at 15:49
  • $\begingroup$ @Mlo27 For rate to be zero, you have to have n = 0, right? So you need to change your Which statement so that n does not become zero. $\endgroup$
    – C. E. ♦
    Commented Apr 14, 2018 at 15:50
  • 1
    $\begingroup$ You could put a If[Abs[rate]<100. $MachineEpsilon,Break[]] before the point of division. $\endgroup$
    – Henrik Schumacher
    Commented Apr 14, 2018 at 16:01
  • $\begingroup$ Have you seen this? mathematica.stackexchange.com/questions/158212/… $\endgroup$
    – OkkesDulgerci
    Commented Apr 14, 2018 at 17:32
  • $\begingroup$ @OkkesDulgerci I have. I don’t want to change code massively though. I just want to add a command so that the simulation stops when I get the error. $\endgroup$
    – Mlo27
    Commented Apr 14, 2018 at 17:49

2 Answers 2

Reset to default
2
$\begingroup$

As @Henrik mentioned in comment, you can break simulation if rate<0.

RESULTS = {};
For[iterator = 1, iterator < 21, iterator++, n = 25; t = 0;
 results = {{t, n}};
 mu = 1; lam = 1;
 While[t < 200, death = lam*n;
  birth = mu*n;
  rate = death + birth;
  If[Abs[rate] < 100. $MachineEpsilon, Break[]];
  deltaT[r_] := -1/r*Log[RandomReal[]];
  t1 = deltaT[rate];
  t = t + t1;
  rand = RandomReal[]*rate;
  Which[rand <= death, n = n - 1, True, n = n + 1];
  results = Append[results, {t, n}];];
 AppendTo[RESULTS, results];]

ListLinePlot[RESULTS, InterpolationOrder -> 0]

enter image description here

Improve this answer
$\endgroup$
0
$\begingroup$

I have tried to make a functional version of the original code that incorporates the observation made in the comments by multiple discussants that n=0 breaks the original code.

This implementation uses NestWhileList to substitute for the innermost While and NestList replaces the outermost For. 

Block[{test, compose},

  With[{no = 25, to = 0, te = 200, mu = 1, lam = 1, iterator = 20, seedo = 123456987},

    (* input {t, n, r}; returns 'True' is t < te (=200) and r \[NotEqual] 0 *)
    test = #1 < te && #3 != 0 &;

    (* input {t, n}; returns {t, n, r} *)
    compose = {##, (lam + mu) #2} &;

    (* get rid of the first entry in the output (= NestList's second argument (check documentation) *)
    data = Rest@NestList[

      (* create the seed to be used in BlockRandom (used for reproducibility) *)
      With[{seed = #[[1]] + 1},

        (* output = {{seed, {{t, n, r},.. }}.. } *)
        {seed, BlockRandom[

          NestWhileList[

            With[{t = #[[1]], n = #[[2]], r = #[[-1]], rands = RandomReal[{0, 1}, 2]},

              (* perhaps redundant With - helps readability *)
              With[{death = lam n},

               (* produce next {t, n, r} from {t, n} *)
               Apply[compose][{
                 t - 1/r Log[rands[[1]]],
                 Which[
                   r rands[[-1]] <= death, n - 1,
                   True, n + 1
                  ]
                }]

               ]
             ] &,

            (* init = {t, n, r} *)
            Apply[compose][{to, no}],

            (* check if t < 200 and r\[NotEqual]0 *)
            Apply[test][#] &], RandomSeeding -> #[[1]]]}

     (* init NestList *)
     ] &, {seedo, {}}, iterator];

  ]

 ]

enter image description here

The plots were produced with:

ListLinePlot[
  Part[#, -1, All, {1, 2}],
  Frame -> True,
  PlotLegends -> Placed[Row[{HoldForm[seed] -> Part[#, 1]}], Bottom],
  InterpolationOrder -> 0,
  PlotRange -> {{0, 200}, Automatic}
 ] & /@ data // Partition[#, 4, 4, {1, 1}, {}] & // Grid
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.