2
$\begingroup$

Question:

A birth and death process is a continuous time Markov chain.

Find an approximative numerical value for the probability

P {max0≤t≤10 X(t) ≥ 10}

for a birth and death process {X(t)}t≥0

with birth intensities λ0 = λ1 = λ2 = . . . = 1

and death intensities µ1 = µ2 = µ3 = . . . = 2

that is in steady-state (that is, which is started according to its stationary distribution).

My attempt at a solution:

rep = 100000

count = 0


for (i in 1:rep) {

   time = 0
   xt = 0
   succ = 0
   start = 0

   u.rand = runif(1)

   # These if statements determine where we start 
   # according to the stationary distribution

   if (u.rand<1/2) {start = 0}                   

   if (1/2<=u.rand && u.rand<3/4) {start = 1}

   if (3/4<=u.rand && u.rand<7/8) {start = 2}

   if (7/8<=u.rand && u.rand<15/16) {start = 3}

   if (15/16<=u.rand && u.rand<31/32) {start = 4}

   if (31/32<=u.rand && u.rand<63/64) {start = 5}

   if (63/64<=u.rand && u.rand<127/128) {start = 6}

   if (127/128<=u.rand && u.rand<255/256) {start = 7}

   if (255/256<=u.rand && u.rand<511/512) {start = 8}

   if (511/512<=u.rand && u.rand<1023/1024) {start = 9}

   if (1023/1024<=u.rand) {start = 10}

   # xt is the position during the while loop
   # and its a success if xt gets bigger than 10 

   xt = start                                          

   while ((time<=10) && (succ==0)) {

    mu = rexp(1, 2)                                   # death rate
    la = rexp(1, 1)                                   # birth rate


    # Here we get a birth and xt increases by 1.
    if (la<mu) {                                      

    xt = xt + 1

    time = time + la

    # Here we get a death as long as xt not is equal to 0

    } else if ((mu<la) && (xt!=0)) {                 

    xt = xt - 1

    time = time + mu

    # at xt=0 we can only get births

     } else if (xt==0) {                               

     xt = xt + 1

     time = time + la }


      # Here we register the successes.

      if (xt>=10) {succ = 1}                              

}

count = count + succ                  # And here we add them up

}

count/rep                                             

# And this frequency should be equal to
# the sought after probability

I get count/rep = 0.00615 but it should be 0.0084865

Something is wrong in my code.

Can anyone help me out?

Im pretty sure im doing something wrong in the while-loop.

$\endgroup$
9
  • $\begingroup$ Please give a full reference. I haven't the faintest idea what "G-S's book" might be. If you'd written G & S I might have guessed Grimmett and Stirzaker, but even then it would only be a guess and there are no doubt other possibilities. Please don't make people guess. We're not mindreaders. $\endgroup$
    – Glen_b
    Dec 13, 2015 at 5:36
  • $\begingroup$ Yes, its Grimmett and Stirzaker, I should have left that part out. I stated the problem as it is written and including Grimmett and Stirzaker was unnecessary. You dont need to know it is Grimmett and Stirzaker. You only need to know its a birth and death process with the specified birth and death rates and which probability im looking to simulate. $\endgroup$
    – JKnecht
    Dec 13, 2015 at 8:26
  • $\begingroup$ I removed that part now. $\endgroup$
    – JKnecht
    Dec 13, 2015 at 8:52
  • $\begingroup$ If you're quoting part of a book, you do need to include the full reference. See the policy on referencing. You would also need to include the self-study tag (and read its tag-wiki, if necessary modifying your question as described). However, it appears as if you're simply asking for help with coding, in which case your post seems to be off-topic, as discussed in the on topic help (see under Programming). ... ctd $\endgroup$
    – Glen_b
    Dec 13, 2015 at 11:22
  • $\begingroup$ ctd ... If you believe it's on topic in spite of that, you'd need to edit to clarify in accordance with that policy. $\endgroup$
    – Glen_b
    Dec 13, 2015 at 11:32

1 Answer 1

0
$\begingroup$

First of all, are you sure that 0.0085 is the correct answer? Looking at here suggests otherwise???

Some comments about your code:

(1) I think your code is essentially doing the right thing. However, the structure of the code within your while loop makes it difficult to read. I would split up the code into two 'if' clauses, i.e. if(xt==0) {...} and if(xt!=0) {...}, reflecting the special nature of the xt==0 situation, so that it is clear what actions are performed for these two different cases.

(2) Be careful as it is possible that 'time' may exceed 10 in your while loop, so you would need to ask that time<=10 as well as xt>=10 in order to register a success.

(3) u.rand>=1023/1024 corresponds to the situation where start>=10. Surely in these cases success is already guaranteed i.e. you should set succ=1.

(4) The statements about u.rand could be combined into a for loop, which would save typing time and reduce the risk of typing error.

(5) Running your code more than once reveals substantial between-run variability. You could reduce this by increasing 'rep'. Even better, you could stratify your sampling of the stationary distribution, as at the moment your use of 'runif' is introducing additional variability into your answers. To stratify, you could start by specifying rep=1024 (for example), then for the first 512 reps set start=0, for the next 256 reps set start=1, and so on.

$\endgroup$
3
  • $\begingroup$ +1 Thanks! Im learning a lot!. Doing the first 3 of your suggestions and increasing rep give me the result P = 0.0064 and that is the solution you are linking to. But i was under the impression the simulation on this page was correct and there the answer is 0.0085. Simulation link $\endgroup$
    – JKnecht
    Dec 13, 2015 at 12:28
  • $\begingroup$ And he is doing something different than me there, i dont understand what tho...and therefore i was thinking my solution was totally wrong...so im a bit confused. (Now im also trying to change my code according to your comments 4 and 5 which surely will improve it(but not change it so it looks like the linked solution though)) $\endgroup$
    – JKnecht
    Dec 13, 2015 at 12:31
  • $\begingroup$ @JKnecht I don't fully understand the code you linked to in your comment - written in Wolfram I think. It's certainly much more difficult to read than your code, but I think that it is basically doing the same thing that your code is doing, and all of the points (1)-(5) in my answer apply to that code as well. I wouldn't worry about trying to get 0.0085 as it's not clear that this is the correct answer. Instead I would focus on understanding your code and making it as efficent as possible. $\endgroup$ Dec 13, 2015 at 15:05

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