#parameter estimation of application that includes marginal 
#parameters and copula parameter and their distaces.
# and BIC criterion.

#remove.packages("rstan")
#if (file.exists(".RData")) file.remove(".RData")
#install.packages("devtools")
#install.packages("usethis")
library(usethis)
library(devtools)

Sys.setenv(PATH = paste("C:\\Rtools\\bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(PATH = paste("C:\\Rtools\\mingw_64\\bin", Sys.getenv("PATH"), sep=";"))
#install.packages("rstan")
library(rstan)



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(rhmc)
library(foreach)
library(doParallel)
library(parallel)
library(tictoc)
library(copula)
library(VineCopula)
library(DEoptim)

#====================== Estimation of marginal parameter =======================
#==================   First component                 ==========================

M=c(0.9,0.95,1,1.05,1.12,1.19,1.27,1.35,1.48,1.64,0,0,0,
    0.9,0.94,0.98,1.03,1.08,1.14,1.21,1.28,1.37,1.47,1.6,0,0,
    0.90,0.94,0.98,1.03,1.08,1.13,1.19,1.26,1.35,1.46,1.58,1.77,0,
    0.90,0.94,0.98,1.03,1.07,1.12,1.19,1.25,1.34,1.43,1.55,1.73,0,
    0.90,0.94,0.98,1.03,1.07,1.12,1.19,1.24,1.34,1.43,1.55,1.71,0,
    0.90,0.94,0.98,1.03,1.07,1.12,1.18,1.23,1.33,1.41,1.51,1.68,0,
    0.90,0.94,0.98,1.02,1.07,1.11,1.17,1.23,1.32,1.41,1.52,1.66,0,
    0.90,0.93,0.97,1.00,1.06,1.11,1.17,1.23,1.30,1.39,1.49,1.62,0,
    0.90,0.92,0.97,1.01,1.05,1.09,1.15,1.21,1.28,1.36,1.44,1.55,1.72,
    0.90,0.92,0.96,1.00,1.04,1.08,1.13,1.19,1.26,1.34,1.42,1.52,1.67,
    0.90,0.93,0.96,1.00,1.04,1.08,1.13,1.18,1.24,1.31,1.39,1.49,1.65,
    0.90,0.93,0.97,1.00,1.03,1.07,1.10,1.16,1.22,1.29,1.37,1.48,1.64,
    0.90,0.92,0.97,0.99,1.03,1.06,1.10,1.14,1.20,1.26,1.31,1.4,1.52,
    0.90,0.93,0.96,1.00,1.03,1.07,1.12,1.16,1.20,1.26,1.30,1.37,1.45,
    0.90,0.92,0.96,0.99,1.03,1.06,1.10,1.16,1.21,1.27,1.33,1.4,1.49,
    0.90,0.92,0.95,0.97,1.00,1.03,1.07,1.11,1.16,1.22,1.26,1.33,1.40,
    0.90,0.93,0.96,0.97,1.00,1.05,1.08,1.11,1.16,1.20,1.24,1.32,1.38,
    0.90,0.92,0.94,0.97,1.01,1.04,1.07,1.09,1.14,1.19,1.23,1.28,1.35,
    0.90,0.92,0.94,0.97,0.99,1.02,1.05,1.08,1.12,1.16,1.20,1.25,1.31,
    0.90,0.92,0.94,0.97,0.99,1.02,1.05,1.08,1.12,1.16,1.19,1.24,1.29)
N=matrix(M,20,13,byrow=T)




dx1=N[1,c(-1,-11,-12,-13)]-N[1,c(-10,-11,-12,-13)]
dx2=N[2,c(-1,-12,-13)]-N[2,c(-11,-12,-13)]
dx3=N[3,c(-1,-13)]-N[3,c(-12,-13)]
dx4=N[4,c(-1,-13)]-N[4,c(-12,-13)]
dx5=N[5,c(-1,-13)]-N[5,c(-12,-13)]
dx6=N[6,c(-1,-13)]-N[6,c(-12,-13)]
dx7=N[7,c(-1,-13)]-N[7,c(-12,-13)]
dx8=N[8,c(-1,-13)]-N[8,c(-12,-13)]
dx9=N[9,-1]-N[9,-13]
dx10=N[10,-1]-N[10,-13]

dy1=N[11,-1]-N[11,-13]
dy2=N[12,-1]-N[12,-13]
dy3=N[13,-1]-N[13,-13]
dy4=N[14,-1]-N[14,-13]
dy5=N[15,-1]-N[15,-13]
dy6=N[16,-1]-N[16,-13]
dy7=N[17,-1]-N[17,-13]
dy8=N[18,-1]-N[18,-13]
dy9=N[19,-1]-N[19,-13]
dy10=N[20,-1]-N[20,-13]


x=matrix(0,10,12)
x[1,]=c(dx1,0,0,0);x[2,]=c(dx2,0,0);x[3,]=c(dx3,0);x[4,]=c(dx4,0);x[5,]=c(dx5,0);x[6,]=c(dx6,0);x[7,]=c(dx7,0);x[8,]=c(dx8,0);x[9,]=dx9;x[10,]=dx10
y=matrix(0,10,12)
y[1,]=dy1;y[2,]=dy2;y[3,]=dy3;y[4,]=dy4;y[5,]=dy5;y[6,]=dy6;y[7,]=dy7;y[8,]=dy8;y[9,]=dy9;y[10,]=dy10




#========================gamma==============================================

stanmodelcodex<- "
data {
 int<lower=0> N1;          // number of increment of fist element
 int<lower=1> N2;          // number of increment of second element
 int<lower=1> N3;          // number of increment of third element
 int<lower=1> N4;          
 int<lower=1> N5;          
 int<lower=1> N6;          
 int<lower=1> N7;          
 int<lower=1> N8;          
 int<lower=1> N9;          
 int<lower=1> N10;          
 int<lower=1> M;          
  real  x1[N1];               
  real  x2[N2];               
  real  x3[N3];               
  real  x4[N4];               
  real  x5[N5];               
  real  x6[N6];               
  real  x7[N7];               // estimated treatment effects
  real  x8[N8];               // estimated treatment effects
  real  x9[N9];               // estimated treatment effects
  real  x10[N10];               // estimated treatment effects
 real t[M];
}
 parameters {
  real a;
   real  b;
   real c;
}
model {
     // likelihood
     for(n in 1:N1){
      real a1;
     a1=c*((t[n+1])^a-(t[n])^a);
x1[n]~gamma(a1,1/b);
 }
for(n in 1:N2){
  real a2;
  a2=c*((t[n+1])^a-(t[n])^a);
  x2[n]~gamma(a2,1/b);
}
for(n in 1:N3){
  real a3;
  a3=c*((t[n+1])^a-(t[n])^a);
  x3[n]~gamma(a3,1/b);
}
for(n in 1:N4){
  real a4;
  a4=c*((t[n+1])^a-(t[n])^a);
 x4[n]~gamma(a4,1/b);
}
for(n in 1:N5){
  real a5;
  a5=c*((t[n+1])^a-(t[n])^a);
  x5[n]~gamma(a5,1/b);
}
for(n in 1:N6){
  real a6;
  a6=c*((t[n+1])^a-(t[n])^a);
  x6[n]~gamma(a6,1/b);
}
for(n in 1:N7){
  real a7;
  a7=c*((t[n+1])^a-(t[n])^a);
  x7[n]~gamma(a7,1/b);
}
for(n in 1:N8){
  real a8;
  a8=c*((t[n+1])^a-(t[n])^a);
  x8[n]~gamma(a8,1/b);
}
for(n in 1:N9){
  real a9;
  a9=c*((t[n+1])^a-(t[n])^a);
  x9[n]~gamma(a9,1/b);
}
for(n in 1:N10){
  real a10;
  a10=c*((t[n+1])^a-(t[n])^a);
  x10[n]~gamma(a10,1/b);
}
}
"


crackdatax<- list(
  N1 =8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=12,
  x1=c(x[1,-c(9:12)]),
  x2=c(x[2,-c(10:12)]),
  x3=c(x[3,-c(11:12)]),
  x4=c(x[4,-c(11:12)]),
  x5=c(x[5,-c(11:12)]),
  x6=c(x[6,-c(11:12)]),
  x7=c(x[7,-c(11:12)]),
  x8=c(x[8,-c(11:12)]),
  x9=c(x[9,-c(12)]),
  x10=c(x[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11)
)

initf1 <- function() {
  list(a=1.5,b=0.0034,c=5685)
}


library(rstan)

fit = stan(model_code=stanmodelcodex, data=crackdatax,iter=1000,cores=3,chains=1)
print(fit, probs=c(0.025, 0.975),digits=4)
mx <- as.matrix(fit)


#==============Wiener==========================================================


stanmodelcodex<- "
data {
  int<lower=1> N1;          // number of increment of fist element
  int<lower=1> N2;          // number of increment of second element
  int<lower=1> N3;          // number of increment of third element
  int<lower=1> N4;          // number of increment of forth  element
  int<lower=1> N5;          // number of increment of 5th element
  int<lower=1> N6;          // number of increment of 5th element
  int<lower=1> N7;          // number of increment of 5th element
  int<lower=1> N8;          // number of increment of 5th element
  int<lower=1> N9;          // number of increment of 5th element
  int<lower=1> N10;          // number of increment of 5th element
  int<lower=1> M;          // number of inspections
  real  x1[N1];               // estimated treatment effects
  real  x2[N2];               // estimated treatment effects
  real  x3[N3];               // estimated treatment effects
  real  x4[N4];               // estimated treatment effects
  real  x5[N5];               // estimated treatment effects
  real  x6[N6];               // estimated treatment effects
  real  x7[N7];               // estimated treatment effects
  real  x8[N8];               // estimated treatment effects
  real  x9[N9];               // estimated treatment effects
  real  x10[N10];               // estimated treatment effects
  real t[M];
}
parameters {
  real a;
   real  b;
   real c;
}
model {
// prior distributions:
    a ~ uniform(0, 10);        
    b ~ uniform(0, 10);
    c ~ uniform(0, 200);
     // likelihood
     for(n in 1:N1){
      real a11;
      real a12;
     a11=c*((t[n+1])^a-(t[n])^a);
     a12=b*sqrt((t[n+1])^a-(t[n])^a);
x1[n]~normal(a11,a12);
 }
for(n in 1:N2){
  real a21;
  real a22;
  a21=c*((t[n+1])^a-(t[n])^a);
  a22=b*sqrt((t[n+1])^a-(t[n])^a);
x2[n]~normal(a21,a22);
}
for(n in 1:N3){
  real a31;
  real a32;
  a31=c*((t[n+1])^a-(t[n])^a);
  a32=b*sqrt((t[n+1])^a-(t[n])^a);
x3[n]~normal(a31,a32);
}
for(n in 1:N4){
  real a41;
  real a42;
  a41=c*((t[n+1])^a-(t[n])^a);
  a42=b*sqrt((t[n+1])^a-(t[n])^a);
x4[n]~normal(a41,a42);
}
for(n in 1:N5){
  real a51;
  real a52;
  a51=c*((t[n+1])^a-(t[n])^a);
  a52=b*sqrt((t[n+1])^a-(t[n])^a);
x5[n]~normal(a51,a52);
}
for(n in 1:N6){
  real a61;
  real a62;
  a61=c*((t[n+1])^a-(t[n])^a);
  a62=b*sqrt((t[n+1])^a-(t[n])^a);
x6[n]~normal(a61,a62);
}
for(n in 1:N7){
 real a71;
  real a72;
  a71=c*((t[n+1])^a-(t[n])^a);
  a72=b*sqrt((t[n+1])^a-(t[n])^a);
x7[n]~normal(a71,a72);
}
for(n in 1:N8){
 real a81;
  real a82;
  a81=c*((t[n+1])^a-(t[n])^a);
  a82=b*sqrt((t[n+1])^a-(t[n])^a);
x8[n]~normal(a81,a82);
}
for(n in 1:N9){
 real a91;
  real a92;
  a91=c*((t[n+1])^a-(t[n])^a);
  a92=b*sqrt((t[n+1])^a-(t[n])^a);
x9[n]~normal(a91,a92);
}
for(n in 1:N10){
 real a101;
  real a102;
  a101=c*((t[n+1])^a-(t[n])^a);
  a102=b*sqrt((t[n+1])^a-(t[n])^a);
x10[n]~normal(a101,a102);
}
}
"






crackdatax<- list(
  N1 =8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=13,
  x1=c(x[1,-c(9:12)]),
  x2=c(x[2,-c(10:12)]),
  x3=c(x[3,-c(11:12)]),
  x4=c(x[4,-c(11:12)]),
  x5=c(x[5,-c(11:12)]),
  x6=c(x[6,-c(11:12)]),
  x7=c(x[7,-c(11:12)]),
  x8=c(x[8,-c(11:12)]),
  x9=c(x[9,-c(12)]),
  x10=c(x[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12)
)


library(rstan)
fit = stan(model_code=stanmodelcodex, data=crackdatax, iter=2000,warmup=1000,cores=4,chains=1)
fit
print(fit, probs=c(0.025, 0.975),digits=5)


traceplot(fit, include = TRUE,inc_warmup = FALSE, window = NULL,
          nrow = NULL, ncol = NULL)



#================================IG process====================================

stanmodelcodex<- "
data {
  int<lower=1> N1;          // number of increment of fist element
 int<lower=1> N2;          // number of increment of second element
  int<lower=1> N3;          // number of increment of second element
    int<lower=1> N4;          // number of increment of second element
    int<lower=1> N5;          // number of increment of second element
    int<lower=1> N6;          // number of increment of second element
    int<lower=1> N7;          // number of increment of second element
    int<lower=1> N8;          // number of increment of second element
    int<lower=1> N9;          // number of increment of second element
     int<lower=1> N10;          // number of increment of second element
  int<lower=1> M;          // number of inspections
  real  x1[N1];               // estimated treatment effects
real  x2[N2];               // estimated treatment effects
real  x3[N3];               // estimated treatment effects
real  x4[N4];               // estimated treatment effects
real  x5[N5];               // estimated treatment effects
real  x6[N6];               // estimated treatment effects
real  x7[N7];               // estimated treatment effects
real  x8[N8];               // estimated treatment effects
real  x9[N9];               // estimated treatment effects
real  x10[N10];               // estimated treatment effects
real t[M];
}
parameters {
  real a;
   real  b;
   real c;
}
model {
// prior distributions:
    a ~ uniform(0, 100);        
    b ~ uniform(0, 300000);
    c ~ uniform(0, 100);
     // likelihood
     for(n in 1:N1){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x1[n])^3))-(b*(x1[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x1[n]);
     }
   for(n in 1:N2){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x2[n])^3))-(b*(x2[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x2[n]);
   }
   for(n in 1:N3){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x3[n])^3))-(b*(x3[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x3[n]);
   }
   for(n in 1:N4){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x4[n])^3))-(b*(x4[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x4[n]);
   }
    for(n in 1:N5){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x5[n])^3))-(b*(x5[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x5[n]);
    }
    for(n in 1:N6){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x6[n])^3))-(b*(x6[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x6[n]);
    }
    for(n in 1:N7){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x7[n])^3))-(b*(x7[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x7[n]);
    }
    for(n in 1:N8){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x8[n])^3))-(b*(x8[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x8[n]);
    }
    for(n in 1:N9){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x9[n])^3))-(b*(x9[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x9[n]);
    }
    for(n in 1:N10){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(x10[n])^3))-(b*(x10[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*x10[n]);
   }
}
"


