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 count 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 = 100000
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)
My lecturer only solved this in Mathematica with the code below but I want to do it in Matlab or R, however I can't seem to translate the code properly.
In[1]:={repetitions, count}={1000000, 0};
For[i=1, i<=repetitions, i++, X=1; time=0;
While[X>0 && time<=4,
If[X==1, time=time+Random[ExponentialDistribution[2]];
If[Random[UniformDistribution[{0,1}]]<1/2, X=0, X=2],If[X==2, time=time+Random[ExponentialDistribution[3]];
If[Random[UniformDistribution[{0,1}]]<2/3, X=1, X=3],If[X==3, time=time+Random[ExponentialDistribution[4]];
If[Random[UniformDistribution[{0,1}]]<3/4, X=2, X=4],If[X==4, time=time+Random[ExponentialDistribution[3]]; X=3]]]]];
If[time>4, count=count+1]];
N[count/repetitions]
Out[1]= 0.117121
Figured it out. The last if-statement is wrong. If x = 4, I should update the time and set X = 3. So it should be
instead of what I've written in my question.