Approximate and generalized pivotal quantities for deriving confidence intervals for the offset between two clocks

Statistical Methodology 6 (2009) 97–107

Contents lists available at ScienceDirect

Statistical Methodology journal homepage: www.elsevier.com/locate/stamet

Approximate and generalized pivotal quantities for deriving confidence intervals for the offset between two clocks Jun Li, Daniel R. Jeske ∗ , Jeff Pettyjohn Department of Statistics, University of California, Riverside, CA 92521, United States

article

info

Article history: Received 26 May 2007 Received in revised form 22 February 2008 Accepted 10 April 2008 Keywords: Relative profile likelihood Non-regular case Generalized confidence intervals Bootstrap calibration Clock offset

a b s t r a c t Development of algorithms that estimate the offset between two clocks has received a lot of attention, with the motivating force being data networking applications that require synchronous communication protocols. Recently, statistical modeling techniques have been used to develop improved estimation algorithms with the focus being obtaining robust estimators in terms of mean squared error. In this paper, we extend the use of statistical modeling techniques to address the construction of confidence intervals for the offset parameter. We consider the case where the distributions of network delays are members of a scale family. Our results include an asymptotic confidence interval and a generalized confidence interval in the sense of [S. Weerahandi, Generalized confidence intervals, Journal of the American Statistical Association 88 (1993) 899–905. Correction in vol. 89, p. 726, 1994]. We compare and contrast the two approaches for obtaining a confidence interval, and illustrate specific applications using exponential, Rayleigh and heavy-tailed Weibull network delays as concrete examples. © 2008 Elsevier B.V. All rights reserved.

1. Introduction A widely used technique for clock synchronization involves each pair of clocks exchanging timing messages with each other. A recent application to wireless sensor networks is described in [1]. Other motivating applications and a formulation of the problem within a statistical modeling context, which we now briefly review, are given in [2]. Fig. 1 shows the ith exchange between two clocks, A and B, that

∗ Corresponding address: Room 2605 Statistics-Computer Building, Department of Statistics, University of California, Riverside CA 92511, United States. Tel.: +1 951 827 3014; fax: +1 951 827 3286. E-mail addresses: [email protected] (J. Li), [email protected] (D.R. Jeske), [email protected] (J. Pettyjohn). 1572-3127/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.stamet.2008.04.002

98

J. Li et al. / Statistical Methodology 6 (2009) 97–107

Fig. 1. Timing message exchange.

begins with A sending a message to B which includes a time stamp Ti0 indicating the time as known at A when the message is sent. Immediately upon receipt of this message, B puts a time stamp Ti1 on the received message. Just before B sends the message back to A (which is not necessarily immediately), it places another time stamp Ti2 on the message. When A receives the returned message from B, it records another time stamp Ti3 on the message. After n such exchanges, A has the sequence of observations {Ti0 , Ti1 , Ti2 , Ti3 }ni=1 from which to synchronize itself with B. The sojourn times for the A → B and B → A timing messages are Xi = Ti1 − Ti0 and Yi = Ti3 − Ti2 , respectively. Define θ to be the (unknown) offset of clock B, relative to clock A. [We note that offset is frequently assumed to be constant, an assumption that is justified if the clocks are synchronized on a relatively frequent basis.] It follows that if at a given instant clock B shows time tB and clock A shows time tA , then tB = tA + θ. Apart from the offset, two other components of the message sojourn times are propagation delay and network delay. Propagation delay corresponds to the amount of time required for a timing message to travel across the transmission medium (e.g., optical fiber) between A and B if there were no other intervening delays. In general, the propagation delay is fixed and quite small, on the order of nanoseconds, for example. It is customarily assumed that the (unknown) propagation delays in the two transmission directions A → B and B → A are the same, and we denote the common value by d. Network delay arises due to congestion in the network that requires messages to sit in one or more queues before being granted the resources needed to continue their journey. Network delay can be on the order of milliseconds, or even seconds, and is the variable component of sojourn times. We let eAB i and eBA denote network delays associated with the A → B and B → A timing messages, respectively. It i follows from the preceding definitions, that Xi = d + θ + eAB and Yi = d − θ + eBA i i . The usual assumption n BA n is that the {eAB } and { e } sequences are independent and identically distributed according to i i i =1 i =1 distributions FAB (u) and FBA (u), respectively. A considerable amount of attention has been given to point estimation of θ under varying assumptions for the network delay distributions. One of the earliest estimators was proposed in [3] and is θˆ = (X1:n − Y1:n )/2, where X1:n and Y1:n denote the minimum of {Xi }ni=1 and {Yi }ni=1 , respectively. It was shown in [4] that θˆ is the maximum likelihood estimator (MLE) of θ under the assumption that the network delays follow exponential distributions. The robustness of alternative estimators of θ, including bootstrap bias-corrected versions of θˆ , using mean squared error as the criterion was investigated in references [2,5,6]. Until now, other inference procedures concerning θ such as construction of confidence intervals or hypothesis tests have not appeared in the literature. In this paper, we show that if one-parameter scale families are used to model network delays it is possible to obtain satisfactory confidence intervals for θ using an approximate pivotal quantity (APQ). Furthermore, we construct a generalized pivotal quantity (GPQ) in the sense of [7] that leads to an exact generalized confidence interval for θ. It is well known that generalized confidence intervals do not, in general, have frequentist coverage probabilities that are equal to the nominal value. However, there is growing practical experience with generalized confidence intervals that suggests they have frequentist coverage probabilities that are adequately close to the nominal value, and we will show that this indeed turns out to be the case for our application. The rest of this paper is organized as follows. In Section 2, we develop an approximate confidence interval for θ based upon the APQ and also alternative generalized confidence intervals based on different GPQs. The important features of each procedure are discussed. In Section 3 we illustrate and contrast the procedures by applying them to three special cases where network delays have exponential, Rayleigh and heavy-tailed Weibull distributions. Included in Section 3 are results from

