The Arrival rate is 2, Service rate is 0.2 and Capacity of the server is 7

Data Generation

The following bunch of codes has been used to generate the overall data, data for the Subsection 2.3 and Section 3 is also generated.

# Fixing the values of the system parameters
lambda=2;mu=0.2;b=7
# Generating data using simulation 
set.seed(5)
arrival<-rexp(60*24*50,lambda)
cumA=cumsum(arrival)
service<-NULL
service[1]=cumA[b]+rexp(1,mu)
sys_size<-NULL
sys_size[1]=sum(as.numeric(cumA<service[1]))-b
for(j in 2:500){
  if(sys_size[j-1]<b){
    service[j]=cumA[j*b]+rexp(1,mu)
  }else{
    service[j]=service[j-1]+rexp(1,mu)
  }
  sys_size[j]=sum(as.numeric(cumA<service[j]))-j*b
}

# Generating data for the Subsection 2.3 and Section 3

x<-NULL;t<-NULL
n=90 # sample size
start=120 # starting time of sample collection
t[1]=start
x[1]=sum(as.numeric(cumA<t[1]))-b*sum(as.numeric(service<t[1]))
l=sum(as.numeric(service<120))
x[2:n]=sys_size[(l+1):(l+n-1)]
t[2:n]=service[(l+1):(l+n-1)]

Computation of Likelihood function

The likelihood function of Subsection 2.3, which is mentioned by the equation (2.4) has been computed using following bunch of codes.

##The likelihood function

## Stationary probability of the first observation

root<-polyroot(c(-(lambda/mu),seq(1,1,len=b)))
r1<-root[abs(Im(root))<10**(-10)&Re(root)<1&Re(root)>0]
r1<-Re(r1)
pi<-function(x){
  if(x[1]<b){p1=((1-r1**(x[1]-1))/b)}else{p1=(lambda/mu)*(1-r1)*r1**(x[1]-b)}
  return(p1)}


## Density function obtained from Semi Markov kernel

pdf_function<-function(t,j,k){if(j>=b){return(mu*((lambda*t)**(k+b-j))*exp(-(lambda+mu)*t)/factorial(k+b-j))}else{
  f1<-function(x){
    (mu*(lambda**(k+b-j))*(x**k)*((t-x)**(b-j-1))*exp(-lambda*t -mu*x))/(factorial(b-j-1) * factorial(k))
  }
  return(integrate(f1,lower=0,upper=t)$value)}
}


## Likelihood computation

L=pi(x[1])
for(i in 2:n){
  L=L*pdf_function(t[i]-t[i-1],x[i-1],x[i])
}


# Hence the Likelihood function is as follows;

likelihood<-function(lambda,mu){
  if(lambda< b*mu && lambda>0){
    root<-polyroot(c(-(lambda/mu),seq(1,1,len=b)))
    r1<-root[abs(Im(root))<10**(-10)&Re(root)<1&Re(root)>0]
    r1<-Re(r1)
    pi<-function(x){
      if(x[1]<b){p1=((1-r1**(x[1]-1))/b)}else{p1=(lambda/mu)*(1-r1)*r1**(x[1]-b)}
      return(p1)}
    pdf_function<-function(t,j,k){if(j>=b){return(mu*((lambda*t)**(k+b-j))*exp(-(lambda+mu)*t)/factorial(k+b-j))}else{
      f1<-function(x){
        (mu*(lambda**(k+b-j))*(x**k)*((t-x)**(b-j-1))*exp(-lambda*t -mu*x))/(factorial(b-j-1) * factorial(k))
      }
      return(integrate(f1,lower=0,upper=t)$value)}
    }
    L=pi(x[1])
    for(i in 2:n){
      L=L*pdf_function(t[i]-t[i-1],x[i-1],x[i])}
    return((L)*as.numeric(identity(lambda<b*mu)))}else{return(0)}
}

Contour plot of the Likelihood function

g<-seq(0.01,3,len=500)
h<-seq(0.01,.5,len=500)
z<-matrix(0,500,500)
for(i in 1:500){
  for(j in 1:500){
    z[i,j]=likelihood(g[i],h[j])}}
