#############################################################################################

#R-code for constructing GCI, AGCI, BPCI, FCI and F- HPDCI 


#############################################################################################

library(statmod)
library(MASS)
library(EnvStats)
library(HDInterval)


CI.SINGIG=function(M,m,n,mu,lambda)	
{						# Start
  
  starttime = proc.time()
  
  B=1000					# number of bootstrap sample
  
  CP.GCI = rep(0,M)	      # coverage probability of CI based on GCI
  Length.GCI = rep(0,M)				# average length of CI based on GCI
  R_theta_GCI = rep(0,m)
  
  CP.AGCI = rep(0,M)				# coverage probability of CI based on AGCI
  Length.AGCI = rep(0,M)				# average length of CI based on AGCI
  R_theta_AGCI = rep(0,m)
  
  CP.BPCI = rep(0,M)				# coverage probability of CI based on BPCI
  Length.BPCI = rep(0,M)				# average length of CI based on BPCI
  theta.star = rep(0,B)
  
  CP.FCI = rep(0,M)				# coverage probability of CI based on Fiducial
  Length.FCI = rep(0,M)				# average length of CI based on Fiducial
  
  CP.FHPD = rep(0,M)				# coverage probability of CI based on Fiducial HPD
  Length.FHPD = rep(0,M)				# average length of CI based on Fiducial HPD
  
  mu.t = c()
  lambda.t = c()
  mu.gibbs = c()
  lambda.gibbs = c()
  
  muhat.t = c()
  lambdahat.t = c()
  ################################ Estimation of parameter of interest ################################
  theta = sqrt(mu/lambda)				# single CV
  #############################Start Iteration##########################
  for (i in 1:M)
  {						# begin loop i
    x 	      = rinvgauss(n, mu, lambda)		# generate x           	  
    muhat 	  = mean(x) 		          	# compute mean of x        		 
    lambdahat_inverse = (1/n)*sum((1/x)-(1/muhat))	 
    lambdahat = 1/lambdahat_inverse		# compute scale of x
    
    ################################Compute CI based on GCI& AGCI################################
    for (j in 1:m)
    {						# begin loop j
      xbar = rinvgauss(1, mu, n*lambda)
      A=rchisq(1,n-1)
      Z=rnorm(1,0,1)
      V=A/(n*lambda)
      R_lambda = A/(n*lambdahat_inverse)		# compute GPQ for scale of x
      R_mu_GCI =xbar/abs(1+Z*(sqrt(xbar/(n*R_lambda)))) # compute GPQ for mean of x
      R_theta_GCI[j]=sqrt(R_mu_GCI/R_lambda) # compute GPQ for single CV
      
      W = (sqrt(n-1)*(xbar/mu-1))/sqrt(V*xbar)
      R_mu_AGCI = xbar/max(0,1+(W*sqrt((V*xbar)/(n-1))))	# compute GPQ for mean of x
      R_theta_AGCI[j] = sqrt(R_mu_AGCI/R_lambda)		# compute GPQ for single CV
      
    }						# end loop j 
    L.GCI  <- quantile( R_theta_GCI,0.025,type=8)  	# compute lower based on GCI of theta     
    U.GCI  <- quantile( R_theta_GCI,0.975,type=8)	# compute upper based on GCI of theta
    
    L.AGCI  <- quantile( R_theta_AGCI,0.025,type=8)  	# compute lower based on AGCI of theta     
    U.AGCI  <- quantile( R_theta_AGCI,0.975,type=8)	# compute upper based on AGCI of theta
    
    CP.GCI[i] 	<- ifelse(  L.GCI <theta&&theta<U.GCI ,1,0)	# compute probability of CI based on GCI of theta
    Length.GCI[i] <- U.GCI  -  L.GCI 				 # compute length of CI based on GCI of theta 
    
    CP.AGCI[i] 	<- ifelse(  L.AGCI <theta&&theta<U.AGCI ,1,0)	# compute probability of CI based on AGCI of theta
    Length.AGCI[i] <- U.AGCI  -  L.AGCI 				 # compute length of CI based on AGCI of theta 
    
    ################################Compute CI based on BPCI################################
    
    for (k in 1:B)
    {							# begin loop k
      x.star= sample(x,n,replace=T)				# generate bootstrap sample of x
      xbar.star= mean(x.star)				# compute mean of bootstrap sample of x
      lambdahat_inv_bt = (1/n) * sum((1/x.star)-(1/xbar.star))	
      lambdahat_bt = 1/lambdahat_inv_bt			# compute scale of bootstrap sample of x
      theta.star[k]= sqrt(xbar.star/lambdahat_bt)		# compute single CV of bootstrap sample	
    } 							# end loop k
    
    
    L.BPCI <- quantile(theta.star,0.025,type=8)       		# compute lower based on PBCI of theta
    U.BPCI <- quantile(theta.star,0.975,type=8)		# compute lower based on PBCI of theta
    CP.BPCI[i] <- ifelse( L.BPCI <theta&&theta<  U.BPCI,1,0)	# compute probability of CI based on PBCI of theta
    Length.BPCI[i] <- U.BPCI-L.BPCI  			# compute length of CI based on PBCI of theta
    
    ################################Compute CI based on FCI################################
    muhat.t[1] =   muhat 					# MLE of  mean x
    lambdahat.t[1] = lambdahat				# MLE of  scale x
    
    ### set t for Gibbs ####
    for (t in 2:20000)
    {							# begin loop t
      muhat.t[t]=muhat.t[t-1] 				
      lambdahat.t[t] =lambdahat.t[t-1]
      lambda.t[t]  = (lambdahat.t[t]/n)*rchisq(1,n-1)
      mu.t[t] = rinvgauss(1, muhat.t[t], n*lambdahat.t[t])	
    }							# end loop t
    
    mu.gibbs = mu.t
    lambda.gibbs = lambda.t
    
    ## Gibbs sampling for fiducial distribution#
    
    thetahat.FCI = sqrt(mu.gibbs/lambda.gibbs)		# compute CV base on MCMC
    theta.FCI =  thetahat.FCI[-1:-1000]			# burn in 1000
    
    L.CI.FCI <- quantile(theta.FCI,0.025,type=8)       		# compute lower based on Fiducial of theta
    U.CI.FCI <- quantile(theta.FCI,0.975,type=8)		# compute upper based on Fiducial of theta
    
    CP.FCI[i] <- ifelse(L.CI.FCI<theta&&theta<  U.CI.FCI,1,0)  # compute probability of CI based on Fiducial of      theta
    Length.FCI[i] <-  U.CI.FCI-L.CI.FCI			# compute length of CI based on Fiducial of theta
    
    ################################Compute CI based on F-HPDCI################################
    FHPD = hdi(theta.FCI,0.95) 
    L.FHPD = FHPD[1]					# compute lower based on Fiducial HPD of theta
    U.FHPD = FHPD[2]					# compute upper based on Fiducial HPD of theta
    CP.FHPD[i] = ifelse(L.FHPD<theta&&theta<U.FHPD,1,0)	# compute probability of CI based on Fiducial HPD of theta
    Length.FHPD[i] <-  U.FHPD-L.FHPD			# compute length of CI based on Fiducial HPD of theta
    
  } # end loop i
  
  elapseTime = proc.time()-starttime
  
  ACP.GCI <- mean(CP.GCI)	 		# compute coverage probability based on GCI of theta   
  Ex.length.GCI <- mean(Length.GCI)	   	# compute expected length based on GCI of theta   
  cat("CP|GCI=", ACP.GCI,"Length|GCI=", Ex.length.GCI ,"\n")
  
  ACP.AGCI <- mean(CP.AGCI)	 		# compute coverage probability based on AGCI of theta   
  Ex.length.AGCI <- mean(Length.AGCI)	   	# compute expected length based on AGCI of theta   
  cat("CP|AGCI=", ACP.AGCI,"Length|AGCI=", Ex.length.AGCI ,"\n")
  
  ACP.BPCI <- mean(CP.BPCI)   			# compute coverage probability based on BPCI of theta   
  Ex.length.BPCI <- mean(Length.BPCI) 		# compute expected length based on BPCI of theta     
  cat("CP|BPCI=", ACP.BPCI,"Length|BPCI=", Ex.length.BPCI,"\n")
  
  ACP.FCI <- mean(CP.FCI)   			# compute coverage probability based on Fiducial of theta   
  Ex.length.FCI <- mean(Length.FCI)   		# compute expected length based on Fiducial of theta     
  cat("CP|FCI=", ACP.FCI,"Length|FCI=", Ex.length.FCI,"\n")
  
  
  ACP.FHPD <- mean(CP.FHPD)   			# compute coverage probability based on Fiducial HPD of theta   
  Ex.length.FHPD <- mean(Length.FHPD)   		# compute expected length based on Fiducial HPD of theta     
  cat("CP|FHPD=", ACP.FHPD,"Length|FHPD=", Ex.length.FHPD,"\n")
  
  cat("Time =", elapseTime["elapsed"])
  
} #end process