J. Li et al. / Statistical Methodology 6 (2009) 97–107

99

a simulation study that compares the confidence interval procedures with respect to the frequentist interpretation of coverage probability and also in terms of the expected width of the interval. Section 4 concludes the paper with a discussion. 2. Confidence intervals for θ 2.1. Classical approach 2.1.1. APQ for exponential network delays We begin by considering the case of exponential network delays, since a result available in this case points the way toward obtaining an APQ for the larger class of scale families. In this setting, we will derive the relative profile likelihood function for θ, say Λ(θ), and show that it is an asymptotic pivot in the sense its limiting distribution does not depend on any unknown parameters. That Λ(θ) is an asymptotic pivot, combined with its easily invertible form, provides a path for obtaining an asymptotic confidence interval for θ. We note here, as a portent of things to come later, that the asymptotic distribution of −2 log Λ(θ) is not the anticipated chi-square distribution with one degree of freedom. Consider now a random sample of n timing message exchanges {(Xi , Yi )}ni=1 . Let X and Y denote the −1 −1 ˜ direction, vectors of {Xi }ni=1 and {Yi }ni=1 , and let γAB and γBA denote the mean network delays˜ in each respectively. The full likelihood function for the four parameters of the problem is ( ) n n X X n n L(θ, d, γAB , γBA ) ∝ γAB γBA exp −γAB (Xi − θ − d) − γBA (Yi + θ − d) , i=1

i=1

defined on {(d, θ, γAB , γBA ) : X(1) ≥ d + θ, Y(1) ≥ d − θ, γAB > 0, γBA > 0}. Fixing θ, it is straightforward to show that the conditional MLEs for the nuisance parameters (d, γAB , γBA ) are d˜ (θ) = min{X1:n − θ, Y1:n + θ} n

γ˜ AB (θ) =

1 P

[Xi − d˜ (θ) − θ]

i=n

γ˜ BA (θ) =

n 1 P

[Yi − d˜ (θ) − θ]

.

i=n

Substituting these expressions into the full likelihood results, we obtain the profile likelihood for θ, " ! !#n 1 1 Lp (θ) ∝ , −Y1:n ≤ θ ≤ X1:n . (1) X − d˜ (θ) − θ

Y − d˜ (θ) + θ

Note that the constraint −Y1:n ≤ θ ≤ X1:n follows from the requirement d˜ (θ) ≥ 0. Recalling that θˆ = (X1:n − Y1:n )/2, the MLE of θ in this setting, an alternative expression for (1) is  n  1 1   , −Y1:n ≤ θ ≤ θˆ  Lp (θ) ∝  Y − Y1:n   X − Y1:n − 2θ n (2) 1 1    , θˆ ≤ θ ≤ X1:n . X − X1:n Y − X1:n + 2θ

ˆ , is The relative profile likelihood for θ, defined as Λ(θ) = Lp (θ)/Lp (θ)  !n  X − X1:n   , −Y1:n ≤ θ ≤ θˆ   X − Y1:n − 2θ !n Λ(θ) =  Y − Y1:n    , θˆ ≤ θ ≤ X1:n .  Y − X1:n + 2θ

(3)

J. Li et al. / Statistical Methodology 6 (2009) 97–107

100

