{'fraction_non_alphanumeric': 0.03401133711237079, 'fraction_numerical': 0.011360929833754108, 'mean_word_length': 4.563402676560223, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Factor Graph Approach to Clock Offset Estimation in Wireless Sensor Networks

The problem of clock offset estimation in a two way timing message exchange regime is considered when the likelihood function of the observation time stamps is Gaussian, exponential or log-normally distributed. A parametrized solution to the maximum likelihood (ML) estimation of clock offset, based on convex optimization, is presented, which differs from the earlier approaches where the likelihood function is maximized graphically. In order to capture the imperfections in node oscillators, which may render a time-varying nature to the clock offset, a novel Bayesian approach to the clock offset estimation is proposed by using a factor graph representation of the posterior density. Message passing using the max-product algorithm yields a closed form expression for the Bayesian inference problem. Several lower bounds on the variance of an estimator are derived for arbitrary exponential family distributed likelihood functions which, while serving as stepping stones to benchmark the performance of the proposed clock offset estimators, can be useful in their own right in classical as well Bayesian parameter estimation theory. To corroborate the theoretical findings, extensive simulation results are discussed for classical as well as Bayesian estimators in various scenarios. It is observed that the performance of the proposed estimators is fairly close to the fundamental limits established by the lower bounds.


I. INTRODUCTION
W IRELESS sensor networks (WSNs) typically consist of a large number of geographically distributed sensor nodes, deployed to observe some phenomenon of interest. The nodes constituting such a network are low cost sensors that have limited abilities of data processing and communication. WSNs envisage tremendous applications in such diverse areas as industrial process control, battlefield surveillance, health monitoring, target localization and tracking, etc., [1]. With the recent advances in digital circuit technology, WSNs are expected to play a pivotal role in future wireless communications.
Clock synchronization in sensor networks is a critical component in data fusion and duty cycling operations, and has gained widespread interest in the past few years. Most of the current methods consider sensor networks exchanging time stamps based on the time at their respective clocks. A survey of the popular approaches employed in practice for timing synchronization is presented in [2] and [3]. The one-way message exchange mechanism involves a reference node broadcasting its timing information to other nodes in a network. The receiver nodes record the arrival of these messages with respect to their own clock. After several such time stamps have been exchanged, the nodes estimate their offsets based on these observations. A particular case of this approach is the flooding time synchronization protocol (FTSP) [4] which uses regression to estimate the clock offset. On the other hand, through a two-way timing exchange process, adjacent nodes aim to achieve pairwise synchronization by communicating their timing information with each other. After a round of N messages, each nodes tries to estimate its own clock parameters. The timing-sync protocol for sensor networks (TPSNs) [5] uses this strategy in two phases to synchronize clocks in a network. The level discovery phase involves a spanning tree based representation of a WSN while nodes attempt to synchronize with their immediate parents using a two-way message exchange process in the synchronization phase. In receiver-receiver synchronization, nodes collect time stamps sent from a common broadcasting node and utilize them to adjust their clocks. The reference broadcast synchronization (RBS) protocol [6] uses reference beacons sent from a master node to establish a common notion of time across a network. An alternative framework for network-wide distributed clock synchronization consists of recasting the problem of agreement on oscillation phases and/or frequencies as a consensus based recursive model in which only local message passing is required among nodes. By assuming a connected network, it is possible to design efficient distributed algorithms by carefully choosing the update function. Under this framework, [9] proposed a Laplacianbased algorithm for establishing agreement on oscillation frequencies all over the network based on standard consensus. A combined agreement over both clock phases and frequencies has been studied in [10], by making use of state-of-the-art fast consensus techniques. Scalable synchronization algorithms for large sensor networks are developed in [11] and [12] inspired by mathematical biology models justifying synchrony in the biological agents.
The clock synchronization problem in a WSN offers a natural statistical signal processing framework whereby, the clock parameters are to be estimated using timing information from various sensors [22]. A model based synchronization approach to arrest the clock drifts is explored in [8]. The impairments in message transmission arise from the various delays experienced by the messages as they travel through the transmission medium. Therefore, a crucial component of In this work, considering the classic two-way message exchange mechanism, a unified framework for the clock offset estimation problem is presented when the likelihood function of the observations is Gaussian, exponential or log-normally distributed. A parameterized solution is proposed for the ML estimation of clock offset by recasting the likelihood maximization as an instance of convex optimization. In order to incorporate the effect of time variations in the clock offset between sensor nodes, a Bayesian inference approach is also studied based on a factor graph representation of the posterior density. The major contributions of this work can be summarized as follows.
1) A unified framework for ML estimation of clock offset, based on convex optimization, is presented when the likelihood function of the observations is Gaussian, exponential and log-normally distributed. The proposed framework recovers the already known results for Gaussian and exponentially distributed likelihood functions and determines the ML estimate in case of log-normal distribution. Hence, the proposed convex optimization based approach represents a simpler alternative, and a more general derivation of ML estimator, which bypasses the graphical analysis used in [14] to maximize the likelihood function. 2) In order to capture the time variations in clock offsets due to imperfect oscillators, a Bayesian framework is presented by considering the clock offset as a random Gauss-Markov process. Bayesian inference is performed using factor graphs and the max-product algorithm. The message passing strategy yields a closed form solution for Gaussian, exponential and log-normally distributed likelihood functions. This extends the current literature to cases where the clock offset may not be deterministic, but is in fact a random process. 3) In order to evaluate the performance of the proposed estimators, classical as well as Bayesian bounds are derived for arbitrary exponential family distributed likelihood functions, which is a wide class and contains almost all distributions of interest. While these results aid in comparing various estimators in this work, they can be useful in their own right in classical and Bayesian estimation theory.
This paper is organized as follows. The system model is outlined in Section II. The ML estimation of clock offset based on convex optimization is proposed in Section III. The factor graph based inference algorithm for the synchronization problem in a Bayesian paradigm is detailed in Section IV and a closed form solution is obtained. Section V presents several theoretical lower bounds on the variance of an estimator evaluated in the classical as well as Bayesian regime. Simulation studies are discussed in Section VI which corroborate the earlier results. Finally, the paper is concluded in Section VII along with some directions for future research.

