install.packages("copula")
install.packages("VineCopula")
install.packages("DEoptim")
install.packages("foreach")
install.packages("doParallel")
install.packages("parallel")
install.packages("tictoc")
install.packages("statmod")
library(statmod)

library(LaplacesDemon)
library(foreach)
library(doParallel)
library(parallel)
library(tictoc)
library(copula)
library(VineCopula)
library(DEoptim)



#================= system reliability ========================================

fx=function(x,t,a1,b1){
  b1^(-exp(a1*t))*gamma(exp(a1*t))^(-1)*integrate(function(z) exp(-z/b1)*z^(exp(a1*t)-1),0,x)$value
}


fy=function(y,t,a2,b2){
  b2^(-exp(a2*t))*gamma(exp(a2*t))^(-1)*integrate(function(z) 
    exp(-z/b2)*z^(exp(a2*t)-1),0,y)$value
}



Fg=function(x,y,t,a1,a2,b1,b2,p){
  exp( -((-log(fx(x,t,a1,b1)))^p+(-log(fy(y,t,a2,b2)))^p)^(1/p) )
}
#=====second and easy method for writing copula=================

Fg2=function(x,y,t,a1,a2,b1,b2,p){
  pCopula(c(fx(x,t,a1,b1),fy(y,t,a2,b2)),gumbelCopula(p))
}

Fg(0.5,0.7,0.08,32,27,0.03,0.02,4)
Fg2(0.5,0.7,0.08,32,27,0.03,0.02,4)


Fbg=function(x,y,t,a1,a2,b1,b2,p){
  Fg(x,y,t,a1,a2,b1,b2,p)-fy(y,t,a2,b2)+1-fx(x,t,a1,b1)
}

n=10
k1=2
k2=6

Rg=function(x,y,t,a1,a2,b1,b2,p){
  sum(sapply(k1:n,function(i)sum(sapply(k2:n,function(j) 
    sum(sapply(max(0,n-i-j):min(n-i,n-j),function(k)   
      (factorial(n)/(factorial(k)*factorial(n-k-j)*factorial(n-k-i)*
                       factorial(i+j+k-n)))
      *(Fbg(x,y,t,a1,a2,b1,b2,p))^k*(fx(x,t,a1,b1)-Fg(x,y,t,a1,a2,b1,b2,p))^(n-k-j)
      *(fy(y,t,a2,b2)-Fg(x,y,t,a1,a2,b1,b2,p))^(n-k-i)*
        (Fg(x,y,t,a1,a2,b1,b2,p))^(i+j+k-n) 
      ))
    ))
    ))
}

#qnorm(0.9,0,1)

t = c(0.07,0.08,0.09, 0.095,0.097,0.098,0.1, 0.101,0.102,0.103,0.105,0.106,0.107,0.108,0.109,0.11,0.1105,0.111,0.1115,0.112,0.1125,0.113,0.1135,0.114,0.115,0.116,0.117,0.1175,0.118,0.1185,0.119,0.1195, 0.12,0.123,0.125,0.127,0.130,0.135,0.14)


frg2=function(t,a1,a2,b1,b2,p){Rg(0.7,0.5,t,a1,a2,b1,b2,p)}
length(t)

#=======confidence intervals (CIs) for parameters and system reliability======
#based on bootstrap method====================================================


k1=2
k2=6

