Appendices

Appendix A

Chapter-2

To show that Chakraborti and Li’s CI based on the pivotal quantity (LABEL:MVUE-PQ), it is enough to prove that $T$ is a one-to-one function of $t_m(z_p\sqrt{n})/\sqrt{n}$. The pivotal quantity in (LABEL:MVUE-PQ) is given by

$$T=\frac{{\widehat{\xi }}_p-{\xi }_p}{\sqrt{\widehat{\mathrm{var}}({\widehat{\xi }}_p)}}=\frac{\overline{X}+z_p{c}_{n}S-\mu -z_p\sigma }{\sqrt{\widehat{\mathrm{var}}({\widehat{\xi }}_p)}}.$$

Since the CI based on $T$ and the one based on $T^{\ast }=-T$ are the same, we now show that the CI based on $T^{\ast }$ is the same as the one in (LABEL:NCT-CI). Write

(1) $$T^{\ast }=-T=\frac{{\xi }_p-(\overline{X}+z_p{c}_{n}S)}{\frac{S}{\sqrt{n}}\sqrt{1+n{z}_p^2(c_n^2-1)}}=\frac{({\xi }_p-\overline{X})/S-z_p{c}_n}{\sqrt{1/n+z_p^2(c_n^2-1)}}.$$

Thus, using the result in (LABEL:nct-piv), we see that

(2) $$T^{\ast }\stackrel{d}{=}\frac{\frac{1}{\sqrt{n}}{t}_m(z_p\sqrt{n})-z_p{c}_n}{\sqrt{1/n+z_p^2(c_n^2-1)}}$$

and the percentiles of $T^{\ast }$ can be obtained from the percentiles of $t_m(z_p\sqrt{n})$. Specifically, for any , the 100$\alpha$ percentile of $T^{\ast }$ is given by

$$T_{\alpha }^{\ast }=\frac{\frac{1}{\sqrt{n}}{t}_{m;\alpha }(z_p\sqrt{n})-z_p{c}_n}{\sqrt{1/n+z_p^2(c_n^2-1)}}\mathrm{and}{T}_{\alpha }=-T_{\alpha }^{\ast }.$$

It can be easily checked that the 100$(1-2\alpha$)% CI for ${\xi }_p$ that can be obtained by solving the inequality $T_{\alpha /2}^{\ast }\le T^{\ast }\le T_{1-\alpha /2}^{\ast }$ for ${\xi }_p$, or equivalently, solving the inequality

$$\frac{\frac{1}{\sqrt{n}}{t}_{m;\alpha }(z_p\sqrt{n})-z_p{c}_n}{\sqrt{1/n+z_p^2(c_n^2-1)}}\le \frac{({\xi }_p-\overline{X})/S-z_p{c}_n}{\sqrt{1/n+z_p^2(c_n^2-1)}}\le \frac{\frac{1}{\sqrt{n}}{t}_{m;1-\alpha }(z_p\sqrt{n})-z_p{c}_n}{\sqrt{1/n+z_p^2(c_n^2-1)}}$$

for ${\xi }_p$, is the same as the classical CI in (LABEL:NCT-CI). Thus, we prove that the classical CI based on the noncentral $t$ distribution and Chakraborti and Li’s CI are the same.

Appendix B

NCT Confidence Intervals: We shall now show that the coverage probabilities and expected widths of CIs for $\mu +z_p\sigma$ are the same as those of the CIs for $\mu +z_{1-p}\sigma$. Recall that $t_m(\delta )=\frac{Z+\delta }{U}$, where $U^2\sim \frac{{\chi }_m^2}{m}$ independently of $Z\sim N(0,1)$. Using this representation, it can be readily checked that

(3) $$t_{m;\alpha }(-\delta )=-t_{m;1-\alpha }(\delta )\mathrm{and}{t}_{m;1-\alpha }(-\delta )=-t_{m;\alpha }(\delta ).$$

Using the above relations, it can be easily verified that the widths of the CIs

$$\begin{array}{c}(\overline{X}+t_{m;\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}},\overline{X}+t_{m;1-\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}})\hfill \\ \hfill \text{and}(\overline{X}+t_{m;\alpha }(z_{1-p}\sqrt{n})\frac{S}{\sqrt{n}},\overline{X}+t_{m;1-\alpha }(z_{1-p}\sqrt{n})\frac{S}{\sqrt{n}})\end{array}$$