crackdatax<- list(
  N1 =8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=13,
  x1=c(x[1,-c(9:12)]),
  x2=c(x[2,-c(10:12)]),
  x3=c(x[3,-c(11:12)]),
  x4=c(x[4,-c(11:12)]),
  x5=c(x[5,-c(11:12)]),
  x6=c(x[6,-c(11:12)]),
  x7=c(x[7,-c(11:12)]),
  x8=c(x[8,-c(11:12)]),
  x9=c(x[9,-c(12)]),
  x10=c(x[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12)
)
library(rstan)
fit = stan(model_code=stanmodelcodex, data=crackdatax, iter=4000,warmup=2000,cores=4)
print(fit, probs=c(0.025, 0.975),digits=8)



#=========================  Second component========================
#========gamma================================

stanmodelcodey<- "
data {
  int<lower=1> N1;          // number of increment of fist element
int<lower=1> N2;          // number of increment of second element
int<lower=1> N3;          // number of increment of third element
int<lower=1> N4;          // number of increment of forth  element
int<lower=1> N5;          // number of increment of 5th element
int<lower=1> N6;          // number of increment of 5th element
int<lower=1> N7;          // number of increment of 5th element
int<lower=1> N8;          // number of increment of 5th element
int<lower=1> N9;          // number of increment of 5th element
int<lower=1> N10;          // number of increment of 5th element
  int<lower=1> M;          // number of inspections
  real  y1[N1];               // estimated treatment effects
real  y2[N2];               // estimated treatment effects
real  y3[N3];               // estimated treatment effects
real  y4[N4];               // estimated treatment effects
real  y5[N5];               // estimated treatment effects
real  y6[N6];               // estimated treatment effects
real  y7[N7];               // estimated treatment effects
real  y8[N8];               // estimated treatment effects
real  y9[N9];               // estimated treatment effects
real  y10[N10];               // estimated treatment effects
  real t[M];
}
parameters {
   real a;
   real  b;
   real<lower=0>  c;
}
model {
// prior distributions:
    a ~ uniform(0, 5);        
    b ~ uniform(0, 2);
    c ~ uniform(0, 8000);
     // likelihood
     for(n in 1:N1){
      real a1;
     a1=c*(t[n+1])^a-c*(t[n])^a;
y1[n]~gamma(a1,1/b);
 }
for(n in 1:N2){
  real a2;
  a2=c*(t[n+1])^a-c*(t[n])^a;
  y2[n]~gamma(a2,1/b);
}
for(n in 1:N3){
  real a3;
  a3=c*(t[n+1])^a-c*(t[n])^a;
  y3[n]~gamma(a3,1/b);
}
for(n in 1:N4){
  real a4;
  a4=c*(t[n+1])^a-c*(t[n])^a;
  y4[n]~gamma(a4,1/b);
}
for(n in 1:N5){
  real a5;
  a5=c*(t[n+1])^a-c*(t[n])^a;
  y5[n]~gamma(a5,1/b);
}
for(n in 1:N6){
  real a6;
  a6=c*(t[n+1])^a-c*(t[n])^a;
  y6[n]~gamma(a6,1/b);
}
for(n in 1:N7){
  real a7;
  a7=c*(t[n+1])^a-c*(t[n])^a;
  y7[n]~gamma(a7,1/b);
}
for(n in 1:N8){
  real a8;
  a8=c*(t[n+1])^a-c*(t[n])^a;
  y8[n]~gamma(a8,1/b);
}
for(n in 1:N9){
  real a9;
  a9=c*(t[n+1])^a-c*(t[n])^a;
  y9[n]~gamma(a9,1/b);
}
for(n in 1:N10){
  real a10;
  a10=c*(t[n+1])^a-c*(t[n])^a;
  y10[n]~gamma(a10,1/b);
}
}
"


initf1 <- function() {
  list(a=1,b=0.01,c=4000)
}



crackdatay<- list(
  N1 =8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=13,
  y1=c(y[1,-c(9:12)]),
  y2=c(y[2,-c(10:12)]),
  y3=c(y[3,-c(11:12)]),
  y4=c(y[4,-c(11:12)]),
  y5=c(y[5,-c(11:12)]),
  y6=c(y[6,-c(11:12)]),
  y7=c(y[7,-c(11:12)]),
  y8=c(y[8,-c(11:12)]),
  y9=c(y[9,-c(12)]),
  y10=c(y[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12)
)



library(rstan)
fit = stan(model_code=stanmodelcodey,init=initf1, data=crackdatay, iter=1000,cores=3)
fit
print(fit, probs=c(0.025, 0.975),digits=5)
my <- as.matrix(fit)

#================================wiener========================================

stanmodelcodey<- "
data {
  int<lower=1> N1;          // number of increment of fist element
int<lower=1> N2;          // number of increment of second element
int<lower=1> N3;          // number of increment of third element
int<lower=1> N4;          // number of increment of forth  element
int<lower=1> N5;          // number of increment of 5th element
int<lower=1> N6;          // number of increment of 5th element
int<lower=1> N7;          // number of increment of 5th element
int<lower=1> N8;          // number of increment of 5th element
int<lower=1> N9;          // number of increment of 5th element
int<lower=1> N10;          // number of increment of 5th element
  int<lower=1> M;          // number of inspections
  real  y1[N1];               // estimated treatment effects
real  y2[N2];               // estimated treatment effects
real  y3[N3];               // estimated treatment effects
real  y4[N4];               // estimated treatment effects
real  y5[N5];               // estimated treatment effects
real  y6[N6];               // estimated treatment effects
real  y7[N7];               // estimated treatment effects
real  y8[N8];               // estimated treatment effects
real  y9[N9];               // estimated treatment effects
real  y10[N10];               // estimated treatment effects
  real t[M];
}
parameters {
   real a;
   real  b;
   real<lower=0>  c;
}
model {
// prior distributions:
    a ~ uniform(0, 10);        
    b ~ uniform(0, 10);
    c ~ uniform(0, 200);
     // likelihood
     for(n in 1:N1){
      real a11;
      real a12;
     a11=c*((t[n+1])^a-(t[n])^a);
     a12=b*sqrt((t[n+1])^a-(t[n])^a);
y1[n]~normal(a11,a12);
 }
for(n in 1:N2){
  real a21;
  real a22;
  a21=c*((t[n+1])^a-(t[n])^a);
  a22=b*sqrt((t[n+1])^a-(t[n])^a);
y2[n]~normal(a21,a22);
}
for(n in 1:N3){
  real a31;
  real a32;
  a31=c*((t[n+1])^a-(t[n])^a);
  a32=b*sqrt((t[n+1])^a-(t[n])^a);
y3[n]~normal(a31,a32);
}
for(n in 1:N4){
  real a41;
  real a42;
  a41=c*((t[n+1])^a-(t[n])^a);
  a42=b*sqrt((t[n+1])^a-(t[n])^a);
y4[n]~normal(a41,a42);
}
for(n in 1:N5){
  real a51;
  real a52;
  a51=c*((t[n+1])^a-(t[n])^a);
  a52=b*sqrt((t[n+1])^a-(t[n])^a);
y5[n]~normal(a51,a52);
}
for(n in 1:N6){
  real a61;
  real a62;
  a61=c*((t[n+1])^a-(t[n])^a);
  a62=b*sqrt((t[n+1])^a-(t[n])^a);
y6[n]~normal(a61,a62);
}
for(n in 1:N7){
 real a71;
  real a72;
  a71=c*((t[n+1])^a-(t[n])^a);
  a72=b*sqrt((t[n+1])^a-(t[n])^a);
y7[n]~normal(a71,a72);
}
for(n in 1:N8){
 real a81;
  real a82;
  a81=c*((t[n+1])^a-(t[n])^a);
  a82=b*sqrt((t[n+1])^a-(t[n])^a);
y8[n]~normal(a81,a82);
}
for(n in 1:N9){
 real a91;
  real a92;
  a91=c*((t[n+1])^a-(t[n])^a);
  a92=b*sqrt((t[n+1])^a-(t[n])^a);
y9[n]~normal(a91,a92);
}
for(n in 1:N10){
 real a101;
  real a102;
  a101=c*((t[n+1])^a-(t[n])^a);
  a102=b*sqrt((t[n+1])^a-(t[n])^a);
y10[n]~normal(a101,a102);
}
}
"


initf1 <- function() {
  list(a=1,b=0.01,c=4000)
}



crackdatay<- list(
  N1 =8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=13,
  y1=c(y[1,-c(9:12)]),
  y2=c(y[2,-c(10:12)]),
  y3=c(y[3,-c(11:12)]),
  y4=c(y[4,-c(11:12)]),
  y5=c(y[5,-c(11:12)]),
  y6=c(y[6,-c(11:12)]),
  y7=c(y[7,-c(11:12)]),
  y8=c(y[8,-c(11:12)]),
  y9=c(y[9,-c(12)]),
  y10=c(y[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12)
)



library(rstan)
fit = stan(model_code=stanmodelcodey, data=crackdatay, iter=4000,warmup=2000,cores=4)
fit
print(fit, probs=c(0.025, 0.975),digits=5)


#================================IG process====================================

stanmodelcodey<- "
data {
  int<lower=1> N1;          // number of increment of fist element
 int<lower=1> N2;          // number of increment of second element
  int<lower=1> N3;          // number of increment of second element
    int<lower=1> N4;          // number of increment of second element
    int<lower=1> N5;          // number of increment of second element
    int<lower=1> N6;          // number of increment of second element
    int<lower=1> N7;          // number of increment of second element
    int<lower=1> N8;          // number of increment of second element
    int<lower=1> N9;          // number of increment of second element
     int<lower=1> N10;          // number of increment of second element
  int<lower=1> M;          // number of inspections
  real  y1[N1];               // estimated treatment effects
real  y2[N2];               // estimated treatment effects
real  y3[N3];               // estimated treatment effects
real  y4[N4];               // estimated treatment effects
real  y5[N5];               // estimated treatment effects
real  y6[N6];               // estimated treatment effects
real  y7[N7];               // estimated treatment effects
real  y8[N8];               // estimated treatment effects
real  y9[N9];               // estimated treatment effects
real  y10[N10];               // estimated treatment effects
real t[M];
}
parameters {
  real a;
   real  b;
   real c;
}
model {
// prior distributions:
    a ~ uniform(0, 100);        
    b ~ uniform(0, 300000);
    c ~ uniform(0, 100);
     // likelihood
     for(n in 1:N1){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y1[n])^3))-(b*(y1[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y1[n]);
     }
   for(n in 1:N2){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y2[n])^3))-(b*(y2[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y2[n]);
   }
   for(n in 1:N3){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y3[n])^3))-(b*(y3[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y3[n]);
   }
   for(n in 1:N4){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y4[n])^3))-(b*(y4[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y4[n]);
   }
    for(n in 1:N5){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y5[n])^3))-(b*(y5[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y5[n]);
    }
    for(n in 1:N6){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y6[n])^3))-(b*(y6[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y6[n]);
    }
    for(n in 1:N7){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y7[n])^3))-(b*(y7[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y7[n]);
    }
    for(n in 1:N8){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y8[n])^3))-(b*(y8[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y8[n]);
    }
    for(n in 1:N9){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y9[n])^3))-(b*(y9[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y9[n]);
    }
    for(n in 1:N10){
    target+=0.5*(log(b)+2*log(((t[n+1])^a-(t[n])^a))-log(2*3.141593*(y10[n])^3))-(b*(y10[n]-c*((t[n+1])^a-(t[n])^a))^2)/(2*c^2*y10[n]);
   }
}
"