−1 Theorem 1. If the A → B network delays have an exponential distribution with mean γAB and the B → −1 A network delays have an exponential distribution with mean γBA , then the asymptotic distribution of − log Λ(θ) is exponential with mean equal to 1.

Proof. By direct substitution and straightforward manipulation of the generating equations for

{(Xi , Yi )}ni=1 , it follows from (3) that Λ(θ) ,

1−

!n       nW1:n nW1:n − ρnV1:n n nW1:n I ≤ρ + 1− I >ρ , nV1:n nV1:n n(V − W1:n /ρ) n(W − ρV1:n )

nV1:n − nW1:n /ρ

where {Vi }ni=1 and {Wi }ni=1 are independent sequences of random variables from exponential distribution with mean equal to 1, ρ = γBA /γAB and the symbol , signifies “equal in distribution to”. Since E1 = nV1:n and E2 = nW1:n are again independently distributed as standard exponential random variables, we can write !n    n   E2 E2 − ρE1 E2 E1 − E2 /ρ I ≤ρ + 1− I >ρ . Λ(θ) , 1 − n(V − W1:n /ρ)

E1

n(W − ρV1:n )

E1

Using the facts that V and W converge almost surely to 1, V1:n and W1:n converge in probability to 0, and the uniform convergence of the limiting form of the exponential function, it follows that d

as n → ∞, Λ(θ) → exp(Z /ρ)I(Z ≤ 0) + exp(−Z )I(Z > 0), where Z = E2 − ρE1 . Equivalently, d

