M/M/3/4 queue in R, translating code from mathematica.

98 Views Asked by At

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
1

There are 1 best solutions below

0
On BEST ANSWER

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

if (X == 4){
  time = time + rexp(1,3)
  X = 3
}

instead of what I've written in my question.