crackdatay<- list(
  N1 =8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=13,
  y1=c(y[1,-c(9:12)]),
  y2=c(y[2,-c(10:12)]),
  y3=c(y[3,-c(11:12)]),
  y4=c(y[4,-c(11:12)]),
  y5=c(y[5,-c(11:12)]),
  y6=c(y[6,-c(11:12)]),
  y7=c(y[7,-c(11:12)]),
  y8=c(y[8,-c(11:12)]),
  y9=c(y[9,-c(12)]),
  y10=c(y[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12)
)
library(rstan)
fit = stan(model_code=stanmodelcodey, data=crackdatay, iter=4000,warmup=2000,cores=4)
print(fit, probs=c(0.025, 0.975),digits=8)


#============================== Estimation of copula  parameter================

M=c(0.9,0.95,1,1.05,1.12,1.19,1.27,1.35,1.48,1.64,0,0,0,
    0.9,0.94,0.98,1.03,1.08,1.14,1.21,1.28,1.37,1.47,1.6,0,0,
    0.90,0.94,0.98,1.03,1.08,1.13,1.19,1.26,1.35,1.46,1.58,1.77,0,
    0.90,0.94,0.98,1.03,1.07,1.12,1.19,1.25,1.34,1.43,1.55,1.73,0,
    0.90,0.94,0.98,1.03,1.07,1.12,1.19,1.24,1.34,1.43,1.55,1.71,0,
    0.90,0.94,0.98,1.03,1.07,1.12,1.18,1.23,1.33,1.41,1.51,1.68,0,
    0.90,0.94,0.98,1.02,1.07,1.11,1.17,1.23,1.32,1.41,1.52,1.66,0,
    0.90,0.93,0.97,1.00,1.06,1.11,1.17,1.23,1.30,1.39,1.49,1.62,0,
    0.90,0.92,0.97,1.01,1.05,1.09,1.15,1.21,1.28,1.36,1.44,1.55,1.72,
    0.90,0.92,0.96,1.00,1.04,1.08,1.13,1.19,1.26,1.34,1.42,1.52,1.67,
    0.90,0.93,0.96,1.00,1.04,1.08,1.13,1.18,1.24,1.31,1.39,1.49,1.65,
    0.90,0.93,0.97,1.00,1.03,1.07,1.10,1.16,1.22,1.29,1.37,1.48,1.64,
    0.90,0.92,0.97,0.99,1.03,1.06,1.10,1.14,1.20,1.26,1.31,1.4,1.52,
    0.90,0.93,0.96,1.00,1.03,1.07,1.12,1.16,1.20,1.26,1.30,1.37,1.45,
    0.90,0.92,0.96,0.99,1.03,1.06,1.10,1.16,1.21,1.27,1.33,1.4,1.49,
    0.90,0.92,0.95,0.97,1.00,1.03,1.07,1.11,1.16,1.22,1.26,1.33,1.40,
    0.90,0.93,0.96,0.97,1.00,1.05,1.08,1.11,1.16,1.20,1.24,1.32,1.38,
    0.90,0.92,0.94,0.97,1.01,1.04,1.07,1.09,1.14,1.19,1.23,1.28,1.35,
    0.90,0.92,0.94,0.97,0.99,1.02,1.05,1.08,1.12,1.16,1.20,1.25,1.31,
    0.90,0.92,0.94,0.97,0.99,1.02,1.05,1.08,1.12,1.16,1.19,1.24,1.29)
N=matrix(M,20,13,byrow=T)




dx1=N[1,c(-1,-11,-12,-13)]-N[1,c(-10,-11,-12,-13)]
dx2=N[2,c(-1,-12,-13)]-N[2,c(-11,-12,-13)]
dx3=N[3,c(-1,-13)]-N[3,c(-12,-13)]
dx4=N[4,c(-1,-13)]-N[4,c(-12,-13)]
dx5=N[5,c(-1,-13)]-N[5,c(-12,-13)]
dx6=N[6,c(-1,-13)]-N[6,c(-12,-13)]
dx7=N[7,c(-1,-13)]-N[7,c(-12,-13)]
dx8=N[8,c(-1,-13)]-N[8,c(-12,-13)]
dx9=N[9,-1]-N[9,-13]
dx10=N[10,-1]-N[10,-13]

dy1=N[11,-1]-N[11,-13]
dy2=N[12,-1]-N[12,-13]
dy3=N[13,-1]-N[13,-13]
dy4=N[14,-1]-N[14,-13]
dy5=N[15,-1]-N[15,-13]
dy6=N[16,-1]-N[16,-13]
dy7=N[17,-1]-N[17,-13]
dy8=N[18,-1]-N[18,-13]
dy9=N[19,-1]-N[19,-13]
dy10=N[20,-1]-N[20,-13]


x=matrix(0,10,12)
x[1,]=c(dx1,0,0,0);x[2,]=c(dx2,0,0);x[3,]=c(dx3,0);x[4,]=c(dx4,0);x[5,]=c(dx5,0);x[6,]=c(dx6,0);x[7,]=c(dx7,0);x[8,]=c(dx8,0);x[9,]=dx9;x[10,]=dx10
y=matrix(0,10,12)
y[1,]=dy1;y[2,]=dy2;y[3,]=dy3;y[4,]=dy4;y[5,]=dy5;y[6,]=dy6;y[7,]=dy7;y[8,]=dy8;y[9,]=dy9;y[10,]=dy10




#========================  c t^ a Frank==================================================

modelcodecop="
data {
  int<lower=1> N1;          // number of increment of fist element
  int<lower=1> N2;          // number of increment of second element
  int<lower=1> N3;          // number of increment of third element
  int<lower=1> N4;          // number of increment of forth  element
  int<lower=1> N5;          // number of increment of 5th element
  int<lower=1> N6;          // number of increment of 5th element
  int<lower=1> N7;          // number of increment of 5th element
  int<lower=1> N8;          // number of increment of 5th element
  int<lower=1> N9;          // number of increment of 5th element
  int<lower=1> N10;          // number of increment of 5th element
  int<lower=1> M;          // number of inspections
  real  y1[N1];               // estimated treatment effects
  real  y2[N2];               // estimated treatment effects
  real  y3[N3];               // estimated treatment effects
  real  y4[N4];               // estimated treatment effects
  real  y5[N5];               // estimated treatment effects
  real  y6[N6];               // estimated treatment effects
  real  y7[N7];               // estimated treatment effects
  real  y8[N8];               // estimated treatment effects
  real  y9[N9];               // estimated treatment effects
  real  y10[N10];               // estimated treatment effects
  real  x1[N1];               // estimated treatment effects
  real  x2[N2];               // estimated treatment effects
  real  x3[N3];               // estimated treatment effects
  real  x4[N4];               // estimated treatment effects
  real  x5[N5];               // estimated treatment effects
  real  x6[N6];               // estimated treatment effects
  real  x7[N7];               // estimated treatment effects
  real  x8[N8];               // estimated treatment effects
  real  x9[N9];               // estimated treatment effects
  real  x10[N10];  
  real t[M];
  }
parameters{
 real  alpha;
  }
model{
 //priors
 alpha ~ uniform(0.01,14);
     //likelihood
      if(alpha!=0){
   for(n in 1:N1){
    real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
      p[1]= gamma_cdf(x1[n],a11,1/0.0034);
      p[2]= gamma_cdf(y1[n],a12,1/0.0027);
      target+=log(alpha)
      +alpha*(p[1]+p[2])
      +log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N2){
     real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
      p[1]= gamma_cdf(x2[n],a11,1/0.0034);
      p[2]= gamma_cdf(y2[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N3){
    real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.9*(t[n])^1.25;
     real  p[2];
      p[1]= gamma_cdf(x3[n],a11,1/0.0034);
      p[2]= gamma_cdf(y3[n],a12,1/0.0027);
     target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N4){
     real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
       p[1]= gamma_cdf(x4[n],a11,1/0.0034);
      p[2]= gamma_cdf(y4[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N5){
    real a11=4675.2*(t[n+1])^1.39-4675.2*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
      p[1]= gamma_cdf(x5[n],a11,1/0.0034);
      p[2]= gamma_cdf(y5[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N6){
     real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
       p[1]= gamma_cdf(x6[n],a11,1/0.0034);
      p[2]= gamma_cdf(y6[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N7){
    real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
       p[1]= gamma_cdf(x7[n],a11,1/0.0034);
      p[2]= gamma_cdf(y7[n],a12,1/0.0027);
   target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N8){
    real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
      p[1]= gamma_cdf(x8[n],a11,1/0.0034);
      p[2]= gamma_cdf(y8[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N9){
     real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
     p[1]= gamma_cdf(x9[n],a11,1/0.0034);
      p[2]= gamma_cdf(y9[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    for(n in 1:N10){
    real a11=4675.25*(t[n+1])^1.39-4675.25*(t[n])^1.39;
    real a12=2501.98*(t[n+1])^1.25-2501.98*(t[n])^1.25;
     real  p[2];
      p[1]= gamma_cdf(x10[n],a11,1/0.0034);
      p[2]= gamma_cdf(y10[n],a12,1/0.0027);
    target+=log(alpha)+alpha*(p[1]+p[2])+log(exp(alpha)-1)
      -2*log(exp(alpha)-exp(alpha*(p[1]))-exp(alpha*(p[2]))+exp(alpha*(p[1]+p[2])));
         }
    }
}
"


dat<- list(
  N1=8,
  N2=9,
  N3=10,
  N4=10,
  N5=10,
  N6=10,
  N7=10,
  N8=10,
  N9=11,
  N10=11,
  M=13,
  x1=c(x[1,-c(9:12)]),
  x2=c(x[2,-c(10:12)]),
  x3=c(x[3,-c(11:12)]),
  x4=c(x[4,-c(11:12)]),
  x5=c(x[5,-c(11:12)]),
  x6=c(x[6,-c(11:12)]),
  x7=c(x[7,-c(11:12)]),
  x8=c(x[8,-c(11:12)]),
  x9=c(x[9,-c(12)]),
  x10=c(x[10,-c(12)]),
  y1=c(y[1,-c(9:12)]),
  y2=c(y[2,-c(10:12)]),
  y3=c(y[3,-c(11:12)]),
  y4=c(y[4,-c(11:12)]),
  y5=c(y[5,-c(11:12)]),
  y6=c(y[6,-c(11:12)]),
  y7=c(y[7,-c(11:12)]),
  y8=c(y[8,-c(11:12)]),
  y9=c(y[9,-c(12)]),
  y10=c(y[10,-c(12)]),
  t=c(0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12)
)

initf1 <- function() {
  list(alpha=5)
}

fit=stan(model_code=modelcodecop,data=dat,init=initf1,iter=4000,thin=10,chains=3,cores=4)
alpha_draws= extract(fit)$alpha
print(fit, probs=c(0.025, 0.975),digits=5)



#===================================BIC=========================================

n=10
pdx=matrix(0,n,11)
pdy=matrix(0,n,11)

t=c(seq(0,0.11,0.01))
for(h in 1:n){
  for(i in 1:11){
    pdx[h,i]=pgamma(x[h,i],shape=4675.25*(t[i+1])^1.4-4675.25*(t[i])^1.4,scale=0.0034)
    pdy[h,i]=pgamma(y[h,i],shape=(2501.98*(t[i+1])^1.25-2501.98*(t[i])^1.25),scale=0.0027)
  }
}



Lg=function(p){
  sum(sapply(1:8,function(i)   (-log (dCopula(c(pdx[1,i],pdy[1,i]), gumbelCopula(p)))
  )    ) )
  +sum(sapply(1:9,function(i)   (-log (dCopula(c(pdx[2,i],pdy[2,i]), gumbelCopula(p)))
  )    ) )
  + sum(sapply(3:8, function(h)      sum(sapply(1:10,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), gumbelCopula(p))))     ) )    ))
  + sum(sapply(9:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), gumbelCopula(p))))     ) )    ))
  }
2*Lg(1.387)+log(99)

Lf=function(p){
  sum(sapply(1:8,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), frankCopula(p)))
  )    ) )
  +sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), frankCopula(p)))
  )    ) )
  + sum(sapply(3:8, function(h)      sum(sapply(1:10,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), frankCopula(p))))     ) )    ))
  + sum(sapply(9:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), frankCopula(p))))     ) )    ))
  }