ff=function(n){
  a1=30.1
  a2=27.06
  b1=0.028
  b2=0.025
  p=0.93
  u=array(0,c(12,2,n))
  DX=matrix(0,12,n)
  DY=matrix(0,12,n)
  t=seq(0,0.2,0.01)
  
  for(j in 1:n){
    for(i in 1:12){
      library(copula)
      library(VineCopula)
      library(DEoptim)
      u[i,,j]=rCopula(1,normalCopula(p))
      DX[i,j]=qgamma(u[i,1,j] ,shape=exp(t[i+1]*a1)-exp(t[i]*a1), scale = b1)
      DY[i,j]=qgamma(u[i,2,j] ,shape=exp(t[i+1]*a2)-exp(t[i]*a2), scale = b2)
      
    }
  }
  
  X=matrix(0,12,n)
  Y=matrix(0,12,n)
  for(j in 1:n){
    for(i in 1:12){
      X[i,j]=sum(DX[1:i,j])
      Y[i,j]=sum(DY[1:i,j])
    }
  }
  
  
  d1=0.7
  d2=0.5
  tX=array(12,n)
  tY=array(12,n)
  for(j in 1:n){
    if(X[12,j]<d1) tX[j]=12 else tX[j]=which.max(X[,j]>d1)
    if(Y[12,j]<d2) tY[j]=12 else tY[j]=which.max(Y[,j]>d2)
  }
  
  m1=sort(tX)
  m2=sort(tY)
  jst=min(m1[n-k1+1],m2[n-k2+1])
  
  L=function(tt) -sum(sapply(1:n,function(j) sum(sapply(1:min(tX[j],tY[j],jst),function(i) log(dCopula(u[i,,j], normalCopula(tt))) ))))
  
  
  Lx=function(aa){
    a=aa[1]
    b=aa[2]
    sum(sapply(1:n,function(j) sum(sapply(1:min(tX[j],jst),function(i)  (exp(a*t[i+1])-exp(a*t[i]))*log(b)+DX[i,j]/b-(exp(a*t[i+1])-exp(a*t[i])-1)*log(DX[i,j])+log(gamma(exp(a*t[i+1])-exp(a*t[i]))) ))))
  }
  
  
  
  Ly=function(cc){
    c=cc[1]
    d=cc[2]
    sum(sapply(1:n,function(j) sum(sapply(1:min(tY[j],jst),function(i)  (exp(c*t[i+1])-exp(c*t[i]))*log(d)+DY[i,j]/d-(exp(c*t[i+1])-exp(c*t[i])-1)*log(DY[i,j])+log(gamma(exp(c*t[i+1])-exp(c*t[i]))) ))))
  }
  
  
  theta=DEoptim(L,lower=0.002,upper = 1,control=list(itermax=10))$optim$bestmem
  theta1=DEoptim(Lx,lower=c(16,0.01),upper=c(35,0.1),control=list(itermax=10))$optim$bestmem
  theta2=DEoptim(Ly,lower=c(17,0.01),upper=c(32,0.1),control=list(itermax=10))$optim$bestmem
  return(c(theta,theta1,theta2))
}



p.quan=a1.quan=a2.quan=b1.quan=b2.quan=c()
A<- matrix(0,150,5)
registerDoParallel(makeCluster(detectCores()))
tic()
A=foreach(j=1:150,.combine='rbind') %dopar% c(ff(10))
B<- matrix(0,150,39)
for(i in 1:150){
  for(j in 1:39){
    B[i,j]=frg2(t[j],A[i,2],A[i,4],A[i,3],A[i,5],A[i,1])  
  }
}
toc()
stopImplicitCluster()
stopCluster(makeCluster(detectCores()))


reli=matrix(0,39,4)
p.quan=quantile(A[,1],  probs = c(2.5,5,95,97.5)/100)
a1.quan=quantile(A[,2],  probs = c(2.5,5,95,97.5)/100)
b1.quan=quantile(A[,3],  probs = c(2.5,5,95,97.5)/100)
a2.quan=quantile(A[,4],  probs = c(2.5,5,95,97.5)/100)
b2.quan=quantile(A[,5],  probs = c(2.5,5,95,97.5)/100)

