Chapter 5 R Codes- Additional Computation
Kuntal Bakuli
5 May 2022
The Arrival rate is 5, Service rate is 1 and Capacity of the server is 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=5;mu=1;b=10
# 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(1,10,len=100)
h<-seq(0.01,3,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(1,10),ylim=c(0.01,3))
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':
##
## deGA<-ga(type = "real-valued", fitness = function(x) (10^5)*likelihood(x[1],x[2]),lower=c(1,.01),upper=c(10,3),popSize = 100,maxiter = 1000)
fig<-image(g,h,z,xaxs="i",yaxs="i",xlab="Arrival Rate",ylab="Service Rate",xlim=c(1,10),ylim=c(0.01,3))
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,] 4.986404 0.8527041
R codes for the MCMC algorithm used in Case-1 of Section 3:
arrival_dens<-function(x){dgamma(x,shape=260,rate=60)}
service_dens<-function(x,y){
dnorm(x,mean=9/20,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(2,1.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.1chain1<-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 42 34560 3746 9.23
S 630 523593 3746 140.00
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 63 404019 28811 14.0
S 355 2784975 28811 96.7
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 60 662400 38415 17.2
S 156 1791855 38415 46.6
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 65 498693 28811 17.3
S 52 464893 28811 16.1
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 42 38064 3746 10.20
S 25 25045 3746 6.69
Autocorrelation plot and Trace plot
par(mfrow=c(2,2))
plot.ts(A,xlab="iterations",ylab="Arrival rate",ylim=c(0,10),col=9)
plot.ts(S,xlab="iterations",ylab="Service rate",ylim=c(0,3),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/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(5,2,1,.5),upper=c(7,5,3,1),popSize = 100,maxiter = 5000)
plot(GA2)
GA2@solution x1 x2 x3 x4[1,] 6.133661 4.998727 1.000559 0.6133661
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,2),xlab="Inter arrival time",main="Histogram of Interarrival time")
lines(density(arrival_post),col=9)
curve(dexp(x,5),col=6,lwd=1.5,lty=2,add=T)
hist(service_post,freq=F,ylim=c(0,1),xlab="Service time",main="Histogram of Service time")
lines(density(service_post),col=9)
curve(dexp(x,1),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,1.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
}
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,10),col=9)
plot.ts(S2,xlab="iterations",ylab="Service rate",ylim=c(0,3),col=10)
acf(A2[-c(1:1000)],lag.max = 500,main="For Arrival Rate")
acf(S2[-c(1:1000)],lag.max=500,main="For Service Rate")
Histograms of the marginal Posterior distributions
par(mfrow=c(1,2))
hist(A2[seq(from=1500,to=200000,by=250)],xlab=NULL,freq = FALSE,main=paste("Histogram of Arrival rate"))
lines(density(A2[seq(from=1500,to=200000,by=250)]),col=1,lwd=1.5)
hist(S2[seq(from=1500,to=200000,by=250)],xlab=NULL,freq=F,main=paste("Histogram of Service rate"))
lines(density(S2[seq(from=1500,to=200000,by=250)]),col=9,lwd=1.5)
Credible region
modiA2<-A2[seq(from=1000,to=200000,by=250)]
modiS2<-S2[seq(from=1000,to=200000,by=250)]
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(5,2,1,.5),upper=c(7,5,3,1),popSize = 100,maxiter = 5000)
plot(GA22)
GA22@solution x1 x2 x3 x4[1,] 5.94355 4.998625 1.001141 0.5943549
Posterior Predictive distribution
arrival_post<-rexp(length(A2[seq(from=1500,to=200000,by=250)]),c(A2[seq(from=1500,to=200000,by=250)]))
service_post<-rexp(length(S2[seq(from=1500,to=200000,by=250)]),c(S2[seq(from=1500,to=200000,by=250)]))
par(mfrow=c(1,2))
hist(arrival_post,freq=F,ylim=c(0,2),xlab="Inter arrival time",main="Histogram of Interarrival time")
lines(density(arrival_post),col=9)
curve(dexp(x,5),col=6,lwd=1.5,lty=2,add=T)
hist(service_post,freq=F,ylim=c(0,1),xlab="Service time",main="Histogram of Service time")
lines(density(service_post),col=9)
curve(dexp(x,1),col=5,lwd=1.5,lty=2,add=T)