# Trying to simulate a birth death process [closed]

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.

• 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. Dec 13, 2015 at 5:36
• 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. Dec 13, 2015 at 8:26
• I removed that part now. Dec 13, 2015 at 8:52
• 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 Dec 13, 2015 at 11:22
• ctd ... If you believe it's on topic in spite of that, you'd need to edit to clarify in accordance with that policy. Dec 13, 2015 at 11:32

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

(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.

• +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 Dec 13, 2015 at 12:28
• 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)) Dec 13, 2015 at 12:31
• @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. Dec 13, 2015 at 15:05