are the same. The coverage probability of the NCT CI for $\mu +z_p\sigma$ can be expressed as

$$P\left(\overline{X}+t_{m;\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}\le \mu +z_p\sigma \le \overline{X}+t_{m;1-\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}\right).$$

Subtracting $\mu$ from all the terms and letting ${\overline{X}}^{\ast }=(\overline{X}-\mu )\sim N(0,{\sigma }^2/n)$, we can express the above coverage probability as

(4) $$P\left({\overline{X}}^{\ast }+t_{m;\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}\le z_p\sigma \le {\overline{X}}^{\ast }+t_{m;1-\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}\right).$$

Similarly, the coverage probability of the NCT CI for $\mu +z_{1-p}\sigma$ can be expressed as

$$P\left({\overline{X}}^{\ast }+t_{m;\alpha }(z_{1-p}\sqrt{n})\frac{S}{\sqrt{n}}\le z_{1-p}\sigma \le {\overline{X}}^{\ast }+t_{m;1-\alpha }(z_{1-p}\sqrt{n})\frac{S}{\sqrt{n}}\right),$$

or equivalently,

$$P\left(-{\overline{X}}^{\ast }-t_{m;\alpha }(-z_p\sqrt{n})\frac{S}{\sqrt{n}}\ge z_p\sigma \ge -{\overline{X}}^{\ast }-t_{m;1-\alpha }(-z_p\sqrt{n})\frac{S}{\sqrt{n}}\right),$$

Note that ${\overline{X}}^{\ast }=\overline{X}-\mu$ and $-{\overline{X}}^{\ast }=\mu -\overline{X}$ are identically distributed. Using this relation and the properties of the NCT percentiles in (3), the above probability can be expressed as

$$P\left({\overline{X}}^{\ast }+t_{m;1-\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}\ge z_p\sigma \ge {\overline{X}}^{\ast }+t_{m;\alpha }(z_p\sqrt{n})\frac{S}{\sqrt{n}}\right),$$

which is the same as the one in (4).

Normal Approximate Confidence Intervals: To show that the coverage probabilities and expected widths of normal approximate CIs for $\mu +z_p\sigma$ are the same as those of the CIs for $\mu +z_{1-p}\sigma$, we first note that

(5) $$h(z_{\alpha };n,1-p)=-h(z_{1-\alpha };n,p)\mathrm{and}h(z_{1-\alpha };n,1-p)=-h(z_{\alpha };n,p),$$

where $h(z_{\alpha };n,p)$ is given in (LABEL:hzalnp). The coverage probabilities of the normal approximate CI$(\overline{X}+h(z_{\alpha };n,p)S,\overline{X}+h(z_{1-\alpha };n,p)S)$ can be expressed as

(6) $$P\left({\overline{X}}^{\ast }+h(z_{\alpha };n,p)S\le z_p\sigma \le {\overline{X}}^{\ast }+h(z_{1-\alpha };n,p)S\right).$$

The coverage probability of the CI for $\mu +z_{1-p}\sigma$ is given by

$$P\left(\overline{X}+h(z_{\alpha };n,1-p)S\le \mu +z_{1-p}\sigma \le \overline{X}+h(z_{1-\alpha };n,1-p)S\right).$$

After subtracting $\mu$ from all the terms and noting that $z_{1-p}=z_p$, the above probability can be written as

$$P\left({\overline{X}}^{\ast }+h(z_{\alpha };n,1-p)S\le -z_p\sigma \le {\overline{X}}^{\ast }+h(z_{1-\alpha };n,1-p)S\right),$$

or equivalently,

$$P\left(-{\overline{X}}^{\ast }-h(z_{\alpha };n,1-p)S\ge z_p\sigma \ge -{\overline{X}}^{\ast }-h(z_{1-\alpha };n,1-p)S\right).$$

Using the result that ${\overline{X}}^{\ast }\stackrel{d}{=}-{\overline{X}}^{\ast }$ and the results in (5), it can be readily verified that the above probability is the same as the one in (6).

Appendix C

Chapter 2 R code