II. SYSTEM MODEL
The process of pairwise synchronization between two nodes S and R is illustrated in Fig. 1. At the jth message exchange, node S sends the information about its current time through a message including time stamp T 1 j . Upon receipt of this message, Node R records the reception time T 2 j according to its own time scale. The two-way timing message exchange process is completed when node R replies with a synchronization packet containing time stamps T 2 j and T 3 j which is received at time T 4 j by node S with respect to its own clock. After N such messages have been exchanged between nodes S and R, node S is equipped with time stamps The impairments in the signaling mechanism occur due to a fixed propagation delay, which accounts for the time required by the message to travel through the transmission medium, and a variable network delay, that arises due to queuing delay experienced by the messages during transmission and reception. By assuming that the respective clocks of nodes S and R are related by C R (t) = θ + C S (t), the two-way timing message exchange model at the jth instant can be represented as where d represents the propagation delay, assumed symmetric in both directions, and θ is offset of the clock at node R relative to the clock at node S. X j and Y j are the independent and identically distributed variable network delays. By the system in (1) can be equivalently expressed as By further defining the model in (2) can be written as The goal is to determine precise estimates of ξ and ψ using observations {U j , V j } N j=1 . An estimate of θ can, in turn, be obtained using (3) as follows Accurate modeling of the variable delays, X j and Y j , has been a topic of interest in recent years. Several distributions have been proposed that aim to capture the random effects caused by the queuing delays [7]. These distributions include exponential, gamma, log-normal and Weibull. In addition, the authors in [17] argued that X j and Y j result from contributions of numerous independent random processes and can, therefore, be assumed to be Gaussian. The ML estimate of d and θ for the case of exponential distribution was determined in [14]. Recently, the minimum variance unbiased estimate (MVUE) of the clock offset under an exponentially distributed network delay was proposed in [20]. In this work, instead of working with a specific distribution, a general framework of the clock synchronization problem is proposed that yields a parameterized solution of the clock offset estimation problem in the classical as well Bayesian regime when the likelihood function of the observations, U j and V j , is Gaussian, exponential or log-normally distributed.
In particular, the general notation used when the likelihood function of the observations U T is Gaussian or log-normally distributed is given below.
Unconstrained Likelihood: where η ξ (U j ) and η ψ (V j ) are sufficient statistics for estimating ξ and ψ, respectively. The log-partition functions φ ξ (.) and φ ψ (.) serve as normalization factors so that f (U; ξ) and f (V; ψ) are valid probability distributions. The likelihood function is called 'unconstrained' since its domain is independent of the parameters ξ and ψ. Similarly, the general notation used for an exponentially distributed likelihood function is given below.
Constrained Likelihood: where the indicator function I(.) is defined as and the roles of η ξ (U j ), η ψ (V j ), φ ξ (.) and φ ψ (.) are similar to (5) and (6). The likelihood function is called constrained since its domain depends on the parameters ξ and ψ. It must be noted that the likelihood functions (5)- (8) are expressed in terms of general exponential family distributions. This approach helps to keep the exposition sufficiently general and also allows us to recover the known results for the ML estimation of clock offset for Gaussian and exponentially distributed likelihood functions [14] [15], and determine the ML estimator of the clock offset in case of log-normally distributed likelihood function, as shown in Section III. The proposed approach will also prove useful in investigating a unified novel framework for clock offset estimation in the Bayesian setting for Gaussian, exponential or log-normally distributed likelihood functions, as will be shown in Section IV.
Some key ingredients of the proposed solution for the clock offset estimation problem, based on the properties of exponential family, can be summarized as follows [25].
1) The mean and variance of the sufficient statistic η ξ (U j ) are expressed as 2) The moment generating function (MGF) of the statistic η ξ (U j ) is given by 3) The non-negativity of the variance σ 2 η ξ in (10) implies that the log-partition function φ ξ (.) is convex. 4) For Gaussian, exponential and log-normally distributed likelihood functions, the log-partition function φ ξ (ξ) can be expressed as a second degree polynomial given by The coefficient a ξ in this approximation can be obtained using the variance of the statistic η ξ (U j ), which is assumed known. Using (10), a ξ is given by If the statistical moment in (10) is not available, the empirical moment can be substituted since it readily follows from the weak law of large numbers that Similar expressions can also be written for η ψ (V j ), M η ψ (h) and φ ψ (ψ), respectively. Using the aforementioned properties of the exponential family, the ML as well as Bayesian estimates of ξ and ψ are to be determined utilizing the data set {U j , V j } N j=1 , based on (12).

