library(HDInterval) 
library(stats) 
library(MASS) 
library(qualityTools) 
library(magick) 

CI.com.mean.wei<-function(M,m,k,T,n1,n2,k1,k2,mu1,mu2)
{								 #Start 
  
  # x = data 
  # y = data 
  # M = number of simulation runs 
  # m = number of generalized computations
  # k = number of sample case
  # T = number of Bayesian computation  
  # n1 = number of sample size of x 
  # n2 = number of sample size of y 
  # a1 = scale parameter of x 
  # a2 = scale parameter of y 
  # k1 = shape parameter of x   
  # k2 = shape parameter of y 
  
  CP.phiGCI      		= rep(0,M)
  Length.phiGCI  		= rep(0,M)
  CP.phiMV       		= rep(0,M)
  Length.phiMV   	= rep(0,M)
  CP.phiBAYE     	= rep(0,M)
  Length.phiBAYE  	= rep(0,M)
  CP.phiBAYE2    	= rep(0,M)
  Length.phiBAYE2	= rep(0,M)
  
  Rphi		= rep(0,m)
  areal      	= rep(0,k)
  kreal       	= rep(0,k)
  amlex 		= rep(0,k)
  kmlexx 	= rep(0,k)
  a.starx		= rep(0,k)
  k.starxx	= rep(0,k)
  V.mhx       	= rep(0,k)
  V.mh        	= rep(0,k)
  V.mh.r      	= rep(0,k)
  MU_real     	= rep(0,k)
  MU.hat.m    	= rep(0,k)
  L.CIx       	= rep(0,k)
  U.CIx      	= rep(0,k)
  xx1        	= matrix(0,k,n1)
  a.calRvar   	= matrix(0,k,1)
  k.calRvar   	= matrix(0,k,1)
  zvalue	    	= qnorm(0.975)
  zvalue2     	= qnorm(0.025)
  alpha	          	= 0.05
  err	        	= 10^(-6)
  a1          	= mu1/gamma(1+(1/k1))
  a2          	= mu2/gamma(1+(1/k2))
  q           	= k/2
  v1          	= 0.1
  v2          	= 0.1
  z1          	= 0.1
  z2          	= 0.1
  kmlex.baye  	= c()
  bmlex.baye  	= c()
  amlex.baye  	= c()
  MU.hat.baye 	= c()
  phi.hat.baye 	= c()
  V.mhx.b 	= c()
  MU.gibbs   	= matrix(0,k,T)
  V.mu.gibbs  	= matrix(0,k,T)
  
  for (i in 1:M) { # begin iteration 
    gen.x <- function(k,n1,n2,a1,a2,k1,k2) 					#begin function for data
    {
      for (a in 1:k) 
      { 
        if(a <= q){ 
          x 		= rweibull(n1,k1,a1)
          zx 	= log(x)
          zbarx 	= mean(zx)
          s2x  	= sum((zx-zbarx)^2)
          sx		= sqrt(s2x)
          ku1 	= (pi/sqrt(6))*(sqrt(n1-1))/sx
          k01	= ku1
          S11	= mean(zx)
          w1	= x^k01	
          S21	= sum(w1)
          S31 	= sum(w1*zx)
          S41	= sum(w1*(zx^2))
          F1		= (1/k01)+S11-(S31/S21)
          f1		= (1/(k01^2)) + ((S21*S41)-(S31^2))/(S21^2)
          kmlex	= k01 + (F1/f1)
          while(abs(F1)>=err)	{
            k01	= kmlex
            S11	= mean(zx)
            w1	= x^k01	
            S21	= sum(w1)
            S31 	= sum(w1*zx)
            S41	= sum(w1*(zx^2))
            F1	= (1/k01)+S11-(S31/S21)
            f1	= (1/(k01^2)) + ((S21*S41)-(S31^2))/(S21^2)
            kmlex	= k01 + (F1/f1)
          }
          areal[a]	= a1
          kreal[a]	= k1
          amlex[a]	= ((1/n1)*S21)^(1/kmlex)
          kmlexx[a] 	= kmlex
          x1 	= rbind(x)
          xx1[a,]	= x1
          
        }else
        { 
          x 	<- rweibull(n2,k2,a2)
          zx 	= log(x)
          zbarx 	= mean(zx)
          s2x  	= sum((zx-zbarx)^2)
          sx		= sqrt(s2x)
          ku1 	= (pi/sqrt(6))*(sqrt(n2-1))/sx
          k01	= ku1
          S11	= mean(zx)
          w1	= x^k01	
          S21	= sum(w1)
          S31 	= sum(w1*zx)
          S41	= sum(w1*(zx^2))
          F1		= (1/k01)+S11-(S31/S21)
          f1		= (1/(k01^2)) + ((S21*S41)-(S31^2))/(S21^2)
          kmlex	= k01 + (F1/f1)
          while(abs(F1)>=err)	{
            k01	= kmlex
            S11	= mean(zx)
            w1	= x^k01	
            S21	= sum(w1)
            S31 	= sum(w1*zx)
            S41	= sum(w1*(zx^2))
            F1	= (1/k01)+S11-(S31/S21)
            f1	= (1/(k01^2)) + ((S21*S41)-(S31^2))/(S21^2)
            kmlex	= k01 + (F1/f1)
          }			
          areal[a]	= a2
          kreal[a]	= k2
          amlex[a] 	= ((1/n2)*S21)^(1/kmlex)
          kmlexx[a] 	= kmlex
          x1 	= rbind(x)
          xx1[a,]	=x1
        } 
      } 
      return(list(amlex,kmlexx,xx1,areal,kreal))
    }
    data1	= gen.x(k,n1,n2,a1,a2,k1,k2) 		#result of function gen.x
    aa	= data1[[1]]	#estimator for scale parameter (MLE)
    kk	= data1[[2]]	#estimator for shape parameter (MLE)
    xx1 	= data1[[3]]	#data
    areal 	= data1[[4]]	#scale parameter
    kreal 	= data1[[5]] 	#shape parameter
    for (d in 1:k)
    {									#find var(mui)
      MU_real[d]		= areal[d]*gamma(1+(1/kreal[d]))
      L1x.r		= (xx1[d,]/areal[d])^kreal[d] 
      L2x.r		= log(xx1[d,]/areal[d])
      L3x.r		= (n1*kreal[d]/(areal[d]^2)) -( (kreal[d]*(kreal[d]+1)/(areal[d]^2)) *sum( L1x.r))
      L4x.r		= - (n1*(kreal[d]^2)) - sum(L1x.r*(L2x.r^2))
      L5x.r		= -(n1/areal[d]) + (kreal[d]/areal[d])*sum(L1x.r*L2x.r) + ((1/areal[d])*sum(L1x.r))
      Lx.r			= matrix(c(-L3x.r,-L5x.r,-L5x.r,-L4x.r),2,2)	
      Linvx.r		= solve(Lx.r)
      A1x.r		= Linvx.r[2,2]					#var(k)
      A2x.r		= areal[d]*digamma(1+(1/kreal[d]))
      A3x.r		= Linvx.r[1,2]
      A4x.r		= gamma(1+(1/kreal[d]))
      A5x.r		= Linvx.r[1,1]					#var(a)
      term1x.r		= A1x.r*(A2x.r^2)
      term2x.r		= 2*A3x.r*A2x.r*A4x.r
      term3x.r		= A5x.r*(A4x.r^2)
      V.mhx.r	    	= term1x.r+term2x.r+term3x.r   
      V.mh.r[d]     		= V.mhx.r					#var(mui)
    } #end loop d
    
    ################################ Estimation of parameter of interest ################################
    
    phi = sum(MU_real/V.mh.r)/sum(1/V.mh.r)				#estimator of common mean
    
    #################################### Compute CI based on GCI ###################################
    
    for (j in 1:m)
    { 						#start loop j
      gen.star <- function(k,n1) # begin function 
      {						#start k,a star
        for (b in 1:k)
        { #start loop b
          x_x 	= rweibull(n1,1,1)
          zzx 	= log(x_x)
          zzbarx 	= mean(zzx)
          s_21  	= sum((zzx-zzbarx)^2)
          s_11	= sqrt(s_21)
          k_u1 	= (pi/sqrt(6))*(sqrt(n1-1))/s_11
          k.01	= k_u1
          S_11	= mean(zzx)
          ww1	= (x_x)^k.01	
          S_21	= sum(ww1)
          S_31  	= sum(ww1*zzx)
          S_41	= sum(ww1*(zzx^2))
          F21	= (1/k.01)+(S_11)-(S_31/S_21)
          f21	= (1/(k.01^2)) + ((S_21*S_41)-(S_31^2))/(S_21^2)
          k.starx	= k.01 + (F21/f21)
          while(abs(F21)>=err)	{
            k.01	= k.starx
            S_11	= mean(zzx)
            ww1	= (x_x)^k.01	
            S_21	= sum(ww1)
            S_31  	= sum(ww1*zzx)
            S_41	= sum(ww1*( zzx^2))
            F21	= (1/k.01)+(S_11)-(S_31/S_21)
            f21	= (1/(k.01^2)) + ((S_21*S_41)-(S_31^2))/(S_21^2)
            k.starx	= k.01 + (F21/f21)
          }	#end while
          a.starx[b] 	= ((1/n1)*S_21)^(1/k.starx)	
          k.starxx[b] 	= k.starx	
        }#end loop b
        return(list(a.starx, k.starxx)) 
      }#end function gen.star
      data2	= gen.star(k,n1) 			#result of function gen.x
      aa.star	= data2[[1]]			#a.star
      kk.star	= data2[[2]]			#k.star
      Gk 		= kk/kk.star
      Ga		= ((1/aa.star)^(kk.star/kk))*aa
      Rmu 	= (Ga*gamma(1+(1/Gk)))               	#R.mui
      a.calRvar	= cbind(Ga)			#Ga for calculate Rvar(mu)
      k.calRvar	= cbind(Gk)			#Gk for calculate Rvar(mu)
      
      for (d in 1:k)
      {
        L1x	= (xx1[d,]/a.calRvar[d,])^k.calRvar[d,]
        L2x	= log(xx1[d,]/a.calRvar[d,])
        L3x	= (n1*k.calRvar[d,]/(a.calRvar[d,]^2)) -( (k.calRvar[d,]*(k.calRvar[d,]+1)/(a.calRvar[d,]^2)) *sum( L1x))
        L4x	= - (n1*(k.calRvar[d,]^2)) - sum(L1x*(L2x^2))
        L5x	= -(n1/a.calRvar[d,]) + (k.calRvar[d,]/a.calRvar[d,])*sum(L1x*L2x) + ((1/a.calRvar[d,])*sum(L1x))
        Lx		= matrix(c(-L3x,-L5x,-L5x,-L4x),2,2)	
        Linvx	= solve(Lx)
        A1x	= Linvx[2,2]				
        A2x	= a.calRvar[d,]*digamma(1+(1/k.calRvar[d,]))
        A3x	= Linvx[1,2]
        A4x	= gamma(1+(1/k.calRvar[d,]))
        A5x	= Linvx[1,1]				
        term1x	= A1x*(A2x^2)
        term2x	= 2*A3x*A2x*A4x
        term3x	= A5x*(A4x^2)
        V.mhx[d]	= term1x+term2x+term3x        
      } #start loop d
      V.mhx.R 	= V.mhx					#Rvar.mu
      Rphi[j] 	= sum(Rmu/V.mhx.R)/sum(1/V.mhx.R) 	#Rmu
    } #start loop j
    
    #GCI
    L.CI1 		= quantile(Rphi,0.025,type=8) # compute lower based on GCI of phi
    U.CI1 		= quantile(Rphi,0.975,type=8) # compute upper based on GCI of phi
    CP.phiGCI[i] 		= ifelse(L.CI1<phi&&phi<U.CI1,1,0) # compute probability of CI based on GCI of phi
    Length.phiGCI[i] 	= U.CI1-L.CI1	# compute length of CI based on GCI of phi
    
    ############################## compute CI based on Adjusted MOVER ##############################    
    
    for (d in 1:k)
    {
      MU.hatx		= aa[d]*gamma(1+(1/kk[d]))
      L1x.m		= (xx1[d,]/aa[d])^kk[d] 
      L2x.m		= log(xx1[d,]/aa[d])
      L3x.m		= (n1*kk[d]/(aa[d]^2)) -( (kk[d]*(kk[d]+1)/(aa[d]^2)) *sum( L1x.m))
      L4x.m		= - (n1*(kk[d]^2)) - sum(L1x.m*(L2x.m^2))
      L5x.m		= -(n1/aa[d]) + (kk[d]/aa[d])*sum(L1x.m*L2x.m) + ((1/aa[d])*sum(L1x.m))
      Lx.m		= matrix(c(-L3x.m,-L5x.m,-L5x.m,-L4x.m),2,2)	
      Linvx.m		= solve(Lx.m)
      A1x.m		= Linvx.m[2,2]				#var(k)
      A2x.m		= aa[d]*digamma(1+(1/kk[d]))
      A3x.m		= Linvx.m[1,2]
      A4x.m		= gamma(1+(1/kk[d]))
      A5x.m		= Linvx.m[1,1]				#var(a)
      term1x.m		= A1x.m*(A2x.m^2)
      term2x.m		= 2*A3x.m*A2x.m*A4x.m
      term3x.m		= A5x.m*(A4x.m^2)
      V.mhx.m	    	= term1x.m+term2x.m+term3x.m   
      V.lmhx		= V.mhx.m/( MU.hatx^2)
      MU.hat.m[d] 		= MU.hatx
      L.CIx[d] 		= exp(log(MU.hatx)-(zvalue*sqrt(V.lmhx)))	#lower limit based on Wald interval
      U.CIx[d]		= exp(log(MU.hatx)+(zvalue*sqrt(V.lmhx))) #upper limit based on Wald interval
      V.mh[d]     		= V.mhx.m
      
      #Adjusted MOVER
      L.CI.M 		= phi.hat-zvalue*sqrt((1/(sum(((zvalue2)^2)/(MU.hat.m-L.CIx)^2))))	# compute lower based on AMVof phi
      U.CI.M  		= phi.hat+zvalue*sqrt((1/(sum(((zvalue2)^2)/(U.CIx-MU.hat.m)^2)))) 	# compute upper based on AMV of phi
      CP.phiMV[i]		= ifelse(L.CI.M<phi&&phi<U.CI.M,1,0) 	# compute probability of CI based on AMV of phi
      Length.phiMV[i] 	= U.CI.M-L.CI.M	# compute length of CI based on AMV of phi
      
      
      ################################## compute CI based on Bayesian #################################       
      
      kmlex.baye[1]  	= kk[d]
      amlex.baye[1]  	= aa[d]
      bmlex.baye[1]  	= (1/aa[d])^kk[d]
      V.mhx.b[1]     	= 1
      MU.hat.baye[1] 	= amlex.baye[1]*gamma(1+(1/kmlex.baye[1]))
      
      #set t for baye
      for (t in 2:T){
        ###Random walk metropolis
        scaleforbe.b		= 1/(z2+sum(xx1[d,]^kmlex.baye[t-1]))
        bmlex.baye[t]	= rgamma(1,shape=n1+v2,scale=scaleforbe.b)
        repeat{
          epx	= rnorm(1,0,sqrt(A1x.m))
          kstar.x	= kmlex.baye[t-1]+epx
          Atop.x	= ((kstar.x*bmlex.baye[t])^n1)*prod((xx1[d,])^(kstar.x-1))*exp((-bmlex.baye[t])*sum(xx1[d,]^kstar.x))*(v1^z1)/gamma(z1)*(kstar.x^(z1-1))*exp((-v1)*kstar.x)
          Abottom.x	= ((kmlex.baye[t-1]*bmlex.baye[t])^n1)*prod((xx1[d,])^(kmlex.baye[t-1]-1))*exp((-bmlex.baye[t])*sum(xx1[d,]^kmlex.baye[t-1]))*(v1^z1)/gamma(z1)*(kmlex.baye[t-1]^(z1-1))*exp((-v1)*kmlex.baye[t-1])
          Atop.x[is.nan(Atop.x)]=0
          Abottom.x[is.nan(Abottom.x)]=0
          if((Atop.x>0)&(Abottom.x>0))
            break} 
        Aofk.x  	= Atop.x/Abottom.x   		# compute Ak 
        cri.x	= min(1,Aofk.x)
        ux		= runif(1,0,1)			# compute u~U(0,1)
        if(ux<cri.x){kmlex.baye[t]=kstar.x}else{kmlex.baye[t]=kmlex.baye[t-1]}
        amlex.baye[t]	= (bmlex.baye[t])^(-1/kmlex.baye[t])
        MU.hat.baye[t]	= amlex.baye[t]*gamma(1+(1/kmlex.baye[t])) 	#Muhat for Bayesian
        
        L1x.b   		= (xx1[d,]/amlex.baye[t])^kmlex.baye[t]
        L2x.b		= log(xx1[d,]/amlex.baye[t])
        L3x.b		= (n1*kmlex.baye[t]/(amlex.baye[t]^2)) -((kmlex.baye[t]*(kmlex.baye[t]+1)/(amlex.baye[t]^2))*sum( L1x.b))
        L4x.b		= - (n1*(kmlex.baye[t]^2)) - sum(L1x.b*(L2x.b^2))
        L5x.b		= -(n1/amlex.baye[t]) + (kmlex.baye[t]/amlex.baye[t])*sum(L1x.b*L2x.b) + ((1/amlex.baye[t])*sum(L1x.b))
        Lx.b		= matrix(c(-L3x.b,-L5x.b,-L5x.b,-L4x.b),2,2)	
        Linvx.b		= solve(Lx.b)
        A1x.b		= Linvx.b[2,2]				
        A2x.b		= amlex.baye[t]*digamma(1+(1/kmlex.baye[t]))
        A3x.b		= Linvx.b[1,2]
        A4x.b		= gamma(1+(1/kmlex.baye[t]))
        A5x.b		= Linvx.b[1,1]				
        term1x.b		= A1x.b*(A2x.b^2)
        term2x.b		= 2*A3x.b*A2x.b*A4x.b
        term3x.b		= A5x.b*(A4x.b^2)
        V.mhx.b[t]		= term1x.b+term2x.b+term3x.b        #Rvar.mu.baye
      }#for t				
      MU.gibbs[d,]		= MU.hat.baye			#mu.baye
      V.mu.gibbs[d,]	= V.mhx.b 			#var.mu.baye
    } #end loop d
    
    phi.hat   	= sum(MU.hat.m/V.mh)/sum(1/V.mh)
    phi.hat.baye	= colSums(MU.gibbs/V.mu.gibbs)/colSums(1/V.mu.gibbs)
    phi.baye	= phi.hat.baye[-1:-1000]			#mu.baye (common mean)
    
    #Baye-MCMC.rwm
    Lphi.baye 		= quantile(phi.baye,0.025,type=8) # compute lower based on BAYE-MCMC of phi
    Uphi.baye 		= quantile(phi.baye,0.975,type=8) # compute upper based on BAYE-MCMC of phi
    CP.phiBAYE[i]	= ifelse(Lphi.baye<phi&&phi<Uphi.baye,1,0) # compute probability of CI based on BAYE-MCMC of phi
    Length.phiBAYE[i] 	= Uphi.baye-Lphi.baye # compute length of CI based on BAYE-MCMC of phi 
    
    #Baye-HPD.rwm
    HPD1.phi		= hdi(phi.baye,0.95) 
    Lphi.baye2		= HPD1.phi [1]	# compute lower based on BAYE-HPD of phi
    Uphi.baye2		= HPD1.phi [2]	# compute upper based on BAYE-HPD of phi
    CP.phiBAYE2[i] 	= ifelse(Lphi.baye2<phi&&phi<Uphi.baye2,1,0) # compute probability of CI based on BAYE-HPD of phi
    Length.phiBAYE2[i] 	= Uphi.baye2- Lphi.baye2 # compute length of CI based on BAYE-HPD of phi
  }	#end loop i
  
  ACP.phiGCI 	= mean(CP.phiGCI)		# compute coverage prob. based on GCI of phi
  Ex.length1 	= mean(Length.phiGCI)		# compute expected length based on GCI of phi
  ACP.phiMV 	= mean(CP.phiMV)		# compute coverage prob. based on GCI of phi
  Ex.length2 	= mean(Length.phiMV)		# compute expected length based on GCI of phi
  ACP.phiBAYE 	= mean(CP.phiBAYE)		# compute coverage prob. based on GCI of phi
  Ex.length3 	= mean(Length.phiBAYE)		# compute expected length based on GCI of phi
  ACP.phiBAYE2 = mean(CP.phiBAYE2)		# compute coverage prob. based on GCI of phi
  Ex.length4 	= mean(Length.phiBAYE2)	# compute expected length based on GCI of phi
  
  cat("COP.Phi|GCI		=",ACP.phiGCI,"\n")
  cat("Ex.length.Phi|GCI		=",Ex.length1,"\n")
  cat("COP.Phi|MOVER		=",ACP.phiMV,"\n")
  cat("Ex.length.Phi|MOVER	=",Ex.length2,"\n")
  cat("COP.Phi|BAYE		=",ACP.phiBAYE,"\n")
  cat("Ex.length.Phi|BAYE		=",Ex.length3,"\n")
  cat("COP.Phi|BAYE2		=",ACP.phiBAYE2,"\n")
  cat("Ex.length.Phi|BAYE2	=",Ex.length4,"\n")
}#end process