− log Λ(θ) →(−Z /ρ)I(Z ≤ 0) + ZI(Z > 0). The distribution of the right-hand side is a standard exponential distribution from the fact that the conditional distributions of Z and −Z /ρ, given Z > 0 and Z ≤ 0, respectively, are easily shown to be both exponential with mean 1.  Note that Theorem 1 equivalently implies that −2 log Λ(θ) has an asymptotic chi-square distribution with 2 degrees of freedom, rather than the “expected” one degree of freedom. Failure of the classic limiting distribution of −2 log Λ(θ) is due to our non-regular context. Another alternative characterization implied by Theorem 1 is that Λ(θ) has an asymptotic uniform distribution on (0, 1), and this characterization is the most direct path toward establishing that an asymptotic 100(1 − α)% confidence interval for θ is [L(X , Y ; α), U (X , Y ; α)], where ˜ ˜ ˜ ˜  ( )  X − X1:n 1 ˆ L(X , Y ; α) = max θ − − 1 , −Y1:n 2 α1/n ˜ ˜ (4) ) (   Y − Y1:n 1 ˆ − 1 , −X1:n . U (X , Y ; α) = min θ + 2 α1/n ˜ ˜ 2.1.2. APQ for one-parameter scale-family network delays We now drop the assumption of exponentially distributed network delays and consider an arbitrary one-parameter scale family of distributions for the network delays. As in Section 2.1.1, we might similarly try to derive asymptotic pivotal quantities based on the respective relative profile likelihood function under individual network delay distributions. However, in many cases, it will not be possible to obtain a closed form for the MLE of θ and moreover, as we can see for the exponential case, its asymptotic normality cannot be taken for granted and the usual asymptotic chi-square distribution associated with the relative profile likelihood function may not hold. In addition, the relative profile likelihood is usually a complicated function of θ and whether or not its inversion produces an interval will depend on the family of distributions that is assumed for the network delays. Due to these difficulties and the relatively easy inversion of Λ(θ) defined in (3), we will subsequently view Λ(θ) in (3) as a proposed approximate pivotal quantity for θ and then characterize its limiting distribution when the network delay distributions are members of a one-parameter scale family. Hereon we assume that the distribution functions of the A → B and B → A network delays can be written as FAB (u) = H(γAB u) and FBA (u) = H(γBA u), respectively, where H(u) is a continuous cumulative distribution function corresponding to a non-negative random variable. In addition, we assume the Gnedenko condition [8] holds for H. That is, there exists a τ > 0 such that

J. Li et al. / Statistical Methodology 6 (2009) 97–107

101

limu→0+ {H(σ u)/H(u)} = σ τ , for all σ . The Gnedenko condition implies the asymptotic distribution of bn S1:n , where bn = 1/H−1 (n−1 ) and S1:n is the minimum of n random observations from H, is a standardized Weibull distribution with shape parameter equal to τ . Reference [9] gives an accessible summary of this result from Gnedenko’s original article which is written in French. We note that the standardized member of a number of scale distribution families satisfies the Gnedenko condition, including exponential and uniform that have τ = 1, Rayleigh and half-normal that have τ = 2, and, for a specified shape parameter, Weibull and gamma that both have τ equal to the shape parameter. Theorem 2 extends what can be said about the limiting distribution of Λ(θ) under this more general model for the network delay distributions. Theorem 2. If the A → B and B → A network delay distributions have cumulative distribution functions of the forms H(γAB u) and H(γBA u), respectively, where H(u) is a continuous cumulative distribution function for a non-negative random variable and satisfies the Gnedenko condition for some τ > 0, then d

−bn log Λ(θ)/n → −Z /(ξρ)I(Z ≤ 0) + (Z /ξ)I(Z > 0), where Z = E2 − ρE1 with E1 and E2 independently and identicallyR distributed as Weibull with shape ∞ parameter equal to τ and scale parameter equal to 1, and where ξ = 0 (1 − H(u))du is the mean of H. Proof. Our starting point is again (3) which we now manipulate to express as !bn   bn bn V1:n − bn W1:n /ρ bn W1:n ≤ρ [Λ(θ)] n , 1 − I bn (V − W1:n /ρ)



+

bn V1:n bn   bn W1:n bn W1:n − ρbn V1:n I >ρ , 1− bn V1:n bn (W − ρV1:n )

where now {Vi }ni=1 and {Wi }ni=1 are independent sequences of random variables from H(u). We now use the Gnedenko condition, the fact that V and W converge almost surely to ξ, the fact that V1:n and W1:n converge in probability to 0, and the uniform convergence of the limiting form of the exponential function to write bn

d

[Λ(θ)] n → exp{−Z /(ξρ)}I(Z ≤ 0) + exp{−Z /ξ}I(Z > 0), where Z = E2 − ρE1 and E1 and E2 are independently and identically distributed as Weibull with shape parameter equal to τ and scale parameter equal to 1. Equivalently, d

−bn log Λ(θ)/n → −Z /(ξρ)I(Z ≤ 0) + (Z /ξ)I(Z > 0).  Remark. In the exponential case, H(u) = 1 − exp(−u) from which it follows that bn ∼ n and the conclusion in Theorem 2 can be seen to be an intermediate result in the proof of Theorem 1. That we were able to go further in Theorem 1 owes to the fact that Z has a tractable (asymmetric) Laplace distribution in the exponential case. More generally, the distribution of Z is not a particularly tractable one, being the difference of independent, though not identically distributed Weibull random variables. Nevertheless, Theorem 2 shows that a limiting distribution for −bn log Λ(θ)/n exists. While the limiting distribution of −bn log Λ(θ)/n does not have a convenient closed form, it is a simple matter to estimate its percentiles for a given value of ρ by simulating values of −Z /(ξρ)I(Z ≤ 0) + (Z /ξ)I(Z > 0). Denote the upper-α percentile of this limiting distribution by zα (ρ). Then if ρ were known, an asymptotic 100(1 − α)% confidence interval for θ would be {θ : −bn log Λ(θ)/n ≤ zα (ρ)} which can be expressed as [L(X , Y ; ρ, α), U (X , Y ; ρ, α)], where ( ˜ ˜ )  ˜ ˜   L(X , Y ; ρ, α) = max θˆ −

˜ ˜

(

U (X , Y ; ρ, α) = min θˆ +

˜ ˜

X − X1:n

2 Y − Y1:n

2

zα (ρ)

exp 



exp

− 1 , −Y1:n

bn zα (ρ) bn





)

(5)

− 1 , −X1:n .

In the exponential case, (5) coincides with (4) since bn = n and zα (ρ) = − log α. More generally, in order to use (5) the unknown ρ can be replaced by a consistent estimator. One choice would

J. Li et al. / Statistical Methodology 6 (2009) 97–107

102

Fig. 2. Bootstrap calibration algorithm.

be the unconditional MLE of ρ, but a computationally simpler approach is motivated by the fact the conditional MLEs of γAB and γBA , given (d, θ), frequently have relatively simple forms. Denoting ˆ γˆ AB (dˆ , θ) ˆ is a consistent the conditional MLEs by γˆ AB (d, θ) and γˆ AB (d, θ), respectively, ρˆ = γˆ BA (dˆ , θ)/ ˆ ˆ estimator of ρ, where d = (X1:n + Y1:n )/2 and as previously defined θ = (X1:n − Y1:n )/2. In the following, ˆ α), U (X , Y ; ρ, ˆ α)] as the APQ interval. we refer to the interval [L(X , Y ; ρ,

˜ ˜

˜ ˜

2.1.3. Calibration While the APQ interval has the correct asymptotic coverage for θ, in small samples the coverage may not be satisfactorily close to the nominal level of 100(1 − α)%. A calibration step may be necessary to tune the actual coverage probability of the APQ interval toward the nominal level. For this purpose, we use the (parametric) bootstrap calibration algorithm (see, for example, Reference [10]) shown in Fig. 2. The key step in the algorithm is step (3) which generates bootstrap samples {(Xi∗ , Yi∗ )}ni=1 of the twoway sojourn times using model parameter estimates obtained from the data sample {(Xi , Yi )}ni=1 . Each of the bootstrap samples is then used to compute B = 1000 confidence intervals I(X ∗ , Y ∗ ; δ), and each of those is checked for inclusion of θˆ computed from the data sample. The fraction˜ of˜the I(X ∗ , Y ∗ ; δ) confidence sets that include θˆ is the bootstrap estimate of the probability that the interval I(X˜ ∗ , Y˜ ∗ ; δ) ˜ ˜ and covers θ. When this fraction is satisfactorily close to the nominal value 1 − α, the algorithm stops the current value of δ is used as the calibration level in step (5) where the ultimate confidence interval for θ is computed. In the following, we refer to this ultimate confidence interval as the calibrated APQ interval. We note that the suggested choice of B = 1000 and ∆ = 0.01 work well for most practical applications.

J. Li et al. / Statistical Methodology 6 (2009) 97–107

103

2.2. Using a GPQ Weerahandi extended the notion of pivotal quantities to GPQs [7] and defined, through inversion of the GPQ, the concept of a generalized confidence set. The subtle difference in how the coverage probability of a generalized confidence set is interpreted compared to the frequentist interpretation is more clearly discussed in Chapter 6 of [11] where it is also mentioned that practitioners can use generalized confidence sets to obtain very good approximate confidence sets in the frequentist sense. As is the case with classical pivotal quantities, GPQs are not unique and the properties such as the expected length of the generalized confidence interval they produce can vary. In a general setting, a GPQ can be defined as follows. Let z be the realized value of a data vector ˜ θ the parameter of primary interest Z , and let ω = (θ, γ 0 )0 be the vector of model parameters with ˜ ˜ ˜ γ a possibly and vector-valued nuisance parameter. A function R = r(Z , z; ω) of Z , z, and ω is called a ˜ ˜ ˜ ˜ distribution ˜ generalized pivotal quantity for θ if it satisfies two conditions: (1) R has˜a ˜probability that is free of unknown parameters, and (2) the observed value of R, defined as R = r(z, z; ω), depends on ω ˜ ˜ ˜ will in general ˜ only through θ. While R has a distribution free of unknown parameters, its distribution depend on z. Denote Cα (z) to be a subset of the sample space of R which satisfies Pr[R ∈ Cα (z)] = 1 − α. ˜ for θ. The set {θ :˜r(z, z; ω) ∈ C˜α (z)} is the defined to be a 100(1 − α)% generalized confidence set ˜ ˜ how ˜ a class ˜ of GPQs can be constructed for our application when the class of network We now show delay distributions is a one-parameter scale family. Define Z = (X 0 , Y 0 )0 and let z = (x0 , y0 )0 denote its ˜ ˜γAB˜ and γBA are ˜ the˜respective ˜ observed value. With ω = (θ, γ 0 )0 , let γ = (d, γAB , γBA )0 where scale ˜ ˜ ˜ parameters associated with the A → B and B → A directions. Our class of GPQs is of the form R = r(Z , z; ω)

˜ ˜ ˜

=

x1:n − y1:n

2

−

1



h(x)

h(y)



˜ (X1:n − d − θ) − ˜ (Y1:n − d + θ) 2 h(X ) h(Y ) ˜ ˜ n

(6)

where h(u) is a homogeneous function defined on R . It is easy to see that (X1:n − d − θ)/h(X ) has ˜ the same ˜distribution as E1:n /h(E), where {Ei }ni=1 is an independent and identically distributed sample ˜ of size n from the standardized member of the scale family. Similarly, (Y1:n − d + θ)/h(Y ) has the ˜ sample same distribution as F1:n /h(F ), where {Fi }ni=1 is another independent and identically distributed ˜ of size n from the standardized member of the scale family. Making these substitutions in (6) gives a canonical form for R from which it is clear that the distribution of R is independent of unknown parameters, verifying that GPQ property (1) is satisfied. In addition, choosing Z = z it is easily seen ˜ that R is a GPQ, that r(z, z, ω) = θ, verifying that GPQ property (2) is satisfied. These observations˜ prove ˜ ˜ ˜ and moreover, that the implied confidence set reduces to an interval of the form [r1−α/2 (z), rα/2 (z)], ˜ ˜ of where the interval endpoints represent the upper- and lower-α/2 percentiles of the distribution R = r(Z , z, ω). Fig.˜3 ˜is a˜ Monte Carlo algorithm to find the required percentiles r1−α/2 (z) and rα/2 (z) using the GPQ ˜ defined by (6). Step 3 is the key part of the algorithm and is where B simulated values˜of the canonical form of r(Z , z, ω) are computed. In step 4, the lower and upper-α/2 percentiles are extracted from the ˜ ˜ to obtain the endpoints of the generalized confidence interval for θ. simulated˜values The general form of R in (6) allows some flexibility through the choice of h(u). Among ˜ function different choices we explored, the three best ones we found are the standard deviation qP n 2 i=1 (ui − u) /n, the inter-quartile range function u0.25n:n − u0.75n:n , and the Gini index of dispersion P P function ni=1 nj=1 ui − uj . We denote the GPQs formed by these choices as GPQ-SD, GPQ-IQ and GPQ-G, respectively. 3. Applications 3.1. Exponential network delays −1 −1 When network delays follow exponential distributions with means γAB and γBA in the A → B and B → A directions, respectively, then the APQ interval (5) reduces to what was derived directly as (4).

J. Li et al. / Statistical Methodology 6 (2009) 97–107

104

Fig. 3. Monte Carlo algorithm for endpoints of generalized confidence interval.

Table 1 Coverage of nominal 90% confidence intervals for θ with exponential network delays Interval

APQ APQc GPQ-SD GPQ-G GPQ-IQ avg. δ

Mean network delays, (A → B, B → A) (1, 1)

(1, 5)

Sample size, n

Sample size, n

(1, 10) Sample size, n

10

20

40

10

20

40

10

20

40

0.870 0.902 0.896 0.898 0.902 0.0743

0.891 0.907 0.900 0.906 0.898 0.0840

0.893 0.906 0.897 0.892 0.894 0.0869

0.875 0.904 0.903 0.908 0.904 0.0744

0.886 0.903 0.894 0.907 0.909 0.0840

0.891 0.906 0.897 0.896 0.897 0.0869

0.871 0.903 0.901 0.900 0.906 0.0744

0.886 0.903 0.898 0.900 0.902 0.0838

0.892 0.907 0.899 0.902 0.900 0.0870

Table 2 Expected width of nominal 90% confidence intervals for θ with exponential network delays Interval

APQc GPQ-SD GPQ-G GPQ-IQ

Mean network delays, (A → B, B → A) (1, 1)

(1, 5)

(1, 10)

Sample size, n

Sample size, n

Sample size, n

10

20

40

10

20

40

10

20

40

0.267 0.305 0.294 0.337

0.126 0.134 0.129 0.138

0.0614 0.0622 0.0607 0.0630

0.806 0.990 0.959 1.088

0.376 0.437 0.428 0.458

0.184 0.205 0.200 0.208

1.476 1.909 1.842 2.088

0.689 0.840 0.816 0.867

0.338 0.395 0.388 0.401

Identify this interval by the name APQ and denote its calibrated version by APQc. Tables 1 and 2 show the results of a simulation study to evaluate the frequentist coverage probability and the expected −1 −1 interval width for APQ and the three GPQ intervals. Following [2], choices for (γAB , γBA ) were (1, 1), (1, 5) and (1, 10) and the sample sizes considered were {10, 20, 40}. The estimated coverage probabilities (Table 1) and expected interval widths (Table 2) are based on 10,000 simulations of each interval. The last row of Table 1 shows the average of the calibrating value of δ obtained by the bootstrap calibration algorithm (Fig. 2). Table 1 illustrates that the coverage probability of the APQ interval is a little less than the nominal value, but that the calibration procedure effectively corrects it. The coverage probability of all three GPQ intervals is very close to the nominal level. Comparing expected widths, APQc is the best edging out GPQ-G. This is, perhaps, not surprising since the APQc interval is a relative profile likelihood interval under the exponential model.

J. Li et al. / Statistical Methodology 6 (2009) 97–107

105

Table 3 Coverage of nominal 90% confidence intervals for θ with Rayleigh network delays Interval

Mean network delays, (A → B, B → A)

APQ APQc GPQ-SD GPQ-G GPQ-IQ avg. δ

(1, 1)

(1, 5)

Sample size, n

Sample size, n

(1, 10) Sample size, n

10

20

40

10

20

40

10

20

40

0.754 0.921 0.901 0.906 0.904 0.0143

0.816 0.911 0.904 0.907 0.903 0.0311

0.839 0.910 0.901 0.897 0.904 0.0468

0.754 0.912 0.905 0.904 0.907 0.0139

0.802 0.902 0.902 0.906 0.906 0.0308

0.836 0.904 0.908 0.901 0.903 0.0463

0.752 0.906 0.901 0.902 0.904 0.0137

0.799 0.904 0.900 0.904 0.905 0.0307

0.833 0.904 0.904 0.903 0.900 0.0461

Table 4 Expected width of nominal 90% confidence intervals for θ with Rayleigh network delays Interval

Mean network delays, (A → B, B → A)

APQc GPQ-SD GPQ-G GPQ-IQ

(1, 1)

(1, 5)

Sample size, n

Sample size, n

(1, 10) Sample size, n

10

20

40

10

20

40

10

20

40

0.558 0.564 0.564 0.714

0.341 0.334 0.335 0.374

0.224 0.215 0.215 0.227

1.828 1.922 1.930 2.381

1.245 1.158 1.155 1.270

0.830 0.753 0.754 0.790

2.784 3.745 3.780 4.655

2.089 2.277 2.278 2.506

1.548 1.482 1.486 1.556

3.2. Rayleigh network delays When network delays follow Rayleigh distributions with scale parameters γAB and γBA , respectively, their distribution functions can be written as FAB (u) = H(γAB u) and FBA (u) = H(γBA u), respectively, where H(u) = 1 − exp(−u2 /2) is the standardized Rayleigh distribution function. As noted earlier, the Gnedenko condition holds for the Rayleigh family with τ = 2. The APQ interval ˆ α), U (X , Y ; ρ, ˆ α)], defined in conjunction with (5), requires the sequence of constants bn [L(X , Y ; ρ, ˜ ˜ of the percentile function zα (ρ) so that zα (ρ) ˆ can be evaluated. and˜ a˜characterization p The constants bn are the reciprocal solutions to 1 − exp(−u2 /2) = 1/n, and thus bn ∼ n/2. Theorem 2 implies zα (ρ) is the upper-α percentile of W (ρ) = −Z /(ξρ)I(Z ≤ 0) + (Z /ξ)I(Z

> 0),

where Z = E2 − ρE1 with {Ei }i=1 independent and identically distributed as p standardized Weibull random variables with shape parameter equal to two, ρ = γ /γ and ξ = π/2. Thus, using the BA AB qP Pn n 2/ 2 , a straightforward way to obtain z (ρ) ( Y − Y ) conditional MLE ρˆ = ( X − X ) i 1 : n i 1 : n α ˆ is to i =1 i =1 2

ˆ . use percentiles from a set of simulated values of W (ρ) For the same simulation set-up used in Section 3.1, Table 3 again shows that calibration of the APQ is both necessary and effective. Table 4 shows that in terms of expected widths, the APQc and the GPQ intervals are considerably more comparable than they were in the exponential case. While APQc retains a narrow advantage in this setting, both GPQ-SD and GPG-G are quite competitive. 3.3. Heavy-tailed network delays

Reference [6] contains references where heavy-tailed distributions have been used for network delay distributions. To investigate the performance of the proposed confidence intervals in these contexts, we considered network delay distributions that were Weibull with shape parameter equal to 1/2 and scale parameters γAB and γBA , respectively. The distribution functions √ can be written as FAB (u) = H(γAB u) and FBA (u) = H(γBA u), respectively, where H(u) = 1 − exp(− u). As noted earlier, the Gnedenko condition holds for this Weibull family with τ = 1/2. To apply Theorem 2 to find

J. Li et al. / Statistical Methodology 6 (2009) 97–107

106

Table 5 Coverage of nominal 90% confidence intervals for θ with heavy-tailed network delays Interval

APQ APQc GPQ-SD GPQ-G GPQ-IQ avg. δ

Mean network delays, (A → B, B → A) (1, 1)

(1, 5)

(1, 10)

Sample size, n

Sample size, n

Sample size, n

10

20

40

10

20

40

10

20

40

0.873 0.908 0.897 0.899 0.893 0.0284

0.883 0.904 0.901 0.903 ,901 0.0188

0.889 0.903 0.893 0.894 0.895 0.0146

0.869 0.902 0.893 0.894 0.895 0.0277

0.0.881 0.903 0.896 0.899 0.905 0.0185

0.889 0.905 0.898 0.899 0.901 0.0144

0.871 0.906 0.893 0.893 0.903 0.0273

0.884 0.906 0.901 0.901 0.896 0.0184

0.877 0.897 0.892 0.894 0.894 0.0144

Table 6 Expected width of nominal 90% confidence intervals for θ with heavy-tailed network delays Interval

APQc GPQ-SD GPQ-G GPQ-IQ

Mean network delays, (A → B, B → A) (1, 1)

(1, 5)

Sample size, n

Sample size, n

(1, 10) Sample size, n

10

20

40

10

20

40

10

20

40

0.0493 0.0658 0.0582 0.0570

0.0113 0.0137 0.0122 0.0121

0.00273 0.00299 0.00273 0.00272

0.141 0.201 0.178 0.174

0.0323 0.0430 0.0382 0.0377

0.00781 0.00947 0.00863 0.00852

0.247 0.386 0.341 0.333

0.0570 0.0793 0.0706 0.0693

0.0138 0.0176 0.0161 0.0159

the APQ interval, we obtain bn ∼ n2 , ξ = 2. Using ρˆ =

P n

i =1

(Xi − X1:n )/

p

Pn

i =1

2

(Yi − Y1:n ) , the

p

ˆ α), U (X , Y ; ρ, ˆ α)] can thus be constructed in the same way as described in the APQ interval [L(X , Y ; ρ, ˜ ˜ previous section.˜ ˜ Using the same simulation set-up as before, Tables 5 and 6 show the coverage probability and expected interval width for the APQ interval and the three GPQ intervals. From Table 5 we see the calibration of APQ is both necessary and effective and that all of the GPQ intervals achieve satisfactory coverage probability. From Table 6 we see that APQc again has the narrowest width, edging out both GPQ-G and GPQ-IQ which offer comparable widths. 4. Summary and discussion Two approaches for obtaining a confidence interval for the offset parameter θ have been described for cases where the network delay distributions are members of a one-parameter scale family. The APQ interval is a relative profile likelihood confidence interval when network delays are exponentially distributed. Use of profile likelihood techniques with other scale families can be difficult in the sense that customization of the general asymptotic theory will be necessary, and moreover, inversion of the relative profile likelihood will not necessarily produce an interval. These observations prompted the consideration of using the APQ interval more generally when network delays are modeled with distributions in other scale families. Theorem 2 provides the necessary insight into how to modify the interval to make it suitable for use in these other contexts. The GPQ intervals, applicable to all scale families, are comparatively simpler from a computational point of view — though not necessarily from a conceptual point of view. The coverage probability in a frequentist sense of generalized confidence intervals cannot be taken for granted, though there is a growing literature of evidence that they work well when viewed this way. Examination of their performance within the context of our application adds to that evidence, though it also demonstrates that the selection of the generalized pivot can have a significant impact on the expected width of the interval. Future work includes consideration of shape-scale families and seeking a robust confidence interval for θ.

J. Li et al. / Statistical Methodology 6 (2009) 97–107

107

References [1] K.-L. Noh, Q. Chaudhari, E. Serpedin, B. Suter, Power-efficient clock synchronization using two-way timing message exchanges in wireless networks, in: Proceedings of the Military Communications Conference, 2006. [2] D.R. Jeske, A. Sampath, Estimation of clock offset using bootstrap bias-correction techniques, Technometrics 45 (3) (2003) 256–261. [3] V. Paxson, On calibrating measurements of packet transit times, in: Proceedings of SIGMETRICS ’98, 1998. [4] D.R. Jeske, On the maximum likelihood estimation of clock offset, IEEE Transactions on Communications 53 (2005) 53–54. [5] D.R. Jeske, A. Sampath, A real example that illustrates interesting properties of bootstrap bias-corrected estimators, The American Statistician 57 (1) (2003) 62–65. [6] D.R. Jeske, A. Chakravartty, Effectiveness of bootstrap bias correction in the context of clock offset estimators, Technometrics 48 (2006) 530–538. [7] S. Weerahandi, Generalized confidence intervals, Journal of the American Statistical Association 88 (1993) 899–905. Correction in vol. 89, p. 726, 1994. [8] B. Gnedenko, Sur La Distribution Limite Du Terme Maximum D’Une Serie Aleatoire, Annals of Mathematical Statistics 44 (1943) 423–453. [9] N.R. Patel, R.L. Smith, The asymptotic extreme value distribution of the sample minimum of a concave function under linear constraints, Operations Research 31 (1983) 789–794. [10] W. Loh, Calibrating confidence coefficients, Journal of the American Statistical Association 82 (1987) 155–162. [11] S. Weerahandi, Exact Statistical Methods for Data Analysis, Springer-Verlag, New York NY, 1995.