III. MAXIMUM LIKELIHOOD ESTIMATION
In this section, the ML estimates of θ are obtained by recasting the likelihood maximization as an instance of convex optimization. This approach differs from the graphical arguments used to maximize the likelihood in [14]. The specific cases of unconstrained and constrained likelihood functions are considered separately.

A. Unconstrained Likelihood
Using (5), (6) and (12), the unconstrained likelihood functions are given by The ML estimates of ξ and ψ can now be expressed aŝ Theorem 1: The likelihood maximization problems (15) and (16) are strictly concave and the ML estimates are given byξ (17b) Hence, the ML estimatorθ ML for the clock offset can be written asθ Proof: The ML estimate of ξ in (15) can be equivalently determined by maximizing the exponent in the likelihood function. It can be easily verified that the exponent in (15), which is a quadratic function of ξ, is strictly concave [29]. A similar explanation applies to the ML estimate of ψ. The ML estimates in (17) can be obtained by setting the first derivative of the exponent with respect to ξ (resp. ψ) to zero. The estimateθ ML in (18) can be obtained by invoking the invariance principle [28]. Hence, the proof readily follows.
1) Gaussian Distributed Likelihood Function: A particular application of Theorem 1 is the case when the likelihood which can be rearranged as By comparing with (13), we have and the ML estimate using (17a) is given bŷ By a similar reasoning, the ML estimate for ψ (17b) is given byψ Using (18), the ML estimate for the offset θ can be expressed asθ The above estimate coincides exactly with the one reported in [15].
2) Log-Normally Distributed Likelihood Function: Another application of Theorem 1 is when the samples U j and V j are log-normally distributed. In this case, we have A comparison with (13) yields Following a similar line of reasoning, The ML estimator for θ can obtained from (18) using (17a) and (17b), and is given bŷ

B. Constrained Likelihood
Using (7), (8) and (12), the constrained likelihood functions are given by The resulting ML estimates of ξ and ψ can be obtained aŝ Theorem 2: The likelihood maximization problems (27) and (28) are strictly concave and the ML estimates can be expressed asξ where U (1) and V (1) denote the first order statistics of the samples U j and V j , respectively. The ML estimatorθ ML for the clock offset is given bŷ Proof: By using arguments similar to Theorem 1 and noting that the N constrains U j ≥ ξ (resp. V j ≥ ψ) are linear functions of ξ (resp. ψ), the proof of concavity readily follows. Also, the likelihood maximization problems (27) and (28) can be equivalently expressed aŝ The unconstrained maximizerξ of the objective function in (27) is given by (17a). Ifξ ≤ U (1) , thenξ ML =ξ. On the other hand, if U (1) ≤ξ, then the ML estimate is given byξ ML = U (1) using concavity of the objective function. Combining the two cases, the ML estimate in (29a) is obtained. A similar explanation applies to the ML estimateψ ML in (29b) andθ ML follows from the invariance principle.