image(g,h,z,xaxs="i",yaxs="i",xlab="Arrival Rate",ylab="Service Rate",xlim=c(1.5,2.5),ylim=c(0.15,.35))
contour(g,h,z,add=T,xlim=range(x,finite=TRUE),ylim=range(y,finite=TRUE))

optimization of likelihood function

We have used genetic algorithm to optimize the likelihood function as follows.

library(GA)

## Loading required package: foreach

## Loading required package: iterators

## Package 'GA' version 3.2.1
## Type 'citation("GA")' for citing this R package in publications.

## 
## Attaching package: 'GA'

## The following object is masked from 'package:utils':
## 
##     de

GA<-ga(type = "real-valued", fitness = function(x) (10^5)*likelihood(x[1],x[2]),lower=c(.5,.01),upper=c(5,1),popSize = 100,maxiter = 3000)
fig<-image(g,h,z,xaxs="i",yaxs="i",xlab="Arrival Rate",ylab="Service Rate",xlim=c(1.5,2.5),ylim=c(.15,.35))
contour(g,h,z,add=T,xlim=range(x,finite=TRUE),ylim=range(y,finite=TRUE),nlevels = 20)
points(GA@solution,col=5,pch=2,lwd=2)

plot(GA)

GA@solution

       x1        x2

[1,] 1.964106 0.2847781

R codes for the MCMC algorithm used in Case-1 of Section 3:

arrival_dens<-function(x){dgamma(x,shape=120,rate=60)}
service_dens<-function(x,y){
  dnorm(x,mean=0.33,sd=1) }


g<-function(lambda,mu){
  likelihood(lambda,mu)*arrival_dens(lambda)*service_dens(mu,lambda)
}             # g function




update_mu<-function(mu_old){
  g_old1<-g(lambda_old,mu_old)
  e1<-(mu_old-(lambda_old)/b)/2
  mu_prop<-runif(1,min=mu_old-e1,max=mu_old+e1)
  g_prop1<-g(lambda_old,mu_prop)
  if(g_prop1>=g_old1){mu_old=mu_prop}else{
    log_alpha1=log(g_prop1)-log(g_old1)
    alpha1=exp(log_alpha1)
    a<-runif(1)
    if(a<alpha1){mu_old=mu_prop}
  }
  return(mu_old)
}



update_lambda<-function(lambda_old){
  g_old2<-g(lambda_old,mu_new)
  e2<-min(lambda_old,(b*mu_new-lambda_old))/1.25
  lambda_prop<-runif(1,min=lambda_old-e2,max=lambda_old+e2)
  g_prop2<-g(lambda_prop,mu_new)
  if(g_prop2>=g_old2){lambda_old=lambda_prop}else{
    log_alpha2=log(g_prop2)-log(g_old2)
    alpha2=exp(log_alpha2)
    a<-runif(1)
    if(a<alpha2){lambda_old=lambda_prop}
  }
  return(lambda_old)
}


init<-c(1,0.3)
arrival_rate<-NULL;service_rate<-NULL
arrival_rate[1]=init[1];service_rate[1]<-init[2]
for(i in 2:500000){
  mu_old=service_rate[i-1];lambda_old=arrival_rate[i-1]
  mu_new<-update_mu(mu_old)
  lambda_new<-update_lambda(lambda_old)
  arrival_rate[i]=lambda_new
  service_rate[i]=mu_new
}
A<-as.ts(arrival_rate);S<-as.ts(service_rate)

Raftery and Lewi’s daignostic

library(coda)

## Warning: package 'coda' was built under R version 4.1.1

chain1<-as.mcmc(cbind.data.frame(A,S))
raftery.diag(chain1,q=0.025,s=0.95,r=0.005,converge.eps = 0.0001)

Quantile (q) = 0.025 Accuracy (r) = +/- 0.005 Probability (s) = 0.95

Burn-in Total Lower bound Dependence (M) (N) (Nmin) factor (I) A 4320 3319360 3746 886
S 5202 4256256 3746 1140

raftery.diag(chain1,q=0.25,s=0.95,r=0.005,converge.eps = 0.0001)

