Arrival rate is 4, Service rate is 0.5 and capacity of the server in 10

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=1.5;mu=0.4;b=5
# Generating data using simulation 
set.seed(2)
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=30 # 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(.5,2.5,len=100)
h<-seq(0.05,1,len=100)
z<-matrix(0,100,100)
for(i in 1:100){
  for(j in 1:100){
    z[i,j]=likelihood(g[i],h[j])}}
image(g,h,z,xaxs="i",yaxs="i",xlab="Arrival Rate",ylab="Service Rate",xlim=c(0.5,2.5),ylim=c(0.05,1))
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(1,.25),upper=c(2,1),popSize = 100,maxiter = 1000)
fig<-image(g,h,z,xaxs="i",yaxs="i",xlab="Arrival Rate",ylab="Service Rate",xlim=c(1,2),ylim=c(0.05,1))
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.481039 0.6549585

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=13/30,sd=0.25) }


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(0.2,0.5)
arrival_rate<-NULL;service_rate<-NULL
arrival_rate[1]=init[1];service_rate[1]<-init[2]
for(i in 2:200000){
  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 714 635103 3746 170
S 814 662112 3746 177

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 48 356328 28811 12.4
S 360 2771424 28811 96.2

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 48 481896 38415 12.5
S 148 1734782 38415 45.2

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 50 380130 28811 13.2
S 56 484694 28811 16.8

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 45 42435 3746 11.30
S 24 21872 3746 5.84

Autocorrelation plot and Trace plot

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

Histograms of the marginal Posterior distributions

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

Credible region

We have used Genetic algorithm to find credible region

modiA<-A[seq(from=1500,to=200000,by=250)]
modiS<-S[seq(from=1500,to=200000,by=250)]

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**9}
  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**10000000000000}
}


fit_credi<-function(a,b,c,d){
  (10**2/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/(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(1,.5,.5,.1),upper=c(2,1.5,1,.6),popSize = 100,maxiter = 5000)
plot(GA2)

GA2@solution

       x1       x2        x3        x4

[1,] 1.709918 1.499724 0.5002435 0.3419836

Posterior Predictive distributions

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

par(mfrow=c(1,2))
hist(arrival_post,freq=F,ylim=c(0,1.5),xlab="Inter arrival time",main="Histogram of Interarrival time")
lines(density(arrival_post),col=9)
curve(dexp(x,1.5),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,0.4),col=5,lwd=1.5,lty=2,add=T)

Chapter 5 R Codes

Kuntal Bakuli

20 April 2022