1) Exponentially Distributed Likelihood Function:
For the case when the likelihood functions are exponentially distributed, the density function of the samples U j can be written as [14] f where λ −1 ξ is the mean of the delays X j . The density function can be rearranged as Comparing the above formulation with (25), Using (29a), the ML estimate is given bŷ Employing a similar reasoning, Using (30), the ML estimate of θ is given bŷ which coincides exactly with the one reported in [14], where it is derived using graphical arguments.
Remark 1: The ML estimation method outlined above differs from the previous work [14] in that it is based on convex optimization, while [14] maximized the likelihood graphically. Hence, this approach presents an alternative view of the ML estimation of the clock offset. It also allows us the determine the ML estimator of θ when the likelihood function is lognormally distributed. In addition, Theorem 2 will also be useful in Section IV where estimation of θ in the Bayesian regime is discussed .

IV. A FACTOR GRAPH APPROACH
The imperfections introduced by environmental conditions in the quartz oscillator in sensor nodes results in a time-varying clock offset between nodes in a WSN. To cater for such a temporal variation, in this section, a Bayesian approach to the clock synchronization problem is adopted by representing the a-posteriori density as a factor graph. The inference is performed on the resulting factor graph by message passing using max-product algorithm. To ensure completeness, a brief description of factor graphs and the max-product algorithm is provided below.
A factor graph is a bipartite graph that represents a factorization of a global function as a product of local functions called factors, each factor being dependent on a subset of variables. Factor graphs are often used to produce a graphical model depicting the various inter-dependencies between a collection of interacting variables. Each factor is represented by a factor node and each variable has an edge or a half-edge. An edge connects a particular variable to a factor node if and only if it is an argument of the factor expressed by the factor node [24].
Inference can be performed by passing messages (sometimes called beliefs) along the edges of a factor graph. In particular, max-product algorithm is used to compute the messages exchanged between variables and factor nodes. These messages can be summarized as follows variable to factor node : factor node to variable : where Z = n (f ) is the set of arguments of the local function f . The marginal distributions associated with each variable can be obtained by the product of all incoming messages on the variable. In order to sufficiently capture the temporal variations, the parameters ξ and ψ are assumed to evolve through a Gauss-Markov process given by where w k and v k are i.i.d such that w k , v k ∼ N (0, σ 2 ). The posterior pdf can be expressed as where uniform priors f (ξ 0 ) and f (ψ 0 ) are assumed. Define where the likelihood functions are given by based on (12). The resulting factor graph representation of the posterior pdf is shown in Fig. 2.

Remark 2:
A few important observations of this representation can be summarized below.
• Notice that the substitution in (3) renders a cycle-free nature to the factor graph. Therefore, inference by message passing on such a factor graph is indeed optimal [24]. • The two factor graphs shown in Fig. 2 have a similar structure and hence, message computations will only be shown for the estimateξ N . Clearly, similar expressions will apply toψ N . In addition, only the case of constrained likelihood will be considered, since the case of an unconstrained likelihood is subsumed, as will be shown shortly. The clock offset estimator θ N can be obtained fromξ N andψ N using (4).
, the constrained likelihood function for the samples U k can be written as The message passing strategy starts by sending a message from the factor node f N to the variable ξ N . The variable ξ N relays this message to the factor node δ N N −1 . The factor node computes the product of this message with the factor δ N N −1 and sends the resulting message to the variable ξ N −1 after 'summarizing' over the variable ξ N . In the max-product algorithm, a 'max' function is used as a summary propagation operator (cf. (36)). These messages are computed as which can be rearranged as Letξ N be the unconstrained maximizer of the exponent in the objective function above, i.e., This implies thatξ Following a line of reasoning similar to Theorem 2, it follows thatξ However,ξ N depends on ξ N −1 , which is undetermined at this stage. Hence, we need to further traverse the chain backwards. Assuming thatξ N ≤ U N ,ξ N from (42) can be plugged back in (40) which after some simplification yields The message passed from the variable ξ N −1 to the factor node δ N −1 N −2 is the product of the message (43) and the message received from the factor node f N −1 , i.e., Upon receipt of this message, the factor node δ N −1 N −2 delivers a product of this message and the factor δ N −1 N −2 to the variable node ξ N −2 after summarizing over ξ N −1 . This message can be expressed as After some algebraic steps, the message above can be compactly represented as Proceeding as before, the unconstrained maximizerξ N −1 of the objective function above is given bȳ and the solution to the maximization problem (44) is expressed asξ Again,ξ N −1 depends on ξ N −2 and therefore, the solution demands another traversal backwards on the factor graph representation in Fig. 2. By pluggingξ N −1 back in (44), it follows that which has a form similar to (43). It is clear that one can keep traversing back in the graph yielding messages similar to (43) and (45). In general, for i = 1, . . . , N − 1 we can write Using (47) and (48) with i = N − 1, it follows that Similarly, by observing the form of (43) and (45), it follows that (51) The estimateξ 0 can be obtained by maximizing the received message in (51). It can be noticed from the structure of the factor graph that this maximization is inherently unconstrained i.e.,ξ The estimate in (52) can now be substituted in (49) to yield ξ 1 , which can then be used to solve forξ 1 in (50). Clearly, this chain of calculations can be continued using recursions (47) and (48). Define A key property of the function g ξ,k (.), which proves useful in the quest for a closed form solution, can be summarized in the following lemma. Lemma 1: For real numbers a and b, the function g ξ,k (.) defined in (53) satisfies g ξ,k (min(a, b)) = min (g ξ,k (a), g ξ,k (b)) .
The estimateξ 3 can, therefore, be compactly represented aŝ Hence, one can keep estimatingξ k at each stage using this strategy. Note that the estimator only depends on functions of data and can be readily evaluated.
Generalizing this framework, the closed form expression for the clock offset estimateθ N is given by the following theorem.
Theorem 3: The state estimatesξ N andψ N for the posterior pdf in (37) can be expressed aŝ ψ 0 and the factor graph based clock offset estimate (FGE)θ N is given byθ Proof: The proof follows from the discussion above and using (4).