Quantile (q) = 0.25 Accuracy (r) = +/- 0.005 Probability (s) = 0.95

Burn-in Total Lower bound Dependence (M) (N) (Nmin) factor (I) A 4920 31183206 28811 1080
S 5187 33339033 28811 1160

raftery.diag(chain1,q=0.5,s=0.95,r=0.005,converge.eps = 0.0001)

Quantile (q) = 0.5 Accuracy (r) = +/- 0.005 Probability (s) = 0.95

Burn-in Total Lower bound Dependence (M) (N) (Nmin) factor (I) A 4862 42792178 38415 1110
S 4074 37207842 38415 969

raftery.diag(chain1,q=0.75,s=0.95,r=0.005,converge.eps = 0.0001)

Quantile (q) = 0.75 Accuracy (r) = +/- 0.005 Probability (s) = 0.95

Burn-in Total Lower bound Dependence (M) (N) (Nmin) factor (I) A 3196 20754636 28811 720
S 3240 19935882 28811 692

raftery.diag(chain1,q=0.975,s=0.95,r=0.005,converge.eps = 0.0001)

Quantile (q) = 0.975 Accuracy (r) = +/- 0.005 Probability (s) = 0.95

Burn-in Total Lower bound Dependence (M) (N) (Nmin) factor (I) A 1092 942552 3746 252
S 1024 795520 3746 212

Autocorrelation plot and Trace plot

par(mfrow=c(2,2))
plot.ts(A,xlab="iterations",ylab="Arrival rate",ylim=c(0,4),col=9)
plot.ts(S,xlab="iterations",ylab="Service rate",ylim=c(0,.4),col=10)
acf(A[-c(1:1000)],lag.max = 1000,main="For Arrival Rate")
acf(S[-c(1:1000)],lag.max=1000,main="For Service Rate")

Histograms of the marginal Posterior distributions

par(mfrow=c(1,2))
hist(A[seq(from=1000,to=500000,by=1000)],xlab=NULL,freq = FALSE,main=paste("Histogram of Arrival rate"))
lines(density(A[seq(from=1000,to=500000,by=1000)]),col=1,lwd=1.5)
hist(S[seq(from=1000,to=500000,by=1000)],xlab=NULL,freq=F,main=paste("Histogram of Service rate"))
lines(density(S[seq(from=1000,to=500000,by=1000)]),col=9,lwd=1.5)

Credible region

We have used Genetic algorithm to find credible region

modiA<-A[seq(from=1000,to=500000,by=1000)]
modiS<-S[seq(from=1000,to=500000,by=1000)]

credible<-function(lamhi,lamlo,muhi,mulo){
  vec=as.numeric(modiS>mulo)*as.numeric(modiS<muhi)*as.numeric(modiA>lamlo)*as.numeric(modiA<lamhi)
  return(sum(vec)/length(modiA))
}
# Area of the restricted parameter space
Area<-function(l2,l1,mu2,mu1){
  if(l2>l1 & mu2>mu1){
    if(l1>b*mu2){10**5}
    else if(l1>b*mu1 & l1<b*mu2 & l2>b*mu2){(.5/b)*(b*mu2-l1)^2}
    else if(l1 <b*mu1 & l2 > b*mu1 & l2>b*mu2){(mu2-mu1)*((mu2+mu1)*b-2*l1)*.5}
    else if(l1> b*mu1 & l2 <b*mu2){(.5/b)*((b*mu2)**2-l1**2-(b*mu2-l2)**2)}
    else if(l1<b*mu1 & l2>b*mu1 & l2<b*mu2){(l2-l1)*(mu2-mu1)*(.5/b)*(l2-b*mu1)**2}
    else if(l2<b*mu1){(l2-l1)*(mu2-mu1)}}else{10**1000}
}


fit_credi<-function(a,b,c,d){
  (10/abs(credible(a,b,c,d)-.95))*as.numeric(I(b>0))*as.numeric(I(a>b))*as.numeric(I(d>0))*as.numeric(I(c>d))+1/(10*Area(a,b,c,d))
}
library(GA)
GA2<-ga(type = "real-valued", fitness = function(x) fit_credi(x[1],x[2],x[3],x[4]),lower=c(2,1,0.2,0.05),upper=c(3,1.9999,.3,.19999),popSize = 50,maxiter = 5000)
plot(GA2)

