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?

$\endgroup$
  • $\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 Apr 14 '18 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. Apr 14 '18 at 15:50
  • 1
    $\begingroup$ You could put a If[Abs[rate]<100. $MachineEpsilon,Break[]] before the point of division. $\endgroup$ – Henrik Schumacher Apr 14 '18 at 16:01
  • $\begingroup$ Have you seen this? mathematica.stackexchange.com/questions/158212/… $\endgroup$ – OkkesDulgerci Apr 14 '18 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 Apr 14 '18 at 17:49
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

$\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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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