Appendix A
In this section, we have included the code to analyze data set 1, using the software R. The codes for the graphs in the data analysis are also plotted.
Appendix A.1. Data Sets
Data set 1
(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43, 3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10, 5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
Data set 2
(12, 15, 22, 24, 24, 32, 32, 33, 34, 38, 38, 43, 44, 48, 52, 53, 54, 54, 55, 56, 57, 58, 58,59, 60, 60, 60, 60, 61, 62, 63, 65, 65, 67, 68, 70, 70, 72, 73, 75, 76, 76, 81, 83, 84, 85, 87, 91, 95, 96,98, 99, 109, 110, 121, 127, 129, 131, 143, 146, 146, 175, 175, 211, 233, 258, 258, 263, 297, 341, 341, 376)
Data set 3
(0.96, 1.06, 1.09, 1.16, 1.19, 1.20, 1.32, 1.33, 1.40, 1.42, 1.46, 1.49, 1.51, 1.52, 1.54, 1.57, 1.59, 1.68, 1.70, 1.70, 1.76, 1.76, 1.77, 1.80, 1.81, 1.86, 1.89, 1.89, 1.94, 2.20, 2.20, 2.22, 2.36, 2.36, 2.39, 2.41, 2.45, 2.69, 2.71, 2.73, 2.77, 2.80, 2.83, 2.87, 2.94, 2.98, 3.03, 3.04, 3.19, 3.31, 3.57, 3.73, 4.17, 4.27, 4.30, 4.36, 4.45, 4.79, 4.85, 4.97, 5.26, 5.33, 5.53, 5.55, 5.91, 6.25, 6.31, 7.62, 7.84, 8.49, 8.63, 8.99, 9.94, 10.43, 10.86, 11.18)
Appendix A.2. Graphics for the PDF of S-ML Distribution
x= 0:10
f= function(x,p) #defining the pdf of S-ML model
{
(((pi/2)*(p/(1+p))*exp(-2*p*x))*(((1+p)*exp(p*x))+(2*p*x)-1)*sin((pi/2)*
(1+(exp(-p*x)*((p*x)/(1+p))))*exp(-p*x)))
}
curve(f(x,p= 5),col="yellow",xlab="x", ylim=c(0,1),ylab="pdf",lwd=2 )
curve(f(x,p= 15),col="pink", lwd=2, add= TRUE)
curve(f(x,p= 50),col="purple", lwd=2, add= TRUE)
curve(f(x,p=100),col="orange", lwd=2, add= TRUE)
legend("topright",legend=c(expression(paste(beta," = ",5)),
expression(paste(beta," = ",15)),
expression(paste(beta, " = ",50)),
expression(paste(beta, " = ",100))),
ncol=1, col=c("yellow","pink","purple", "orange"),
lwd=c(2,2,2,2), cex=c(1,1,1,1),text.width = 0.1, inset=0.011, bty ="n")
### In the same way, we can plot the pdf for other beta~values.
Appendix A.3. Parameter Estimate along with GOF Metrics
install.packages (c("EstimationTools", "MASS", "plyr" ))
library(EstimationTools)
library(MASS)
library(plyr)
# Data set 1
st = c(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43,
3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10,
5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
# S-ML distribution
dSml = function(x, p, log = FALSE) #log of pdf of S-ML model
{
n=count(x)
loglik <- (log(pi/2)+log(p)-log(1+p)-(2*p*x)+log((1+p)*exp(x*p)*(2*p*x)-1)+
log(sin((pi/2)*(1+exp(-p*x)*(p*x)/(1+p))*exp(-x*p))))
if ( log == FALSE)
density <- exp(loglik)
else density <- loglik
return(density)
}
theta <- maxlogL(x =st, dist = "dSml",start = 0.59)
summary(theta)
# S-Lindley~distribution
dSL = function(x, p, log = FALSE) #log of pdf of S-Lindley model
{
n=count(x)
loglik <- (log(pi/2)+log(p^2)+log(1+x)-log(1+p)-(p*x)+
log(sin((pi/2)*(1+((x*p)/(1+p)))*exp(-x*p))))
if ( log == FALSE)
density <- exp(loglik)
else density <- loglik
return(density)
}
theta <- maxlogL(x =st, dist = "dSL",start = 0.59)
summary(theta)
# SE distribution
dSE = function(x, p, log = FALSE) #log of pdf of SE model
{
n=count(x)
loglik <- (log(pi/2)+log(p)-(p*x)+log(sin((pi/2)*exp(-x*p))))
if ( log == FALSE)
density <- exp(loglik)
else density <- loglik
return(density)
}
theta <- maxlogL(x =st, dist = "dSE",start = 0.59)
summary(theta)
#IL distribution
dIL = function(x,p, log = FALSE)
{
n=count(x)
loglik <- (2*log(p))-log(1+p)-(p/x)+log(1+x)-(3*log(x))
if ( log == FALSE)
density <- exp(loglik)
else density <- loglik
return(density)
}
theta <- maxlogL(x =st, dist = "dIL",start = 0.65)
summary(theta)
#Exp distribution
dE = function(x,p, log = FALSE) #log of pdf of Expo model
{
n=count(x)
loglik <- log(dexp(x, p, log = FALSE))
if ( log == FALSE)
density <- exp(loglik)
else density <- loglik
return(density)
}
st = c(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43,
3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10,
5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
theta <- maxlogL(x =st, dist = "dE",start = 0.6)
summary(theta)
Appendix A.4. KS Test Statistic, p-Value and Other Test Statistics
install.packages ("goftest")
library(goftest)
y = c(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43,
3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10,
5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
# S-ML~CDF
pSml = function(x,p)
{
p = 0.14870
cos((pi/2)*(1+(exp(-p*x)*(x*p)/(1+p)))*exp(-p*x))
}
ks1=ad1=cvm1=NULL
ks1=ks.test(y,pSml)
ad1=ad.test(y,pSml)
cvm1=cvm.test(y,pSml)
result1=c(ks1$statistic,ks1$p.value,ad1$statistic,ad1$p.value,
cvm1$statistic,cvm1$p.value)
# S-Lindley~CDF
pSL = function(x,p)
{
p = 0.22660
cos((pi/2)*(1+((x*p)/(1+p)))*exp(-p*x))
}
ks1=ad1=cvm1=NULL
ks1=ks.test(y,pSL)
ad1=ad.test(y,pSL)
cvm1=cvm.test(y,pSL)
result1=c(ks1$statistic,ks1$p.value,ad1$statistic,ad1$p.value,
cvm1$statistic,cvm1$p.value)
# S-Expo CDF
pSe = function(x,p)
{
p = 0.10780
cos((pi/2)*exp(-p*x))
}
ks1=ad1=cvm1=NULL
ks1=ks.test(y,pSe)
ad1=ad.test(y,pSe)
cvm1=cvm.test(y,pSe)
result1=c(ks1$statistic,ks1$p.value,ad1$statistic,ad1$p.value,
cvm1$statistic,cvm1$p.value)
#IL distribution
pIL = function(x,p)
{
p = 4.9096
(1 + (p/((1+p)*x)))*exp(-p/x)
}
ks1=ad1=cvm1=NULL
ks1=ks.test(y,pIL)
ad1=ad.test(y,pIL)
cvm1=cvm.test(y,pIL)
result1=c(ks1$statistic,ks1$p.value,ad1$statistic,ad1$p.value,
cvm1$statistic,cvm1$p.value)
# Expo~CDF
pEx = function(x,p)
{
p = 0.18848
pexp(x,p)
}
ks1=ad1=cvm1=NULL
ks1=ks.test(y,pEx)
ad1=ad.test(y,pEx)
cvm1=cvm.test(y,pEx)
result1=c(ks1$statistic,ks1$p.value,ad1$statistic,ad1$p.value,
cvm1$statistic,cvm1$p.value)
Appendix A.5. Graphics—To Plot the EPDF for the First Data Set
x = c(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43,
3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10,
5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
hist(x,prob=T,main="Histogram and estimated PDFs",
col="pink", ylab = "PDF",ylim=c(0,0.05), bty ="n")
p = 0.14870 ## parameter estimate of S-ML model
curve(((pi/2)*(p/(1+p))*(exp(-2*p*x))*((1+p)*exp(p*x)+(2*x*p)-1)*
(sin((pi/2)*(1+exp(-p*x)*((x*p)/(1+p)))*exp(-x*p)))),col="blue",lwd=3,
add=T)
p = 0.22660 # parameter estimate of S-Lindley
curve(((pi/2)*((p^2)/(1+p))*(1+x)*exp(-p*x)*(sin((pi/2)*
(1+((x*p)/(1+p)))*exp(-x*p)))), col="green",lwd = 3, add=T)
p = 0.10780 # parameter estimate of S-Expo
curve(((pi/2)*(p*exp(-p*x))*(sin((pi/2)*exp(-x*p)))), col="orange",
lwd = 3, add=T)
p = 4.9096 # parameter estimate of IL
curve((((p^2)/(1+p))*((1+x)/x^3)*exp(-p/x)),col="red",lwd = 3, add=T)
p = 0.18848 # parameter estimate of Expo
curve(dexp(x,p), col="yellow", lwd = 3, add = T)
legend("topright",legend = c("S-ML","S-Lindley","S-Expo","IL","Expo"),
ncol = 1,
col= c("blue","green","orange","red","yellow"),lty =1,lwd=3,
text.width = 2.5 , inset= 0.00005, bty ="n")
Appendix A.6. Graphics—To Plot the ECDF for the First Data Set
y = c(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43,
3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10,
5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
plot(ecdf(y) , verticals=TRUE, main="Empirical and estimated CDFs",
ylab="CDF", xlab="x", bty ="n")
p = 0.14870 # parameter estimate of S-ML
curve((cos((pi/2)*(1+(exp(-p*x)*((x*p)/(1+p))))*exp(-p*x))),col="blue",
lwd=3, add=T)
p = 0.22660 # parameter estimate of S-Lindley
curve((cos((pi/2)*(1+((x*p)/(1+p)))*exp(-p*x))), col="green",lwd = 3,
add=T)
p = 0.010780 # parameter estimate of S-Expo
curve((cos((pi/2)*exp(-p*x))),col="orange",lwd = 3, add=T)
p = 4.9096 # parameter estimate of IL
curve(((1 + (p/((1+p)*x)))*exp(-p/x)),col="red",lwd = 3, add=T)
p = 0.18848 # parameter estimate of Expo
curve(pexp(x,p), col="yellow", lwd = 3, add = T)
legend("topleft",legend = c("S-ML","S-Lindley","S-Expo","IL","Expo"),
ncol = 1,
col= c("blue","green","orange","red","yellow"),lty =1,lwd=3,
text.width = 2.5 , inset= 0.00005, bty ="n")
#######
Appendix A.7. Graphics: Bar Plot and TTT Plot for First Data Set
###Bar plot
x = c(2.15, 2.20, 2.55, 2.56, 2.63, 2.74, 2.81, 2.90, 3.05, 3.41, 3.43,
3.43, 3.84, 4.16, 4.18, 4.36, 4.42, 4.51, 4.60, 4.61, 4.75, 5.03, 5.10,
5.44, 5.90, 5.96, 6.77, 7.82, 8.00, 8.16, 8.21, 8.72, 10.40, 13.20, 13.70)
boxplot(x,main = "Tumour size of lung cancer patients",
col = "orange", border="brown",horizontal = TRUE,notch = TRUE)
###TTT plot
install.packages ("AdequacyModel")
library(AdequacyModel)
TTT(x, lwd = 2, lty 2, col = "red", grid=FALSE)
############################