Remark 3:
The closed form expressions in Theorem 3 enable the estimation of the clock offset, when it may be timevarying and the likelihood functions, f (U k |ξ k ) and f (V k |ψ k ), have a Gaussian, exponential or log-normal distribution.

A. Gaussian Distributed Likelihood Function
A particular case of the Bayesian framework described above occurs when likelihood functions f (U k |ξ k ) and The aforementioned Gaussian distribution constitutes an unconstrained likelihood function, i.e., the domain of the pdf is independent of the unknown parameter ξ k . It is clear from the message passing approach that at each stage k of the factor graph, the unconstrained maximizerξ k is the actual solution to the likelihood maximization problem i.e.,ξ k =ξ k ∀k = 1, . . . , N . Hence, the unconstrained likelihood maximization problem is subsumed in the message passing framework for constrained likelihood maximization. It follows from Theorem 3 thatξ N for Gaussian distributed observations U k in (57) is given bŷ ξ 0 whereξ 0 and G N ξ,1 (.) are defined in (52) and (55), respectively. Evaluatingξ N requires to determine the constants in (41) and (46). By comparing (57) with (39), we have Using these values for α ξ,k and β ξ,k , (41) and (46) can be written as for i = 1, . . . , N −1. Using similar arguments, it can be shown that the estimateψ N is given bŷ It follows from (4) that the FGE,θ N , can be expressed aŝ The behavior ofθ N can be further investigated for the case when the noise variance σ 2 in the Gauss-Markov model goes to zero. Consider where the constants A ξ,N , B ξ,N , C ξ,N and D ξ,N are given by (59). After some algebraic steps, we have g ξ,N (ξ) = σ 2 ξ ξ + σ 2 U N σ 2 ξ + σ 2 . As σ 2 → 0, g ξ,N (ξ) → ξ. Similarly, it can be shown that g ξ,N −1 (ξ) → ξ as σ 2 → 0. Hence, it follows that in the low system noise regime, as σ 2 → 0

Similarly, it can be shown that
which can be proven to be equal to the ML estimator (22).

B. Log-Normally Distributed Likelihood Function
The log-normally distributed likelihood function in the Bayesian regime can be expressed as Clearly, the only difference here with the Gaussian distribution is a redefinition of β ξ,k . The expression ofξ N in this case is againξ where G N ξ,1 (.) andξ 0 are given by (55) and (52), respectively. The recursively evaluated constants in (41) and (46) can be written as for i = 1, . . . , N − 1. Similar arguments apply to ψ. Hence, the FGEθ N can be expressed aŝ Again, as the Gauss-Markov system noise σ 2 → 0, the above estimator approaches its ML counterpart (24).