GA2@solution

       x1       x2        x3       x4

[1,] 2.478775 1.735349 0.2919946 0.137777

Posterior Predictive distributions

arrival_post<-rexp(length(A[seq(from=1500,to=500000,by=1000)]),c(A[seq(from=1500,to=500000,by=1000)]))
service_post<-rexp(length(S[seq(from=1500,to=500000,by=1000)]),c(S[seq(from=1500,to=500000,by=1000)]))

par(mfrow=c(1,2))
hist(arrival_post,freq=F,ylim=c(0,5),xlab="Inter arrival time",main="Histogram of Interarrival time")
lines(density(arrival_post),col=9)
curve(dexp(x,2),col=6,lwd=1.5,lty=2,add=T)
hist(service_post,freq=F,ylim=c(0,.5),xlab="Service time",main="Histogram of Service time")
lines(density(service_post),col=9)
curve(dexp(x,.2),col=5,lwd=1.5,lty=2,add=T)

R codes for the MCMC algorithm used in Case-2 of Section 3

a1=sum(as.numeric(cumA<120))-sum(as.numeric(cumA<90))
a2=sum(as.numeric(service<120))-sum(as.numeric(service<90))
# prior densities
arrival_dens<-function(x){dgamma(x,shape=a1,rate=30)}
service_dens<-function(x,y){
  dgamma(x,shape=a2,rate=30)
}
# g function
g<-function(lambda,mu){
  likelihood(lambda,mu)*arrival_dens(lambda)*service_dens(mu,lambda)
}         
init<-c(2,.3)
arrival_rate<-NULL;service_rate<-NULL
arrival_rate[1]=init[1];service_rate[1]<-init[2]
for(i in 2:500000){
  mu_old=service_rate[i-1];lambda_old=arrival_rate[i-1]
  mu_new<-update_mu(mu_old)
  lambda_new<-update_lambda(lambda_old)
  arrival_rate[i]=lambda_new
  service_rate[i]=mu_new
}
A2<-as.ts(arrival_rate);S2<-as.ts(service_rate)

Raftery and Lewi’s daignostic

library(coda)
chain2<-as.mcmc(cbind.data.frame(A2,S2))
plot(chain2)

raftery.diag(chain2,q=0.025,s=0.95,r=0.005,converge.eps = 0.0001)
raftery.diag(chain2,q=0.25,s=0.95,r=0.005,converge.eps = 0.0001)
raftery.diag(chain2,q=0.5,s=0.95,r=0.005,converge.eps = 0.0001)
raftery.diag(chain2,q=0.75,s=0.95,r=0.005,converge.eps = 0.0001)
raftery.diag(chain2,q=0.975,s=0.95,r=0.005,converge.eps = 0.0001)

Autocorrelation plot and Trace plot

par(mfrow=c(2,2))
plot.ts(A2,xlab="iterations",ylab="Arrival rate",ylim=c(0,4),col=9)
plot.ts(S2,xlab="iterations",ylab="Service rate",ylim=c(0,.4),col=10)
acf(A2[-c(1:1000)],lag.max = 1000,main="For Arrival Rate")
acf(S2[-c(1:1000)],lag.max=1000,main="For Service Rate")

Histograms of the marginal Posterior distributions

par(mfrow=c(1,2))
hist(A2[seq(from=1500,to=500000,by=1200)],xlab=NULL,freq = FALSE,main=paste("Histogram of Arrival rate"))
lines(density(A2[seq(from=1500,to=500000,by=1200)]),col=1,lwd=1.5)
hist(S2[seq(from=1500,to=500000,by=1200)],xlab=NULL,freq=F,main=paste("Histogram of Service rate"))
lines(density(S2[seq(from=1500,to=500000,by=1200)]),col=9,lwd=1.5)

Credible region

modiA2<-A2[seq(from=1000,to=500000,by=1200)]
modiS2<-S2[seq(from=1000,to=500000,by=1200)]