R code to compute the CIs for difference/ratio of two normal percentiles
************************************************************************
ci.norm.quantiles = function(n, xb, sq, p, cl){
  al = (1-cl)/2; m = n-1; zp = qnorm(p)
  sqn = sqrt(n); del = zp*sqn;
  s = sqrt(sq)
# CI for the Ratio; NCT
  t1 = qt(c(al,0.5,1-al), m[1], del[1])/sqn[1]
  t2 = qt(c(al,0.5,1-al), m[2], del[2])/sqn[2];
  Q1m = xb[1] + t1[2]*s[1];
  Q1l = xb[1] + t1[1]*s[1]; Q1u = xb[1] + t1[3]*s[1]
  Q2m = xb[2] + t2[2]*s[2];
  Q2l = xb[2] + t2[1]*s[2]; Q2u = xb[2] + t2[3]*s[2]
  prd = Q1m*Q2m;
  nr11 = Q2m^2-(Q2m-Q2u)^2; nr12 = Q1m^2-(Q1m-Q1l)^2
  Lowx = (prd- sqrt(prd^2-nr11*nr12))/nr11
  nr21 = Q2m^2-(Q2m-Q2l)^2; nr22 = Q1m^2-(Q1m-Q1u)^2
  Uppx = (prd + sqrt(prd^2-nr21*nr22))/nr21
  cat("NCT fiducial CI for Ratio:", c(Lowx, Uppx), "\n")
#CI for the Ratio; Normal Approximation
  cm = 1+0.25/m; zq = qnorm(c(al, 0.5, 1-al))
  l1 = (cm[1]*zp+zq[1]*sqrt(cm[1]^2/n[1]+
     +.5*(zp^2-zq[1]^2/n[1])/m[1]))/(cm[1]^2-zq[1]^2/2/m[1])
  u1 = (cm[1]*zp+zq[3]*sqrt(cm[1]^2/n[1]+
     +.5*(zp^2-zq[3]^2/n[1])/m[1]))/(cm[1]^2-zq[3]^2/2/m[1])
  md = zp/cm
  l2 = (cm[2]*zp+zq[1]*sqrt(cm[2]^2/n[2]+
     +.5*(zp^2-zq[1]^2/n[2])/m[2]))/(cm[2]^2-zq[1]^2/2/m[2])
  u2 = (cm[2]*zp+zq[3]*sqrt(cm[2]^2/n[2]+
     +.5*(zp^2-zq[3]^2/n[2])/m[2]))/(cm[2]^2-zq[3]^2/2/m[2])
  Ls = s[1]*md[1]-s[2]*md[2]-sqrt(sq[1]*(md[1]-l1)^2+sq[2]*(md[2]-u2)^2)
  Us = s[1]*md[1]-s[2]*md[2]+sqrt(sq[1]*(md[1]-u1)^2+sq[2]*(md[2]-l2)^2)
  xl1 = xb[1] + l1*s[1]; xu1 = xb[1] + u1*s[1]; xm1 = xb[1] + md[1]*s[1]
  xl2 = xb[2] + l2*s[2]; xu2 = xb[2] + u2*s[2]; xm2 = xb[2] + md[2]*s[2]
  prd = xm1*xm2; nr11 = xm2^2-(xm2-xu2)^2; nr12 = xm1^2-(xm1-xl1)^2
  nr21 = xm2^2-(xm2-xl2)^2; nr22 = xm1^2-(xm1-xu1)^2
  Lowx = (prd-sqrt(prd^2-nr11*nr12))/nr11
  Uppx = (prd+sqrt(prd^2-nr21*nr22))/nr21
  cat("Normal Apprx CI for Ratio:", c(Lowx, Uppx), "\n")
# CI for the Difference; NCT
  td1 = t1*sqn[1]; td2 = t2*sqn[2]
  xbd = xb[1]- xb[2]; w = s/sqn
  cent = w[1]*td1[2]-w[2]*td2[2]
  Dl = cent-sqrt(w[1]^2*(td1[2]-td1[1])^2+w[2]^2*(td2[2]-td2[3])^2)
  Du = cent+sqrt(w[1]^2*(td1[2]-td1[3])^2+w[2]^2*(td2[2]-td2[1])^2)
  Lowx = xbd + Dl
  Uppx = xbd + Du
  cat("NCT fiducial CI for Difference:", c(Lowx, Uppx), "\n")
# CI for the Difference; Normal Approx.
  l1 <- (cm[1] * zp +
            zq[1] * sqrt(cm[1]^2 / n[1] +
                .5 * (zp^2 - zq[1]^2 / n[1]) / m[1])
         ) / (cm[1]^2 - zq[1]^2 / 2 / m[1] )
  u1 = (cm[1]*zp +
            zq[3] * sqrt(cm[1]^2 / n[1] +
                .5 * (zp^2 - zq[3]^2 / n[1]) / m[1])
        ) / (cm[1]^2 - zq[3]^2 / 2 / m[1])
  l2 = (cm[2] * zp +
            zq[1] * sqrt(cm[2]^2 / n[2]+
                .5 * (zp^2 - zq[1]^2 / n[2]) / m[2])
        ) / (cm[2]^2 - zq[1]^2 / 2 / m[2])
  u2 = (cm[1] * zp +
            zq[3] * sqrt(cm[2]^2 / n[2] +
                .5 * (zp^2 - zq[3]^2 / n[2]) / m[2])
        ) / (cm[2]^2 -zq[3]^2 / 2 / m[2])
  Ls = s[1]*md[1]-s[2]*md[2]-sqrt(sq[1]*(md[1]-l1)^2+sq[2]*(md[2]-u2)^2)
  Us = s[1]*md[1]-s[2]*md[2]+sqrt(sq[1]*(md[1]-u1)^2+sq[2]*(md[2]-l2)^2)
  Low = xbd + Ls; Upp = xbd + Us
  cat("Normal Apprx CI for Difference:", c(Low, Upp), "\n")
}
#
n = c(107,100); xb = c(4840.3,7144.9); sq = c(2354470,2414487)
ci.norm.quantiles(n, xb, sq, .75,.95)
#
NCT fiducial CI for Ratio: 0.6693633 0.7689852
Normal Apprx CI for Ratio: 0.6694232 0.7686759
NCT fiducial CI for Difference: -2793.328 -1846.042
Normal Apprx CI for Difference: -2790.792 -1847.67