C. Exponential Distribution
Theorem 3 can also be used to derive a Bayesian estimator ξ N for the exponentially distributed likelihood case considered in [14]. In this case, we have where λ −1 ξ is the mean network delay of X k . A comparison of (63) with (39) reveals that For these values of α ξ,k and β ξ,k , the constants A ξ,k , B ξ,k , C ξ,k and D ξ,k are given by Similarly it can be shown that and so on. The estimatorξ 0 at the last step can be evaluated asξ Using (64) and Theorem 3, it readily follows that Similarly, for the pdf [14] f the estimateψ N is given bŷ and the estimateθ N can be obtained using (56), (65) and (66) asθ As the Gauss-Markov system noise σ 2 → 0, (67) yieldŝ which is the ML estimator given by (34).

V. CLASSICAL AND BAYESIAN BOUNDS
To evaluate the performance of the estimators derived in the preceding sections, classical as well as Bayesian lower bounds on the variance of the estimators are discussed. The placement of a lower bound allows one to compare estimators by plotting their performance against the bound. It must be emphasized here that the results in this section assume no specific form of the log-partition function and are therefore, valid for arbitrary distributions from the exponential family, which is a wide class and contains almost all distributions of interest. Hence, these results are fairly general and can be useful in their own right in classical as well as Bayesian parameter estimation theory, and at the same time will be used as a stepping stone towards comparing the estimators developed thus far.
The likelihood function of the data is considered an arbitrary member of the exponential family of distributions. In addition, depending on whether the domain of the likelihood depends on the parameter to be estimated, both cases of unconstrained as well as constrained likelihood functions are discussed to maintain full generality. The general expressions for the unconstrained and constrained likelihood functions for observations Z given by Unconstrained Likelihood: Constrained Likelihood: where ρ is the scalar parameter to be estimated. The goal is to derive lower bounds on the variance of estimators of ρ. For the case of classical estimation, the Cramer-Rao and Chapman-Robbins bounds are considered, while the Bayesian Cramer-Rao bound and a Bayesian version of the Chapman-Robbins bound are derived for the Bayesian paradigm.

A. Cramer-Rao Lower Bound
The Cramer-Rao lower bound (CRB) is a lower bound on the variance of an unbiased estimator of a deterministic parameter. It is useful primarily because it is relatively simple to compute. However, it relies on certain 'regularity conditions' which are not satisfied by constrained likelihood functions when the domain of the likelihood depends on the unknown parameter (cf. (69)). Hence, CRB is determined for the case of unconstrained likelihood functions only.
In particular, CRB states that the variance of an unbiased estimator of ρ is lower bounded by Theorem 4: The CRB for ρ in the unconstrained likelihood function in (68) is given by where Proof: The Fisher information for the likelihood function is given by and the proof readily follows.

B. Chapman-Robbins Bound
Chapman-Robbins bound (CHRB), proposed in [26], sets a lower bound on the variance of an estimator of a deterministic parameter. The CHRB does not make any assumptions on the differentiability of the likelihood function and the regularity conditions that often constrain the use of CRB, and is substantially tighter than the CRB in many situations. Hence, CHRB is employed to determine a lower bound on the variance of an unbiased estimator of ρ for constrained likelihood functions.
In general for a parameter ρ, the CHRB is given by which can be evaluated as shown below. Theorem 5: The CHRB for the parameter ρ given the likelihood function (69) can be expressed as where M η (h) is the MGF of the statistic η(Z j ) and with the expectation taken with respect to any Z j . Proof: The details of the proof are relegated to Appendix A.

C. Bayesian Cramer-Rao Lower Bound
The Bayesian Cramer-Rao bound (BCRB) is a lower bound on the variance of an unbiased estimator when the parameter assumes a prior density. It requires the same regularity conditions to be satisfied as its classical counterpart.
For an estimatorρ k of ρ k , the BCRB states that the variance of the estimator is bounded below by the lower-right submatrix of the inverse of the Bayesian information matrix, J −1 CR (k) [30], i.e., where the Bayesian information matrix is given by where the expectation is taken with respect to the joint pdf and It is assumed that the parameter ρ k evolves through a Gauss-Markov model given by A recursive formula to evaluate the Bayesian sub-matrix, derived in [27], is given by and the expectation is again with respect to the joint pdf.
Theorem 6: For the Bayesian framework in (76) and (77), the recursive Bayesian information matrix in (78) is given by with J CR (0) = 0. Proof: For the density functions, f (ρ k |ρ k−1 ) and f (Z k |ρ k ) in (76), it can be verified that The proof follows by plugging these quantities in (78).