credible<-function(lamhi,lamlo,muhi,mulo){
  vec=as.numeric(modiS2>mulo)*as.numeric(modiS2<muhi)*as.numeric(modiA2>lamlo)*as.numeric(modiA2<lamhi)
  return(sum(vec)/length(modiA2))
}

library(GA)
GA22<-ga(type = "real-valued", fitness = function(x) fit_credi(x[1],x[2],x[3],x[4]),lower=c(1.7,1,.2,0.15),upper=c(3,1.9,0.3,0.25),popSize = 100,maxiter = 5000)
plot(GA22)

GA22@solution

        x1       x2        x3        x4

[1,] 2.365327 1.587313 0.2267592 0.1935528 [2,] 2.505374 1.587313 0.2267592 0.1931404 [3,] 2.396349 1.587313 0.2267592 0.1909488 [4,] 2.366593 1.587313 0.2267592 0.1911565 [5,] 2.536722 1.587313 0.2267592 0.1903096 [6,] 2.382117 1.587313 0.2267592 0.1915198 [7,] 2.379552 1.587313 0.2267592 0.1914036 [8,] 2.280732 1.587313 0.2267592 0.1949583 [9,] 2.363345 1.587313 0.2267592 0.1900399 [10,] 2.483611 1.587313 0.2267592 0.1933296 [11,] 2.410092 1.587313 0.2267592 0.1895195 [12,] 2.363074 1.587313 0.2267592 0.1901943 [13,] 2.421948 1.587313 0.2267592 0.1907171 [14,] 2.425752 1.587313 0.2267592 0.1888217 [15,] 2.392647 1.587313 0.2267592 0.1902179 [16,] 2.412475 1.587313 0.2267592 0.1906829 [17,] 2.400424 1.587313 0.2267592 0.1905707 [18,] 2.370712 1.587313 0.2267592 0.1903219 [19,] 2.393429 1.587313 0.2267592 0.1897448 [20,] 2.414601 1.587313 0.2267592 0.1915158 [21,] 2.410444 1.587313 0.2267592 0.1813989 [22,] 2.354929 1.587313 0.2267592 0.1915782 [23,] 2.410965 1.587313 0.2267592 0.1895159 [24,] 2.386016 1.587313 0.2267592 0.1928843 [25,] 2.410062 1.587313 0.2267592 0.1902322 [26,] 2.130886 1.587313 0.2267592 0.1914714 [27,] 2.375848 1.587313 0.2267592 0.1908565 [28,] 2.418524 1.587313 0.2267592 0.1897943 [29,] 2.398937 1.587313 0.2267592 0.1906445 [30,] 2.388334 1.587313 0.2267592 0.1911771 [31,] 2.484433 1.587313 0.2267592 0.1914309 [32,] 2.409576 1.587313 0.2267592 0.1909450 [33,] 2.430866 1.587313 0.2267592 0.1918028 [34,] 2.492169 1.587313 0.2267592 0.1908691 [35,] 2.323272 1.587313 0.2267592 0.1951438 [36,] 2.337947 1.587313 0.2267592 0.1946341 [37,] 2.404636 1.587313 0.2267592 0.1877303 [38,] 2.383146 1.587313 0.2267592 0.1914288 [39,] 2.279863 1.587313 0.2267592 0.1876009 [40,] 2.320965 1.587313 0.2267592 0.1924548 [41,] 2.361712 1.587313 0.2267592 0.1951467 [42,] 2.401227 1.587313 0.2267592 0.1854992 [43,] 2.404088 1.587313 0.2267592 0.1903020 [44,] 2.365612 1.587313 0.2267592 0.1912698 [45,] 2.397547 1.587313 0.2267592 0.1916616 [46,] 2.390822 1.587313 0.2267592 0.1897233 [47,] 2.368740 1.587313 0.2267592 0.1927046 [48,] 2.334114 1.587313 0.2267592 0.1892824 [49,] 2.343067 1.587313 0.2267592 0.1923385 [50,] 2.379642 1.587313 0.2267592 0.1914870 [51,] 2.421511 1.587313 0.2267592 0.1905601 [52,] 2.383372 1.587313 0.2267592 0.1903331 [53,] 2.410282 1.587313 0.2267592 0.1907898 [54,] 2.363328 1.587313 0.2267592 0.1906417 [55,] 2.423047 1.587313 0.2267592 0.1912058 [56,] 2.407575 1.587313 0.2267592 0.1915185 [57,] 2.396086 1.587313 0.2267592 0.1910991 [58,] 2.413253 1.587313 0.2267592 0.1906393 [59,] 2.502996 1.587313 0.2267592 0.1916614 [60,] 2.434740 1.587313 0.2267592 0.1914666 [61,] 2.410737 1.587313 0.2267592 0.1902620 [62,] 2.341499 1.587313 0.2267592 0.1890917 [63,] 2.382804 1.587313 0.2267592 0.2158170 [64,] 2.348822 1.587313 0.2267592 0.1911050 [65,] 2.371394 1.587313 0.2267592 0.1925267 [66,] 2.334304 1.587313 0.2267592 0.1930520 [67,] 2.395539 1.587313 0.2267592 0.2011779 [68,] 2.408329 1.587313 0.2267592 0.1908302 [69,] 2.429034 1.587313 0.2267592 0.1872870 [70,] 2.404852 1.587313 0.2267592 0.1929044 [71,] 2.381403 1.587313 0.2267592 0.1923255 [72,] 2.320916 1.587313 0.2267592 0.1892252 [73,] 2.504407 1.587313 0.2267592 0.1909727 [74,] 2.424092 1.587313 0.2267592 0.1906083 [75,] 2.379849 1.587313 0.2267592 0.1944446 [76,] 2.418460 1.587313 0.2267592 0.1910142 [77,] 2.380721 1.587313 0.2267592 0.1899432 [78,] 2.316251 1.587313 0.2267592 0.1925607 [79,] 2.391862 1.587313 0.2267592 0.1920746 [80,] 2.393437 1.587313 0.2267592 0.1910213 [81,] 2.395270 1.587313 0.2267592 0.1885845 [82,] 2.371260 1.587313 0.2267592 0.1944961 [83,] 2.050310 1.587313 0.2267592 0.1910712 [84,] 2.386984 1.587313 0.2267592 0.1905123 [85,] 2.316887 1.587313 0.2267592 0.1960256 [86,] 2.330440 1.587313 0.2267592 0.1893787 [87,] 2.386700 1.587313 0.2267592 0.1910799 [88,] 2.426943 1.587313 0.2267592 0.1908816 [89,] 2.391641 1.587313 0.2267592 0.1935499