Appendix D

Chapter 3 R code

#Constrained mles and LRT
mles_q <- function(xbar, vsq, n_vals, p, reps, tol = 1e-7) {
  q_old <- 0
  diff <- Inf
  iterations <- 0
  zp <- qnorm(p)
  k <- length(xbar)
  sig_sq_old <- rep(1, k)

  calculate_terms <- function(sig, q) {
    # Update sigma values based on vsq and eta
    sig_sq_new <- ((zp*(xbar - q) + sqrt(zp^2 * (xbar - q)^2 +
                                           4*(vsq + (xbar-q)^2))) / 2)^2
    # Calculate the full term for eta
    denom_term <- n_vals / sig_sq_new
    sum_denom <- sum(denom_term)
    q_new <- sum((n_vals / sig_sq_new)* (xbar + zp * sqrt(sig_sq_new))
                 ) / sum_denom
    list(q_new = q_new, sig_sq_new = sig_sq_new)
  }
  while (diff > tol) {
    iterations <- iterations + 1
    terms <- calculate_terms(sig_sq_old, q_old)
    q_new <- terms$q_new
    sig_sq_new <- terms$sig_sq_new
    diff <- abs(q_old - q_new)
    # Update eta_old and sig_sq_old
    q_old <- q_new
    sig_sq_old <- sig_sq_new
  }

  sig_sq_new <- sig_sq_new
  q_new <- q_new
  Lam <- sum((n_vals * (vsq + (xbar + zp * sqrt(sig_sq_new) - q_new)^2)
              ) / sig_sq_new) - sum(n_vals * log(vsq / sig_sq_new)
                                    ) - sum(n_vals)
  # Return results silently
  results <- list(
    Q = q_new,
    sigma_sq = sig_sq_new,
    Lam = Lam
  )

}

#DK-MLRT
DK_MLRT <- function(xb, vsq, n, p, reps) {
  zp <- qnorm(p)
  k <- length(n)

  mles <- mles_q(xb, vsq, n, p, reps)
  Q_c <- mles[[1]]
  sig_sq_c <- mles[[2]]
  lam_0 <- mles[[3]]
  mu_c <- Q_c - zp * sqrt(sig_sq_c)
  lam_s <- numeric(reps)
  max_n <- max(n)
  sim_sample <- array(NA, dim = c(max_n, k, reps))

  for (j in 1:k) {
    sim_sample[1:n[j], j, ] <- replicate(reps,
                                         rnorm(n[j],
                                               mu_c[j],
                                               sqrt(sig_sq_c[j])))
  }

  for (i in 1:reps) {
    sample_matrix <- sim_sample[, , i]

    xb_s <- colMeans(sample_matrix, na.rm = TRUE)
    s_sq_s <- apply(sample_matrix, 2, var, na.rm = TRUE)
    vsq_s <- s_sq_s * (n - 1) / n

    lam_s[i] <- mles_q(xb_s, vsq_s, n, p, reps)[[3]]
  }

  # Compute test statistic and p-value
  ml <- mean(lam_s)
  var_l <- var(lam_s)
  v <- 2 * ml^2 / var_l
  mlrt <- v * lam_0 / ml
  p_value <- pchisq(mlrt, v, lower.tail = FALSE)

  return(list(
    "test stat:" = mlrt,
    "p-value:"   = p_value
  ))

}
n=c(16,22,9)
xb <-  c(62.39,72.6,70.13)
s <- c(17.99,23.52,19.67)
s_sq <- s^2
vsq<- (n-1)/n * s_sq
p=0.5
DK_MLRT(xb, vsq, n, p, 100000)
#$‘test stat:‘
#[1] 2.289939