D. Bayesian Chapman-Robbins Bound
A Bayesian version of the Chapman-Robbins bound (BCHRB) can be used to provide a lower bound on the variance of an estimator of ρ k when there are no regularity assumptions on the likelihood. In fact, unlike the BCRB, the BCHRB can be evaluated for constrained likelihood functions where the domain of the likelihood is dependent on the unknown parameter.
BCHRB states that the variance of an estimatorρ k of ρ k is lower bounded as Theorem 7: The BCHRB for the parameter ρ k can be expressed as Proof: See Appendix B for details.
Remark 4: Since the performance bounds are derived for arbitrary exponential family distributions, they can also prove useful in a broad sense in classical as well as Bayesian parameter estimation.

E. Relation to Clock Offset Estimation
The performance bounds derived in the preceding subsections can be used to lower bound the mean square error (MSE) of the clock offset θ. Notice the similarity between the unconstrained and constrained likelihood functions for ξ and ψ in (5)-(8) and the general exponential family likelihood function considered in (68) and (69). Therefore, the bounds derived above for ρ are also applicable to the parameter ξ (and also ψ). The MSE of the clock offset θ can, in turn, be lower bounded using the bounds on ξ and ψ.
Using (4), the following result is immediate. Proposition 1: The MSE of any estimator of θ can be expressed as where b ξ and b ψ are the biases of the estimatorsξ andψ, respectively. The explicit lower bounds on the MSE of any estimator θ for classical as well as Bayesian framework in case of Gaussian and exponentially distributed likelihood functions can be evaluated as shown below. 1) Gaussian Distribution -CRB: If the likelihood function for ξ is Gaussian distributed (19), then using (20) and (71), it is straightforward to see that the CRB for any unbiased estimatorξ is given by and a similar expression is applicable toψ as well. Using Proposition 1, it can be concluded that As a remark, it is evident in this case thatθ ML (22) is efficient in the sense that its MSE achieves (81) with equality (cf. Appendix C-A).

and (74) becomes
so that the statement of the CHRB (73) can be rewritten as and similarly forψ. Using Proposition 1, it follows that 3) Gaussian Distribution -BCRB: In the Bayesian regime, if the likelihood function for ξ is Gaussian distributed (57), by using (58) and (79), it can be seen that with J CR,ξ (0) = 0. A similar line of reasoning can be followed to derive an analogous recursion for J CR,ψ (k). The MSE of θ can be now be lower bounded as 4) Exponential Distribution -BCHRB: If the likelihood for ξ k is exponentially distributed (63), (80) turns out to be In fact, we just have to notice that φ ξ (ξ k ) is a constant function over ξ k and η ξ (U j ) = λ ξ , so that (85) becomes

VI. SIMULATION RESULTS
This section aims to corroborate the theoretical results in preceding sections by conducting simulation studies in various scenarios. The performance of both the classical as well as Bayesian estimators is to be investigated. The measure of fidelity used to rate this performance is the MSE of the estimators for θ and θ N . The parameter choice is σ ξ = σ ψ = 0.1 for both Gaussian and log-normally distributed likelihoods, while λ ξ = λ ψ = 10 for exponentially distributed likelihood functions.

A. Log-normal Distribution
No solution is reported thus far in literature in case the likelihood is log-normally distributed. The proposed estimators, θ ML (24) andθ N (62), can be used to determine an estimate for the clock offset in classical and Bayesian framework, respectively, as shown below.

1) Classical Estimation Framework:
The existing approaches in literature only consider the Gaussian and the exponential cases, therefore (24) is a new result in the stateof-the-art about clock offset estimation. Fig. 3 shows a comparison between the proposed ML estimator (MLE) (24) in case of a log-normally distributed likelihood (23) with MLEs which (wrongly) assume that the likelihood is Gaussian and exponentially distributed, respectively. The plot shows that the latter approaches are not robust with respect to the likelihood distribution, and their performance is extremely poor if their assumptions do not hold. In addition, Fig. 3 also shows that the proposed MLE (24) is efficient since it attains the CRB (as well as the CHRB).
2) Bayesian Estimation Framework: Fig. 4 plots the MSE performance of the FGE (62) as well as the BCRB and BCHRB when the likelihoods are log-normally distributed (61), and σ = 10 −4 . Firstly, it can be seen that the MSE of the proposed FGE coincides with the estimation bounds. Secondly, as in the classical estimation case, if we were to (wrongly) assume a Gaussian or exponential distribution for the likelihoods (38), the resulting FGEs would perform poorly, a fact is evident in Fig. 4 by observing the unboundedness and unpredictability of the dashed curve (Gaussian assumption for the likelihoods) and the dotted curve (likelihoods assumed exponentially distributed). This clearly establishes that the FGE (62), obtained assuming log-normally distributed likelihoods, allows a strong performance improvement with respect to existing estimators if the likelihood functions (38) are actually log-normally distributed.