2*Lf(0.1307)+log(99)


Lj=function(p){
  sum(sapply(1:8,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), joeCopula(p)))
  )    ) )
  +sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), joeCopula(p)))
  )    ) )
  + sum(sapply(3:8, function(h)      sum(sapply(1:10,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), joeCopula(p))))     ) )    ))
  + sum(sapply(9:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), joeCopula(p))))     ) )    ))
  }

2*Lj(1.535)+log(99)

Lc=function(p){
  sum(sapply(1:8,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), claytonCopula(p)))
  )    ) )
  +sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), claytonCopula(p)))
  )    ) )
  + sum(sapply(3:8, function(h)      sum(sapply(1:10,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), claytonCopula(p))))     ) )    ))
  + sum(sapply(9:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), claytonCopula(p))))     ) )    ))
  }

2*Lc(0.492)+log(99)

Lfgm=function(p){
  sum(sapply(1:8,function(i)   (  -log (dCopula(c(pdx[1,i],pdy[1,i]), fgmCopula(p)))
  )    ) )
  +sum(sapply(1:9,function(i)   (  -log (dCopula(c(pdx[2,i],pdy[2,i]), fgmCopula(p)))
  )    ) )
  + sum(sapply(3:8, function(h)      sum(sapply(1:10,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), fgmCopula(p))))     ) )    ))
  + sum(sapply(9:10, function(h)      sum(sapply(1:11,function(i)  (  -log (dCopula(c(pdx[h,i],pdy[h,i]), fgmCopula(p))))     ) )    ))
}


<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>