#.rs.restartR()
#=====================================================================
# Modelled and coded by Ebert Brea
#=====================================================================
cat("\f")
rm(list=ls())
library(xtable)
require(latex2exp)
## Loading required package: latex2exp
## Warning: package 'latex2exp' was built under R version 4.1.2
require(graphics)
library(stargazer)
## Warning: package 'stargazer' was built under R version 4.1.2
##
## Please cite as:
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
simulaMM1K <- function(lambda,mu,cap,ElapsedTime,NoSeed){
# ========================================================
# INPUT
# lambda: arrival rate
# mu: service rate
# ElapsedTime: period time simulation
# NoSeed: number of seed
# =======================================================
# OUTPUT
# DropPerformance:
# 1) Total number of dropped customers
# 2) Total number of finished customers
# 3) Mean of system time
# 4) Mean of queue length
#=======================================================
# For testing function simulaMM1K
#lambda=10; mu=lambda/0.8; cap=5; ElapsedTime=400; NoSeed=20
set.seed(NoSeed)
library(simmer)
library(parallel)
env<-simmer()
BetweenArrival<- function(){rexp(1,lambda)}
ServiceTime<- function(){rexp(1,mu)}
CapacityQueue<-cap
path<-function(){ ifelse( get_queue_count(env, "Server") < CapacityQueue ,1, 2) }
#===============================================================================#%
Start<-trajectory("Start")%>%
set_global("Counter", 1, mod="+",init=-1) %>%
set_attribute("My_ID", function() get_global(env, "Counter") )%>%
set_attribute("QueueLength", function() get_queue_count(env, "Server") )%>%
set_attribute("Time_in",function(){now(env)})
Q1<-trajectory("Queue 1") %>%
seize("Server", amount=1) %>%
timeout(ServiceTime) %>%
release("Server", amount=1)
Q2<-trajectory()%>%
set_global("Leave", 1, mod="+",init=0)
Network<-trajectory()%>%
join(Start)%>%
branch(option = path,
continue=c(TRUE,TRUE),
join(Q1),
join(Q2))
env<-simmer()
env<-env %>%
add_resource("Server", capacity=1, queue_size = CapacityQueue) %>%
add_generator("Customer", Network, BetweenArrival, mon=2) %>%
run(ElapsedTime)
release_all(Network, "Server")
resources <- get_mon_resources(env)
arrivals <- get_mon_arrivals(env,per_resource=TRUE)
attribites <-get_mon_attributes(env)
Arrival.df<-data.frame(get_mon_arrivals(env))
Attribute.df<-data.frame(get_mon_attributes(env))
Ndata<-length(Attribute.df[,4])
tmp<-numeric(Ndata)
QL.Obs<-1
for(ii in 1:Ndata){
if(Attribute.df$key[ii]=="QueueLength"){tmp[QL.Obs]<-Attribute.df$value[ii];QL.Obs<-QL.Obs+1}
}
QL<-numeric(QL.Obs)
for(i in 1:QL.Obs){QL[i]<-tmp[i]}
DropPerformance <- c( sum(Attribute.df$key=="Leave"),
sum(Arrival.df$finished),
mean(Arrival.df$end_time - Arrival.df$start_time),
mean(QL))
return(DropPerformance)
}
#==========================================================
# Multireplications
#==========================================================
NoReplications <- 30
L<-10
Rho<-0.8
#=======================================================================================
# For testing the simulation model
#=======================================================================================
Result <- matrix(data = NA, nrow = NoReplications, ncol = 4, byrow = TRUE)
for(r in 1:NoReplications){
Result[r,] <-simulaMM1K(lambda=L, mu=L/Rho, cap=4, ElapsedTime=400,NoSeed=r)
}
## Warning: package 'simmer' was built under R version 4.1.2
df.Result <- data.frame(Result)
names(df.Result) <- c("Leave","Finish","System.Time","Queue.Length")
summary(df.Result$Leave)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 289.0 342.2 361.0 366.0 392.8 445.0
summary(df.Result$System.Time)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1772 0.1838 0.1856 0.1857 0.1892 0.1933
summary(df.Result$Queue.Length)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.021 1.104 1.159 1.142 1.181 1.241
var(df.Result$Queue.Length)
## [1] 0.002930475
hist(df.Result$Leave,main="Dropped customers",
xlab="Number of customers",
col=blues9,labels=TRUE)
RatioDropped <-df.Result$Leave/df.Result$Finish
df.Result <- data.frame(RatioDropped,Result)
names(df.Result) <- c("Leave/Finish", "Leave","Finish","System.Time","Queue.Length")
summary(RatioDropped)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.07294 0.08540 0.09084 0.09112 0.09761 0.10792
t.test(df.Result$System.Time)
##
## One Sample t-test
##
## data: df.Result$System.Time
## t = 251.22, df = 29, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.1842235 0.1872477
## sample estimates:
## mean of x
## 0.1857356
t.test(df.Result$Queue.Length)
##
## One Sample t-test
##
## data: df.Result$Queue.Length
## t = 115.58, df = 29, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.122161 1.162589
## sample estimates:
## mean of x
## 1.142375
t.test(L*df.Result$System.Time, df.Result$Queue.Length,
paired = TRUE, var.equal = TRUE, conf.level = 0.95)
##
## Paired t-test
##
## data: L * df.Result$System.Time and df.Result$Queue.Length
## t = 107.03, df = 29, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.7013178 0.7286435
## sample estimates:
## mean of the differences
## 0.7149807
#========================================================
stargazer(df.Result, type = "text",
summary.stat = c("min", "p25", "median", "p75", "max", "mean", "sd"),
title="Performance of a queue M/M/1/k")
##
## Performance of a queue M/M/1/k
## ====================================================================
## Statistic Min Pctl(25) Median Pctl(75) Max Mean St. Dev.
## --------------------------------------------------------------------
## Leave/Finish 0.073 0.085 0.091 0.098 0.108 0.091 0.009
## Leave 289 342.2 361 392.8 445 366.033 39.146
## Finish 3,896 3,963.8 3,995 4,062.5 4,146 4,013.400 64.792
## System.Time 0.177 0.184 0.186 0.189 0.193 0.186 0.004
## Queue.Length 1.021 1.104 1.159 1.181 1.241 1.142 0.054
## --------------------------------------------------------------------
stargazer(df.Result, type = "latex",
summary.stat = c("min", "p25", "median", "p75", "max", "mean", "sd"),
title="Performance of a queue M/M/1/k")
##
## % Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
## % Date and time: jue., mar. 24, 2022 - 6:25:26 p. m.
## \begin{table}[!htbp] \centering
## \caption{Performance of a queue M/M/1/k}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lccccccc}
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## Statistic & \multicolumn{1}{c}{Min} & \multicolumn{1}{c}{Pctl(25)} & \multicolumn{1}{c}{Median} & \multicolumn{1}{c}{Pctl(75)} & \multicolumn{1}{c}{Max} & \multicolumn{1}{c}{Mean} & \multicolumn{1}{c}{St. Dev.} \\
## \hline \\[-1.8ex]
## Leave/Finish & 0.073 & 0.085 & 0.091 & 0.098 & 0.108 & 0.091 & 0.009 \\
## Leave & 289 & 342.2 & 361 & 392.8 & 445 & 366.033 & 39.146 \\
## Finish & 3,896 & 3,963.8 & 3,995 & 4,062.5 & 4,146 & 4,013.400 & 64.792 \\
## System.Time & 0.177 & 0.184 & 0.186 & 0.189 & 0.193 & 0.186 & 0.004 \\
## Queue.Length & 1.021 & 1.104 & 1.159 & 1.181 & 1.241 & 1.142 & 0.054 \\
## \hline \\[-1.8ex]
## \end{tabular}
## \end{table}
hist(RatioDropped,
include.lowest = TRUE,
right = TRUE,
breaks = nclass.Sturges,
main=paste("Ratio of dropped customers with","Rho =",Rho ),
xlab="Ratio",
col=rgb(0.5, 0.9, 0),
labels=TRUE,
prob=FALSE)
# =========================================================
# Experiments of dropped Customers
# =========================================================
NoReplications <- 30
mu <- seq(2,5,by=1)
lambda<-seq(2,5,by=1)
k<-4
SimTime<-400
#Nr<-3
N1<-length(lambda)
N2<-length(mu)
#========================================================
# For writing a statistical report
#========================================================
LoDropped <-numeric(N1*N2)
LoAttended <-numeric(N1*N2)
LoRatioDropAtt <-numeric(N1*N2)
LoSystemTime <-numeric(N1*N2)
LoQueueLength <-numeric(N1*N2)
mDropped <-numeric(N1*N2)
mAttended <-numeric(N1*N2)
mRatioDropAtt <-numeric(N1*N2)
mSystemTime <-numeric(N1*N2)
mQueueLength <-numeric(N1*N2)
UpDropped <-numeric(N1*N2)
UpAttended <-numeric(N1*N2)
UpRatioDropAtt <-numeric(N1*N2)
UpSystemTime <-numeric(N1*N2)
UpQueueLength <-numeric(N1*N2)
Mu <- numeric(N1*N2)
Lambda <- numeric(N1*N2)
MultiResult = array(data=NA, dim=c(N1*N2,4,NoReplications))
# MultiResult is a 3D matrix, e.i., num [1:N1*N2, 1:3, 1:NoReplications]
n <- 0
for(i in 1:N1){
for(j in 1:N2){
n <- n+1
for(r in 1:NoReplications){
MultiResult[n,,r] <- simulaMM1K(lambda=lambda[i], mu=mu[j], cap=k, ElapsedTime=SimTime, NoSeed=r)
}
Lambda[n] <- lambda[i]
Mu[n] <- mu[j]
}
}
for(n in 1:(N1*N2)){
LoDropped[n] <-t.test(MultiResult[n,1,])$conf.int[1]
LoAttended[n] <-t.test(MultiResult[n,2,])$conf.int[1]
LoRatioDropAtt[n] <- t.test( (MultiResult[n,2,]-MultiResult[n,1,])/MultiResult[n,2,] )$conf.int[1]
LoSystemTime[n] <- t.test(MultiResult[n,3,])$conf.int[1]
LoQueueLength[n] <- t.test(MultiResult[n,4,])$conf.int[1]
mDropped[n] <-mean(MultiResult[n,1,])
mAttended[n] <-mean(MultiResult[n,2,])
mRatioDropAtt[n] <- mean( (MultiResult[n,2,]-MultiResult[n,1,])/MultiResult[n,2,] )
mSystemTime[n] <- mean(MultiResult[n,3,])
mQueueLength[n] <- mean(MultiResult[n,4,])
UpDropped[n] <-t.test(MultiResult[n,1,])$conf.int[2]
UpAttended[n] <-t.test(MultiResult[n,2,])$conf.int[2]
UpRatioDropAtt[n] <- t.test( (MultiResult[n,2,]-MultiResult[n,1,])/MultiResult[n,2,] )$conf.int[2]
UpSystemTime[n] <- t.test(MultiResult[n,3,])$conf.int[2]
UpQueueLength[n] <- t.test(MultiResult[n,4,])$conf.int[2]
}
df.ST <- data.frame(Lambda,Mu,mDropped,mAttended,mRatioDropAtt,mSystemTime,mQueueLength)
print(df.ST)
## Lambda Mu mDropped mAttended mRatioDropAtt mSystemTime mQueueLength
## 1 2 2 127.533333 793.5000 0.8397499 1.2443727 1.6281808
## 2 2 3 37.800000 795.2667 0.9528720 0.7023446 0.7698396
## 3 2 4 12.633333 802.0333 0.9843660 0.4519107 0.4073608
## 4 2 5 4.466667 801.3333 0.9944395 0.3206152 0.2421428
## 5 3 2 439.766667 1198.9000 0.6335238 1.1841988 2.6144460
## 6 3 3 197.133333 1198.5000 0.8358512 0.8229626 1.6362634
## 7 3 4 85.466667 1199.4667 0.9289775 0.5595536 0.9824714
## 8 3 5 41.200000 1205.2667 0.9659496 0.4021676 0.6351770
## 9 4 2 820.000000 1604.8667 0.4892271 1.0182058 3.1074902
## 10 4 3 484.233333 1599.3333 0.6975237 0.8154854 2.3460515
## 11 4 4 265.733333 1601.8667 0.8343836 0.6194764 1.6496477
## 12 4 5 141.733333 1599.9667 0.9116600 0.4631760 1.1277345
## 13 5 2 1199.866667 2000.0333 0.4002397 0.8656055 3.3474306
## 14 5 3 851.766667 2008.6667 0.5761290 0.7557719 2.8328948
## 15 5 4 547.433333 2008.5333 0.7275063 0.6234108 2.2152214
## 16 5 5 330.166667 2000.6333 0.8352105 0.4950510 1.6466410
print(xtable(df.ST), include.rownames = TRUE)
## % latex table generated in R 4.1.0 by xtable 1.8-4 package
## % Thu Mar 24 18:26:09 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrr}
## \hline
## & Lambda & Mu & mDropped & mAttended & mRatioDropAtt & mSystemTime & mQueueLength \\
## \hline
## 1 & 2.00 & 2.00 & 127.53 & 793.50 & 0.84 & 1.24 & 1.63 \\
## 2 & 2.00 & 3.00 & 37.80 & 795.27 & 0.95 & 0.70 & 0.77 \\
## 3 & 2.00 & 4.00 & 12.63 & 802.03 & 0.98 & 0.45 & 0.41 \\
## 4 & 2.00 & 5.00 & 4.47 & 801.33 & 0.99 & 0.32 & 0.24 \\
## 5 & 3.00 & 2.00 & 439.77 & 1198.90 & 0.63 & 1.18 & 2.61 \\
## 6 & 3.00 & 3.00 & 197.13 & 1198.50 & 0.84 & 0.82 & 1.64 \\
## 7 & 3.00 & 4.00 & 85.47 & 1199.47 & 0.93 & 0.56 & 0.98 \\
## 8 & 3.00 & 5.00 & 41.20 & 1205.27 & 0.97 & 0.40 & 0.64 \\
## 9 & 4.00 & 2.00 & 820.00 & 1604.87 & 0.49 & 1.02 & 3.11 \\
## 10 & 4.00 & 3.00 & 484.23 & 1599.33 & 0.70 & 0.82 & 2.35 \\
## 11 & 4.00 & 4.00 & 265.73 & 1601.87 & 0.83 & 0.62 & 1.65 \\
## 12 & 4.00 & 5.00 & 141.73 & 1599.97 & 0.91 & 0.46 & 1.13 \\
## 13 & 5.00 & 2.00 & 1199.87 & 2000.03 & 0.40 & 0.87 & 3.35 \\
## 14 & 5.00 & 3.00 & 851.77 & 2008.67 & 0.58 & 0.76 & 2.83 \\
## 15 & 5.00 & 4.00 & 547.43 & 2008.53 & 0.73 & 0.62 & 2.22 \\
## 16 & 5.00 & 5.00 & 330.17 & 2000.63 & 0.84 & 0.50 & 1.65 \\
## \hline
## \end{tabular}
## \end{table}
write.csv(x=df.ST, file="dfST.csv")
#========================================================
df.Report <- data.frame(Lambda, Mu,
LoDropped, mDropped, UpDropped,
LoAttended, mAttended, UpAttended,
LoRatioDropAtt, mRatioDropAtt, UpRatioDropAtt,
LoSystemTime, mSystemTime, UpSystemTime,
LoQueueLength, mQueueLength, UpQueueLength)
write.csv(x=df.Report, file="FullReport.csv")
print(df.Report)
## Lambda Mu LoDropped mDropped UpDropped LoAttended mAttended
## 1 2 2 119.423998 127.533333 135.642669 782.8909 793.5000
## 2 2 3 32.590372 37.800000 43.009628 782.8744 795.2667
## 3 2 4 10.428626 12.633333 14.838040 791.1944 802.0333
## 4 2 5 3.482465 4.466667 5.450868 792.7149 801.3333
## 5 3 2 426.814529 439.766667 452.718804 1186.5208 1198.9000
## 6 3 3 186.378855 197.133333 207.887812 1185.9429 1198.5000
## 7 3 4 78.064027 85.466667 92.869306 1186.3405 1199.4667
## 8 3 5 36.900486 41.200000 45.499514 1190.8074 1205.2667
## 9 4 2 805.934447 820.000000 834.065553 1590.3434 1604.8667
## 10 4 3 469.772146 484.233333 498.694521 1583.9229 1599.3333
## 11 4 4 252.928195 265.733333 278.538471 1586.9973 1601.8667
## 12 4 5 132.875773 141.733333 150.590894 1583.6683 1599.9667
## 13 5 2 1181.561823 1199.866667 1218.171510 1982.8823 2000.0333
## 14 5 3 835.611932 851.766667 867.921402 1991.9090 2008.6667
## 15 5 4 536.070626 547.433333 558.796041 1995.9022 2008.5333
## 16 5 5 315.028955 330.166667 345.304378 1984.0591 2000.6333
## UpAttended LoRatioDropAtt mRatioDropAtt UpRatioDropAtt LoSystemTime
## 1 804.1091 0.8310332 0.8397499 0.8484666 1.2212650
## 2 807.6589 0.9467258 0.9528720 0.9590183 0.6832575
## 3 812.8723 0.9817473 0.9843660 0.9869846 0.4429113
## 4 809.9518 0.9932285 0.9944395 0.9956505 0.3142553
## 5 1211.2792 0.6254225 0.6335238 0.6416251 1.1704663
## 6 1211.0571 0.8278789 0.8358512 0.8438235 0.8087987
## 7 1212.5929 0.9233223 0.9289775 0.9346326 0.5476195
## 8 1219.7259 0.9626330 0.9659496 0.9692662 0.3937613
## 9 1619.3899 0.4836858 0.4892271 0.4947685 1.0074115
## 10 1614.7437 0.6905452 0.6975237 0.7045023 0.8061286
## 11 1616.7360 0.8273566 0.8343836 0.8414106 0.6094499
## 12 1616.2651 0.9068009 0.9116600 0.9165191 0.4555435
## 13 2017.1844 0.3946002 0.4002397 0.4058791 0.8578103
## 14 2025.4243 0.5705536 0.5761290 0.5817044 0.7484428
## 15 2021.1645 0.7225486 0.7275063 0.7324640 0.6170405
## 16 2017.2076 0.8285777 0.8352105 0.8418434 0.4894178
## mSystemTime UpSystemTime LoQueueLength mQueueLength UpQueueLength
## 1 1.2443727 1.2674803 1.5820628 1.6281808 1.6742987
## 2 0.7023446 0.7214316 0.7192764 0.7698396 0.8204028
## 3 0.4519107 0.4609100 0.3824268 0.4073608 0.4322947
## 4 0.3206152 0.3269751 0.2278915 0.2421428 0.2563942
## 5 1.1841988 1.1979313 2.5853527 2.6144460 2.6435393
## 6 0.8229626 0.8371265 1.5947234 1.6362634 1.6778034
## 7 0.5595536 0.5714876 0.9413547 0.9824714 1.0235881
## 8 0.4021676 0.4105739 0.6074032 0.6351770 0.6629507
## 9 1.0182058 1.0290002 3.0896379 3.1074902 3.1253426
## 10 0.8154854 0.8248422 2.3115834 2.3460515 2.3805196
## 11 0.6194764 0.6295029 1.6122582 1.6496477 1.6870372
## 12 0.4631760 0.4708085 1.0928528 1.1277345 1.1626162
## 13 0.8656055 0.8734006 3.3315589 3.3474306 3.3633024
## 14 0.7557719 0.7631011 2.8118232 2.8328948 2.8539664
## 15 0.6234108 0.6297811 2.1887953 2.2152214 2.2416476
## 16 0.4950510 0.5006841 1.6132059 1.6466410 1.6800761
#=============================================================
#=============================================================
# Fitting confidence level of the interval
#=============================================================
for(j in 1:4){
for(n in 1:(N1*N2)){
lo<-t.test(MultiResult[n,j,],conf.level = 0.95)$conf.int[1]
cl<-0.95
h<-1e-3
while(lo<0){
cl<-cl-h
lo<-t.test(MultiResult[n,j,],conf.level = cl)$conf.int[1]
}
if(cl<0.95)
{ print("============================================");
print(" confidence interval must be revised! ");
print("============================================");
print(paste("n =",n,"; ","j =",j));print(paste("Conf.Level =",cl));
print(t.test(MultiResult[n,j,],conf.level =cl) );
#LoDropped[n] <- t.test(MultiResult[n,j,],conf.level = cl)$conf.int[1];
#UpDropped[n] <- t.test(MultiResult[n,j,],conf.level = cl)$conf.int[2]
}
}
}
#==============================================================
# For example, if it is reported the case n=4 & j=1,
# we shall then definite a new confidence level
#==============================================================
n<-4
j<-1
tsdf<-NoReplications-1
v<-var(MultiResult[n,j,])
s<-sqrt(v)
plot(function(x) dt((x-mDropped[n])/s, df = tsdf), -8, 10, ylim = c(0, 0.5),
main = "Non-central t - Density", yaxs = "i")
dts<-function(x) dt((x-mDropped[n])/s, df = tsdf)
area <- integrate(dts, LoDropped[n], 0)$value
# New confident interval
alpha<-0.95-area
print(paste("P{ 0 < mu","<",signif( UpDropped[n],digits = 5) ,"} =",signif(0.95-area,digits = 3)))
## [1] "P{ 0 < mu < 5.4509 } = 1.75"