B. Estimator Performance vs Estimation Bounds
It will also be useful to asses the performance of the MLEs in Gaussian and exponential cases derived in Section III against the various benchmark estimation bounds derived in Section V. Similarly, the FGEs for Gaussian and exponential distributions, proposed in Section IV, can also be compared with the Bayesian bounds to study their MSE performance.
1) Classical Estimation Framework: Fig. 5 shows the performance comparison between the MSE of the MLEs (22) and (34) for Gaussian and exponentially distributed likelihood functions against the CRB and the CHRB, respectively. Firstly, it is evident that in the case of Gaussian distribution, the CRB and the CHRB coincide. Moreover, the MSE ofθ ML also coincides with the aforementioned bounds. On the other hand, for an exponentially distributed likelihood function, due to its lack of regularity, the CRB cannot be derived, thus only the CHRB is shown. It can be observed that the MSE ofθ ML is fairly close to CHRB, even though it does not coincide with it. From Fig. 5 the MSE of the MLEs for the Gaussian and exponential distribution case can be also compared. In order to ensure a fair comparison, parameters are chosen in a way to have the same variance of the observations for both distributions. From the MSE curves, one can infer that the MSE in case of an exponentially distributed likelihood is lower than the one for a Gaussian distribution as the number of observations N increases. This behavior is expected since it can be verified by the MSE expressions (86) and (87), reported in Appendix C, that in case of a Gaussian distribution, the MSE decreases proportionally to 1/N , while in the exponential distribution case it decreases proportionally to 1/N 2 .
2) Bayesian Estimation Framework: In Fig. 6, the MSE performance of the FGEsθ N (60) and (67) is compared with BCRB and BCHRB for σ = 10 −4 . As in the classical estimation scenario, it is evident that for Gaussian distributed likelihoods, the MSE using (60) forθ N coincides with the reported bounds. The MSE of the FGE derived assuming    exponentially distributed likelihoods (67) is plotted against the BCHRB as well in Fig. 6. It is clear that the MSE is quite close to BCHRB, although not coinciding with it, as exactly was the case in the classical estimation framework.

C. Comparing Classical and Bayesian frameworks
The estimators proposed in the classical and the Bayesian framework can also be compared with each other based on their MSE performance as the system noise decreases. The aim here is to show that the latter approaches the former as σ → 0. Fig. 7 depicts the MSE for the cases of Gaussian, exponential and log-normal distribution for the likelihoods with N = 25. In the plot, the horizontal lines represent the MSEs in the classical framework, obtained with the MLEs, as shown in (86) and (87) in Appendix C. It can be observed that, for all the three considered distributions, the MSE obtained by using the FGE for estimating θ approaches the MSE of the MLEs as σ < 10 −3 .

VII. CONCLUSIONS AND FUTURE WORK
The clock synchronization problem in sensor networks has received keen attention recently owing to its central role in critical network operations as duty cycling, data fusion, and node localization. Based on a two-way timing message exchange scenario, this work proposes a unified framework for the clock offset estimation problem when the likelihood function of the observation time stamps is Gaussian, exponential and log-normally distributed. A convex optimization based ML estimation approach is presented for clock offsets. The results known thus far for Gaussian and exponentially distributed network delays are subsumed in the general approach while the ML estimator is derived when the likelihood function is lognormally distributed. In order to study the case of a possibly time-varying clock offset, a Bayesian approach is also studied using factor graphs. The novel message passing strategy results in a closed form solution of the time-varying clock offset estimation problem. In order to compare various estimators, several lower bounds on the variance of an estimator have been derived in the classical as well as the Bayesian regime for likelihood functions which are arbitrary members of the exponential family, a wide class containing several distributions of interest. The theoretical findings are corroborated by simulation studies conducted in various scenarios.
In future, it will be useful to incorporate the effect of clock skew in the clock offset estimation model. This can result in further reduction of the re-synchronization periods. In addition, the results about pairwise synchronization can be used to build a framework for network-wide synchronization across a sensor network. With ζ(.) defined in the theorem, the proof is complete.

APPENDIX B PROOF OF THEOREM 7
We have