reli[1,]=quantile(B[,1],  probs = c(2.5,5,95,97.5)/100)
reli[2,]=quantile(B[,2],  probs = c(2.5,5,95,97.5)/100)
reli[3,]=quantile(B[,3],  probs = c(2.5,5,95,97.5)/100)
reli[4,]=quantile(B[,4],  probs = c(2.5,5,95,97.5)/100)
reli[5,]=quantile(B[,5],  probs = c(2.5,5,95,97.5)/100)
reli[6,]=quantile(B[,6],  probs = c(2.5,5,95,97.5)/100)
reli[7,]=quantile(B[,7],  probs = c(2.5,5,95,97.5)/100)
reli[8,]=quantile(B[,8],  probs = c(2.5,5,95,97.5)/100)
reli[9,]=quantile(B[,9],  probs = c(2.5,5,95,97.5)/100)
reli[10,]=quantile(B[,10],  probs = c(2.5,5,95,97.5)/100)
reli[11,]=quantile(B[,11],  probs = c(2.5,5,95,97.5)/100)
reli[12,]=quantile(B[,12],  probs = c(2.5,5,95,97.5)/100)
reli[13,]=quantile(B[,13],  probs = c(2.5,5,95,97.5)/100)
reli[14,]=quantile(B[,14],  probs = c(2.5,5,95,97.5)/100)
reli[15,]=quantile(B[,15],  probs = c(2.5,5,95,97.5)/100)
reli[16,]=quantile(B[,16],  probs = c(2.5,5,95,97.5)/100)
reli[17,]=quantile(B[,17],  probs = c(2.5,5,95,97.5)/100)
reli[18,]=quantile(B[,18],  probs = c(2.5,5,95,97.5)/100)
reli[19,]=quantile(B[,19],  probs = c(2.5,5,95,97.5)/100)
reli[20,]=quantile(B[,20],  probs = c(2.5,5,95,97.5)/100)
reli[21,]=quantile(B[,21],  probs = c(2.5,5,95,97.5)/100)
reli[22,]=quantile(B[,22],  probs = c(2.5,5,95,97.5)/100)
reli[23,]=quantile(B[,23],  probs = c(2.5,5,95,97.5)/100)
reli[24,]=quantile(B[,24],  probs = c(2.5,5,95,97.5)/100)
reli[25,]=quantile(B[,25],  probs = c(2.5,5,95,97.5)/100)
reli[26,]=quantile(B[,26],  probs = c(2.5,5,95,97.5)/100)
reli[27,]=quantile(B[,27],  probs = c(2.5,5,95,97.5)/100)
reli[28,]=quantile(B[,28],  probs = c(2.5,5,95,97.5)/100)
reli[29,]=quantile(B[,29],  probs = c(2.5,5,95,97.5)/100)
reli[30,]=quantile(B[,30],  probs = c(2.5,5,95,97.5)/100)
reli[31,]=quantile(B[,31],  probs = c(2.5,5,95,97.5)/100)
reli[32,]=quantile(B[,32],  probs = c(2.5,5,95,97.5)/100)
reli[33,]=quantile(B[,33],  probs = c(2.5,5,95,97.5)/100)
reli[34,]=quantile(B[,34],  probs = c(2.5,5,95,97.5)/100)
reli[35,]=quantile(B[,35],  probs = c(2.5,5,95,97.5)/100)
reli[36,]=quantile(B[,36],  probs = c(2.5,5,95,97.5)/100)
reli[37,]=quantile(B[,37],  probs = c(2.5,5,95,97.5)/100)
reli[38,]=quantile(B[,38],  probs = c(2.5,5,95,97.5)/100)
reli[39,]=quantile(B[,39],  probs = c(2.5,5,95,97.5)/100)





write(rbind(reli),"C:/Users/NP/Desktop/Result9/ reliability150.txt")



#=====================   plot   ===============================================

f=function(t){Rg(0.7,0.5,t,30.1,27.06,0.028,0.025,3.1)}
r1=c()
for(i in 1:39){
 r1[i]=f(t[i])
}

plot(t[-1],r1[-1],col=1,xlab="time (in million cycles)", ylab=" (2, 6)-out-of-10 System Reliability",type="l")
lines(t[-1],reli[-1,1],col=1,lty=5)
lines(t[-1],reli[-1,4],col=1,lty=5)
grid( ny = 12, col = "darkgray", lty = 1,
     lwd = 0.01, equilogs = TRUE)
<pre>An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler-&gt;renderException()
#3 [internal function]: yii\base\ErrorHandler-&gt;handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application-&gt;end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files-&gt;getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController-&gt;actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction-&gt;runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller-&gt;runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module-&gt;runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application-&gt;handleRequest()
#10 /var/www/html/prof-homepages/frontend/web/index.php(18): yii\base\Application-&gt;run()
#11 {main}</pre>