Posterior Predictive distribution

arrival_post<-rexp(length(A2[seq(from=1500,to=500000,by=1200)]),c(A2[seq(from=1500,to=500000,by=1200)]))
service_post<-rexp(length(S2[seq(from=1500,to=500000,by=1200)]),c(S2[seq(from=1500,to=500000,by=1200)]))

par(mfrow=c(1,2))
hist(arrival_post,freq=F,ylim=c(0,4),xlab="Inter arrival time",main="Histogram of Interarrival time")
lines(density(arrival_post),col=9)
curve(dexp(x,2),col=6,lwd=1.5,lty=2,add=T)
hist(service_post,freq=F,ylim=c(0,.5),xlab="Service time",main="Histogram of Service time")
lines(density(service_post),col=9)
curve(dexp(x,.2),col=5,lwd=1.5,lty=2,add=T)

Add_comp 3

Kuntal Bakuli

5/7/2022

The Arrival rate is 2, Service rate is 0.2 and Capacity of the server is 7

Data Generation

Computation of Likelihood function

Contour plot of the Likelihood function

optimization of likelihood function

R codes for the MCMC algorithm used in Case-1 of Section 3:

Raftery and Lewi’s daignostic

Autocorrelation plot and Trace plot

Histograms of the marginal Posterior distributions

Credible region

Posterior Predictive distributions

R codes for the MCMC algorithm used in Case-2 of Section 3

Raftery and Lewi’s daignostic

Autocorrelation plot and Trace plot

Histograms of the marginal Posterior distributions

Credible region

Posterior Predictive distribution