# R program to simulate an M/M/1 queue with mu=3 and lambda=2 newservicetime = function() { return( rexp(1,3) ) } newinterarrivaltime = function() { return( rexp(1,2) ) } TARGETTIME = 10000; MAXGRAPH = 20; timetoarrive = newinterarrivaltime(); timetoserve = 0; timeremaining = TARGETTIME; thetimes = rep(0,MAXGRAPH+1); Q = 0; while (timeremaining > 0) { if ( (timeremaining <= timetoarrive) && (timeremaining <= timetoserve) ) { # Finish without any further events. if (Q <= MAXGRAPH) thetimes[Q+1] = thetimes[Q+1] + timeremaining; timeremaining = 0; } else if ( (Q==0) || (timetoarrive < timetoserve) ) { # Await next arrival. timeremaining = timeremaining - timetoarrive; timetoserve = timetoserve - timetoarrive; if (Q <= MAXGRAPH) thetimes[Q+1] = thetimes[Q+1] + timetoarrive; if (Q==0) # the new arrival sees an empty queue timetoserve = newservicetime(); Q = Q+1; timetoarrive = newinterarrivaltime(); } else { # Await next service completion. timeremaining = timeremaining - timetoserve; timetoarrive = timetoarrive - timetoserve; if (Q <= MAXGRAPH) thetimes[Q+1] = thetimes[Q+1] + timetoserve; Q = Q-1; timetoserve = newservicetime(); } # cat(Q, timeremaining, timetoserve, timetoarrive, "\n"); } thefracs = thetimes / TARGETTIME; cat("SIMULATED PROBABILITIES:\n", thefracs, "\n"); plot((0:MAXGRAPH), thefracs, type='p', col="blue", xlab="Q", ylab="probs"); # For comparison, in original M/M/1 case only: thetrue = rep(0,MAXGRAPH+1); r=2/3; for (i in (0:MAXGRAPH)) thetrue[i+1] = (1-r)*r^i; cat("TRUE PROBABILITIES:\n", thetrue, "\n"); points((0:MAXGRAPH), thetrue, type='l', col="red");