#$‘p-value:‘
#[1] 0.3161829

Appendix E

Chapter 4 R code

cis.comm.quantile = function(method, n, xb, sq, p, cl){
  dfs = m =n-1; s = sqrt(sq); al = (1-cl)/2
  dels = sqrt(n)*qnorm(p)
  tmd = qt(0.5, m, dels);
  tml = qt(al, m, dels); tmu = qt(1-al, m, dels);
  w = (n/sq)/(sum(n/sq)); sqn = sq/n; rtsqn = sqrt(sqn)
  cent = sum(w*xb)
  cent2 = sum(w*rtsqn*tmd)
  lo = cent2-sqrt(sum(w^2*sqn*(tmd-tml)^2))
  up = cent2+sqrt(sum(w^2*sqn*(tmd-tmu)^2))
  Low1 = cent+lo; Upp1 = cent+up
  Lows = Low1; Upps = Upp1
#
  if(method != "Fiducial"){
  fnR = function(x){
    res = pvalfn(method,n, xb, sq, "R", p, x)-al
    return(res)}
  fnL = function(x){
    res = pvalfn(method,n, xb, sq, "L", p, x)-al
    return(res)}
  LowL = Low1 - .1; LowR = Low1 + .1
  while(fnR(LowL)*fnR(LowR) > 0){
    LowL = LowL - .1; LowR = LowR + .1}
  Lows = uniroot(fnR, c(LowL, LowR))[[1]]
  UppL = Upp1 - .1; UppR = Upp1 + .1
  while(fnL(UppL)*fnL(UppR) > 0){
    UppL = UppL - .1; UppR = UppR + .1}
  Upps = uniroot(fnL, c(UppL, UppR))[[1]]}
#
  cat(method, 100*cl,"%  CI for", 100*p,"th percentile: \n")
  print(c(Lows, Upps), 5)
}

pvalfn = function(method, n, xb, sq, side, p, xi0){
  k = length(n); N = sum(n); s = sqrt(sq)
  dfs = n-1;
  if(method == "Fisher"){
    del = -qnorm(p)*sqrt(n)
  tstat = sqrt(n)*(xb-xi0)/s
  if(side == "R"){
    p = 1-pt(tstat, dfs, del)}
  else{p = pt(tstat, dfs, del)}
  stat = -2*sum(log(p))
  pval = 1-pchisq(stat, 2*k)
  }
  if(method == "Inv Chisq"){
  tstat = sqrt(n)*(xb-xi0)/s
  del = -qnorm(p)*sqrt(n)
  if(side == "R"){
    p = 1-pt(tstat, dfs, del)}
  else{p = pt(tstat, dfs, del)}
  invC = sum(qchisq(1-p, n))
  pval = 1-pchisq(invC, N)
  }
  if(method == "Inv Normal"){
    ws = dfs/sum(dfs)
    dr = sqrt(sum(ws^2))
    del = -qnorm(p)*sqrt(n)
    tstat = sqrt(n)*(xb-xi0)/s
    if(side == "R"){
      p = 1-pt(tstat, dfs, del)}
    else{p = pt(tstat, dfs, del)}
    invN = sum(ws*qnorm(p))/dr
    pval = pnorm(invN)
  }
  return(pval)
}
n = c(8,12,14,8); sq = c(85.711,20.748,2.729,33.640)
xb = c(105.00,109.75,109.50,113.25)

n = c(16, 22, 9); sq = c(17.99^2, 23.52^2, 19.67^2);
xb = c(62.39, 72.60, 70.13)
cis.comm.quantile("Fiducial",n, xb, sq, .5, .95)