I have the following problem:
Consider an M/M/3/4 queuing system with $\lambda=\mu=1$ that is the arrival time is exponentially distributed with parameter $\lambda = 1$ and the service times are exponentially distributed with parameter $\mu=1$. A busy time $B$ for the system is the random time it takes from a customer comes into the system when it is empty until the system gets empty again. Write a computer program that by means of stochastic simulation finds an approximate value of the probability $\Pr(B>4)$.
So the basic idea of the code below is to count the occurrences of $B>4$ and divide by the number of simulations. X
denotes the possible states. My problem is that the last if-statement does not seem to work as it should. If runif(1)<1/2
then X=0
, so the while loop should exit. Now it's very unlikely that time
is greater than 4
so countcount
should not be updated. But it's always updated no matter what. Note that time
is the same thing as $B$. I can't find my error. The answer should be
$$\Pr(B>4)\approx 0.117121$$
nrsim = 100
count = 0
for (k in 1:nrsim) {
time = 0
X = 1
while (X > 0 && time <= 4) {
if (X == 1){
time = time + rexp(1,2)
if (runif(1) < 1/2) {
X = 0
}
else {X = 2}
}
if (X == 2){
time = time + rexp(1,3)
if (runif(1) < 2/3) {
X = 1
}
else {X = 3}
}
if (X == 3){
time = time + rexp(1,4)
if (runif(1) < 3/4) {
X = 2
}
else {X = 4}
}
if (X == 4){
time = time + rexp(1,3)
}
else {X = 3}
}
if (time > 4){
count = count + 1
}
}
print(count/nrsim)