{'fraction_non_alphanumeric': 0.043311382878645346, 'fraction_numerical': 0.020169332079021638, 'mean_word_length': 4.496587383660807, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Moving-Resting Process with Measurement Error in Animal Movement Modeling

Statistical modeling of animal movement is of critical importance. The continuous trajectory of an animal's movements is only observed at discrete, often irregularly spaced time points. Most existing models cannot handle the unequal sampling interval naturally and/or do not allow inactivity periods such as resting or sleeping. The recently proposed moving-resting (MR) model is a Brownian motion governed by a telegraph process, which allows periods of inactivity in one state of the telegraph process. The MR model shows promise in modeling the movements of predators with long inactive periods such as many felids, but the lack of accommodation of measurement errors seriously prohibits its application in practice. Here we incorporate measurement errors in the MR model and derive basic properties of the model. Inferences are based on a composite likelihood using the Markov property of the chain composed by every other observed increments. The performance of the method is validated in finite sample simulation studies. Application to the movement data of a mountain lion in Wyoming illustrates the utility of the method.


Introduction
Statistical modeling of animal movement is of great importance in addressing fundamental questions about space use, movement, resource selection, and behavior in animal ecology (Hooten et al., 2017). The explosion of telemetric data on animal movement from the recent advancements in tracking and observation technologies presents countless opportunities and challenges (Cagnacci et al., 2010;Patterson et al., 2017). Telemetry devices, like Global Positioning System (GPS) receivers, can only determine an animal's position at discrete moments in time so a continuous trajectory is never available. Programming the device to record at a very fast fix rate could approximate a continuous trajectory, but this is seldom done due to battery life limitations: it is very expensive to capture and collar an animal, so a long episodic time record is usually preferable to a highly detailed record, especially for animals that spend a great deal of time not moving around. GPS receivers often produce fixes at irregularly spaced time points, even if researchers program for regular intervals, due to environmental factors such as satellite communication issues or sky occlusion. As a result, discrete-time models such as the state space model (e.g., Jonsen et al., 2005;Patterson et al., 2008;Brett et al., 2012) are not realistic. Continuous-time models based on stochastic differential equations (SDE) (e.g., Preisler et al., 2004;Horne et al., 2007;Brillinger, 2010) can handle the irregular spacing naturally. Nonetheless, most existing work assumes perpetual motion and cannot accommodate periods of inactivity. On the time scale of most telemetry data, most animals alternate between periods of movement (foraging) and periods of inactivity (e.g., prey handling or rest Mashanova et al., 2010;Ueno et al., 2012;Jeschke, 2007). Realistic continuous-time models that accommodate inactive periods are needed.
The recently proposed moving-resting (MR) process by  is a promising model to accommodate inactive periods. The MR process is a Brownian motion governed by a telegraph or on-off process (e.g., Zacks, 2004). Specifically, it allows an animal to alternate between a moving state, during which it moves in a Brownian motion (BM), and a resting state, during which it remains motionless. The switch between the two states is governed by a telegraph process, where the holding time (or duration) of each state is assumed to follow an exponential distribution. The memoryless holding time makes the underlying state process a continuous time Markov Chain. As a consequence, the MR process can be analyzed with the help of hidden Markov model (HMM) tools (Cappé et al., 2005). The MR process is a first step towards more realistic animal movement modeling with discretely observed telemetry data where the trajectories contain evident motionless segments. Implementation of likelihood based inferences for the MR process based on dynamic programming (Pozdnyakov et al., 2019) is publicly available in an R package smam (Hu et al., 2020b).
The MR model does not accommodate the measurement error or noise of telemetric devices, which is a major limitation in applying it to animal movement data. Adding measurement error to a Brownian motion model is not crucial as long as the noise is small in comparison to the total standard deviation of the increments of the Brownian motion between two consecutive time points . In such cases, discarding the noise would not produce significant bias. The impact of the noise on inferences about MR processes, however, is much greater. For a given sequence of hidden states, the likelihood is a product of both densities and probabilities. With perfect instrumentation, if a sequence of observed locations are exactly the same, that is, there is no change in either the easting nor northing coordinates for a period of time resulting in a "flat" piece of trajectory, then the animal is known to be motionless over the time period spanned by the sequence. The likelihood contribution is the probability of staying in the motionless state instead of the density of the increment at zero. Adding even a tiny bit of noise would remove those flat pieces and, hence, cause drastic bias in the likelihood estimator of the parameters. One possible remedy is to round the observed coordinates, which enforces flat pieces. The number of such pieces, however, depends greatly on the rounding level, and there are no obvious rules to aid researchers in choosing best levels. A detailed illustration of the issue is given in Section 2.
Dealing with added noise in an MR process is challenging because it invalidates the Markov property of the joint location-state process. The transition density from one time point to the next can in principle be obtained from convolving the results for the MR process  with normally distributed noises, although computationally very intensive. A lack of Markov property of the joint location-state process means that the likelihood cannot be easily formed by multiplying these transition densities. Because the measurement errors are continuous, the dynamic programming tools of HMM based on a finite number of hidden states (Cappé et al., 2005) are not applicable. The generic simulation based inferences such as iterated filtering (Ionides et al., 2011(Ionides et al., , 2015 or particle Markov chain Monte Carlo (Andrieu et al., 2010), available in R package pomp (King et al., 2016), are not applicable to our investigation due to the complexity of the MR process with measurement error.
Our contribution in this paper is a toolbox for applying the MR process with measurement error to animal movement modeling. First, we show that discarding the measurement error, even tiny ones, causes severe bias in estimation, and that rounding does not provide any satisfactory solution. To make inferences for MR process with measurement error, we establish that, after thinning every other observation, the remaining observations are location-state Markov. This facilitates a composite likelihood which contains two true likelihood components, one based on odd-numbered observations and the other based on even-numbered observations. The true likelihood of each component is computed with dynamic programming. The variance of the maximum composite likelihood estimator can be estimated through parametric bootstrap. The validity of the approach is confirmed in a simulation study. We then apply the approach to model the movement data of a mountain lion in Wyoming, whose trajectory is known to have long inactive periods. Our methods are publicly available in an R package smam (Hu et al., 2020b) with efficient C++ code.

Moving-Resting Process
The MR process is a Brownian motion with an infinitesimal variance that is governed by an alternating renewal process with two different holding times. Let random variables {M i } i≥1 be independent exponential variables with rate λ 1 , and {R i } i≥1 be independent exponential variables with rate λ 0 . These are the holding times. There are two possible alternating sequences of the holding times, (M 1 , R 1 , M 2 , R 2 , . . . ) or (R 1 , M 1 , R 2 , M 2 , . . . ). Which one represents a particular realization depends on an initial distribution. A continuous time state process, S(t), t ≥ 0, takes only two values, 0 and 1, and it is defined by the holding times. In particular, for sequence (M 1 , R 1 , M 2 , R 2 , . . . ), if there exists k ≥ 0 such that then S(t) = 0, otherwise, S(t) = 1. It is well-known that the state process is stationary, if the initial probability of {S(0) = 1} is set as and the initial probability of {S(0) = 0} is set as p 0 = 1 − p 1 .
An MR process X(t), t ≥ 0, is defined by a stochastic differential equation where B(t) is the standard Brownian motion, and σ is a volatility parameter. It is important to note that {X(t)} t≥0 itself is not Markov, but the location-state process {X(t), S(t)} is a Markov process with stationary increments.
Properties and inferences of the MR process have been studied in  and Pozdnyakov et al. (2019). A key element is the distribution of occupation times, that is, the total time spent in the moving state by time t M (t) = t 0 S(s)ds, and the total time spent in the resting state R(t) = t − M (t). Let P i (·) be the conditional probability Pr(·|S(0) = i). Zacks (2004) derived computationally efficient formulas for the following (defective) densities for 0 < w < t: p 11 (w, t)dw = P 1 (M (t) ∈ dw, S(t) = 1), Having this at hand, one can derive the marginal distribution of the increment X(t)−X(0).
Without loss of generality, let X(0) to be 0, and X(t) becomes the increment from time 0 to time t. Then, the joint distribution of the increment X(t) and S(t), t > 0, is where δ 0 (x) is the delta function with an atom at 0, x ∈ R, and h ij (x, t), i, j ∈ {0, 1}, are functions derived in : with φ(·; σ 2 ) being the density function of normal distribution N(0, σ 2 ). The marginal dis- For actual animal movement data, we never observe the exact values of X(t) but only X(t) with added measurement errors. For an MR process, the probability of observing a zero increment is strictly positive. Adding noise makes this probability zero. Rounding can help, but it is not trivial to come up with an appropriate rounding level. Figure 1 (upper left) shows the easting/northing coordinates in the Universal Transverse Mercator (UTM) coordinate system of a female mountain lion in 2012 in the Gros Ventre mountain range, Wyoming. The patterns of resting -places where both lines are more or less flatand moving are readily apparent, which can hardly be captured by any existing model that assumes perpetual movements. The other three panels of Figure 1 show the coordinates of a simulated MR path without noise and with noise of two levels. The pattern is very similar to that in the upper left panel for the female mountain lion. The similarity is obvious, suggesting that an MR process might be a good model, but as shown next, ignoring the noise is disastrous in estimating the model parameters.
We demonstrate the impact of noise on estimation by a simulation study. Consider an MR process with parameters λ 1 = 1 hour −1 , λ 0 = 0.5 hour −1 , σ = 1 km/hour 1/2 . The measurement errors were independent Gaussian noise with standard deviation 0.05 km (1/20 of σ) and 0.01 km (1/100 of σ). The time intervals between consecutive observations was 5 hours. We generated 100 datasets, each with sampling horizon 1000 hours. The maximum likelihood estimates based on the MR process were obtained for dataset with and without noise, where for data with noise, three levels of rounding were considered, 10, 50, and 100 meters. Table 1 summarizes the parameter estimates based on the 100 replicates. When there was no noise, the point estimates were good, recovering the true parameters with high accuracy. For data with noise but no rounding, the optimization did not converge for most replications because of the Nelder-Mead simplex degeneracy (Nelder and Mead, 1965). With the help of various levels of rounding, the convergence percentage increases as the rounding becomes coarser, and the bias decreases but remain notable. This is true even for the cases with a noise standard deviation 0.01 km. It is indeed unclear how to choose an appropriate rounding level. A practical model should handle the measurement errors directly.

Moving-Resting Process with Measurement Error
Suppose the observations are recorded at times t 0 = 0, t 1 , . . . , t n . Let { k } k=0,...,n be independent and identically normally distributed random variables with mean 0 and variance σ 2 . An MR process with measurement error (MRME) Z(t k ), k = 0, 1, . . . , n, is the superimposition of a measurement error and the exact location: where X(·) is an MR process in Equation (1), and k 's are independent N(0, σ 2 ) noises.
Some properties of the process {Z(t k )} k=0,...,n are in order. Obviously, it is not Markov.
Neither is the location-state process {Z(t k ), S(t k )} k=0,...,n . To get a Markov process one might consider including the measurement errors. It is true that {Z(t k ), S(t k ), k } k=0,...,n is Markov. But the cardinality of hidden states (S(t k ), k ) is a continuum in this case. This makes the dynamic programming approach for computing likelihood infeasible. To address this difficulty we suggest considering the process The process is Markov, because the increment of the Brownian motion between times t 2k+1 and t 2k+2 and measurement errors 2k+2 and 2k+1 are independent of observations collected by time t 2k . Moreover, the process has a finite set of hidden states. This is an important property which is used to develop a forward algorithm for efficient computing of a composite likelihood in the next section. The conditional distribution of (Z(t 2k+2 ) − Z(t 2k+1 ), S(t 2k+2 )) (given the observations up to time t 2k ) depends only on state S(t 2k ).
Here we present our derivations for a one-dimensional case. The real-world animal movement data sets are two-dimensional. The formulas from below can be generalized to a ddimensional case. The details are given in Appendix. All our simulations and data analysis are also performed using two-dimensional formulas.

Composite Likelihood Estimation
Since the full likelihood is unavailable, we resort to composite likelihood to estimate the parameters (Lindsay, 1988). A composite likelihood is a weighted product of likelihood where L k is the true likelihood of the k-th data segment with a non-negative weight w k , k = 1, . . . , K, and K is the number of segments depending on the construction of the CL.
The weights are useful, for example, in pairwise likelihood when some pairs with stronger dependence contribute more than other pairs. Suppose that the location-state observations are denoted as The observed data only contain Z. We propose two ways to construct the composite likelihood.
4.1. Two-piece Composite Likelihood. The likelihood of increment-state observations at even numbered time points in which a is the largest integer not greater than a, θ = (λ 1 , λ 0 , σ, σ ) and ν(S(t 0 )) is the initial distribution that is assumed to be stationary. Since S even is not observed, we need to marginalize over all possible state trajectories: L(Z even ; θ) = S(t 0 ),S(t 2 ),...,S(t 2[n/2] )∈{0,1} L (Z even , S even ; θ) .
The cardinality of the set of the state trajectories is 2 [n/2]+1 , which makes the direct summation infeasible for even moderate n. It can, however, be tackled with the help of dynamic programming, specifically, by the forward algorithm.
First, let us define the forward variables by where Z even (t 2k ) = (Z(t 0 ), Z(t 2 ), . . . , Z(t 2k )), k = 1, . . . , [n/2], and the initial forward variable α (Z even (t 0 ), S(t 0 ), θ) = ν(S(t 0 )). Then, it is easy to see that the forward variables satisfy the following recursive relationship: This allows us to compute the likelihood in linear time with respect to n, because Now, when the sample size n is large, the series of multiplications may cause underflow problems where some terms are too small to be distinguished from zero by a computer. A normalized forward algorithm addresses the underflow problem. More specifically, let us introduce the normalized forward variables as α (Z even (t 2k ), S(t 2k ), θ) = α (Z even (t 2k ), S(t 2k ), θ) L (Z even (t 2k ); θ) , and let ρ (Z even (t 2k+2 ); θ) = L (Z even (t 2k+2 ); θ) L (Z even (t 2k ); θ) .
Each piece in (4) is a true log likelihood for about a half of the observations. The maximum composite likelihood estimator (MCLE) of θ is the maximizerθ of (4).

Marginal Composite Likelihood.
The second approach is to use the one-step transition density with the dependence between two consecutive increments ignored. If S were observed, for i, j ∈ {0, 1}, the likelihood of each pair of consecutive location-state observa- Since the dependence among the increments is discarded, the resulting estimator is expected to be less efficient if the dependence is stronger.

Variance Estimation.
To make inferences about θ, we need the variance ofθ. It can be estimated by parametric bootstrap with the time points fixed easily because simulating from the MRME process is simple. The general approach of parametric bootstrap is given as Algorithm 1.
Algorithm 1 Estimating standard error from parametric bootstrap input: Observed data; number of resampling M . · Fit model to get the parameter estimates; for m = 1 to M do · Use the estimated parameters to generate a bootstrap sample on the observed time grids; · Fit model to the bootstrap sample to get bootstrap estimate; end for · Make inference based on empirical distribution from these bootstrap estimators.

Simulation Study
We ran three simulation studies to check the performance of the MCLE based on both the marginal composite likelihood and the two-piece composite likelihood. The objective of the studies were threefold: (1) to see if the procedures successfully recovered the model parameters, (2) to verify that standard errors could be obtained with the help of parametric bootstrap, and (3) to compare the performance of the marginal method to the two-piece one.
In Study 1, we generated movement data using the MRME model described by equation (2). The model parameters were set to be λ 1 = 1, λ 0 = 0.5, σ = 1, and σ ∈ {0.01, 0.05}. This is the same setup that was used for simulations in Section 2. For each configuration, we generated 200 two-dimensional datasets on a time grid from 0 to 1000, with sampling interval 5. The resulting data has length n = 200. of MRME process were: λ 1 = 1, λ 0 = 0.5, σ = 1, and σ = 0.01. Table 2 (upper block) summarizes the results. Once again we can see that both the marginal method and the twopiece method recover the true parameters well. Their empirical standard errors are similar, suggesting that the two methods have comparable efficiency for these setups. Moreover, the standard errors were estimated by the parametric bootstrap procedure with 50 replications.
The estimated standard errors are reasonably close to empirical ones. The coverage rates of the 95% bootstrap confidence intervals are as low as 81% for λ 1 for the case with sampling interval 5 and sampling horizon 200. As the sampling interval decreases and the sampling horizon increases, the coverage rates get reasonably close to the nominal level.
Let us make a few remarks on the influence of sampling horizon and sampling interval on the efficiency of estimation. When the sampling interval is held fixed but the sampling horizon is 2.5 times longer, the ESE (empirical standard error) seems to be √ 2.5(= 1.58) times smaller for most parameter estimates. The longer sampling horizon covers more movingresting cycles (M j + R j ), and provides more information on both mobility and measurement error parameters. If we fix the sampling horizon and increase sampling frequency by reducing the sampling interval, however, the ESE of onlyσ reduces in proportion to the square root of the number of observations. Theoretically, if one can take observations nearly continuously, the mobility parameter σ can be estimated with absolute accuracy. Increasing sampling frequency also improves the estimation of λs to a certain degree, but improves the estimation of σ drastically. The results show the difference between the domain expansion asymptotics and the in-fill asymptotics.  Finally, let us note that the performance of the two methods is similar in both simulation studies described above. That was surprising because the marginal method basically ignores the dependence of the MRME process and treats increments as if they are independent.
One possible explanation is that when the distance between two consecutive observation is relatively long, then the dependence between two consecutive increments of the MRME process is weaker. To illustrate this point, one can calculate the correlation of absolute values of consecutive increments via simulation by employing the auto-correlation function with lag 1 (ACF(1)) and lag 2 (ACF(2)). For example, for the same parameter set as in the above simulations and a long time horizon 100,000, both ACF(1) and ACF(2) for the sampling interval 5 are very close to 0 from a Monte Carlo study. For the sampling interval 0.1, however, they are 0.46 and 0.40, respectively. Our third simulation study was based on this design. We also considered sampling interval 0.8 in the simulation because it has a relatively large ACF(1), 0.23, and a significantly smaller ACF(2), 0.07. The results of Study 3 with these small sampling intervals are presented in Table 2 (lower block). It is clear that the two-piece procedure is preferable for datasets with shorter sampling intervals (more frequent observations).

Movement of a Mountain Lion
The MRME model was applied to GPS data collected on a mature female mountain lion living in the Gros Ventre Mountain Range near Jackson, Wyoming. The data were collected by a code-only GPS wildlife tracking collar from 2009 to 2012. The collar was programmed to record locations every 8 hours, but the actual sampling intervals were irregular. In the Grand Teton and Gros Ventre mountains of Wyoming, deep winter snows ensure that mountain-lion movements differ across seasons (Elbroch et al., 2013). So, we fitted a MRME model to the summer data (from June 1, 2012 to August 31, 2012) and winter data (from December 1, 2011 to February 29, 2012), respectively. These two periods of data were plotted as Figure 3.
The summer data had an average sampling interval of 5.46 hours with standard deviation 5.14 hours, ranging from 0.5 hours to 32 hours. The average sampling interval was 5.58 hours in the winter data, with standard deviation 4.09 hours and range from 0.5 hours to 25 hours. The summer data has 401 observations and winter data has 392 observations.
The estimates based on two-piece composite likelihood for summer data areλ 1 = 2.841 hour −1 ,λ 0 = 0.179 hour −1 ,σ = 1.335 km/hour 1/2 andσ = 0.019 km. On average, this lion was moving 0.352 hours for each 5.587 hours resting, and, if she kept moving for (exactly) one full hour, the average deviation from the initial position is 1.335 km in both directions (northing and easting). Compared to the summer data, the estimates for winter data arê λ 1 = 6.225 hour −1 ,λ 0 = 0.118 hour −1 ,σ = 1.506 km/hour 1/2 andσ = 0.009 km, so, during the winter period, she spent 51.7% more time staying in place and 54.4% less time moving. The estimate of σ (standard deviation of Gaussian noise) indicates that the GPS tracking collar had about 10-to 20-meters of measurement error, and the error is twice as high in summer than in winter, which is consistent with the report by Owari et al. (2009).
However, the variability observed by Owari et al. (2009) was primarily due to sky obstruction from broadleaf tree canopy, which cannot be the case for the conifer-dominated mountains in Wyoming. Instead, these differences in GPS error likely reflect the fact that mountain lions in the study area select thicker vegetation with shade and cover in which to bed in the summer and more open terrain with southern aspects in rugged terrain, which catch sun and provide thermoregulatory benefits in the winter (Kusler et al., 2017). The results from the marginal method are similar (Table 3) except that the estimate of the rate parameter λ 1 of the moving state is noticeably smaller. Based on the comparison between the two methods in the simulation study, our discussion used the results from the two-piece method.

Discussion
Inactive periods and measurement errors are both necessary features of animal movement models, if we are to successfully reconstruct natural behaviors. Handling measurement errors in the MR process is especially critical because, if discarded, a microscopic amount of measurement error causes substantial bias in estimation. Our approach employing composite likelihood is the first to make the MRME model practically feasible. For movement data from predators that are known to have long inactive periods (Jeschke, 2007), the MRME model has great potential in revealing insights for animal ecologists. For example, the clear seasonal patterns in movement patterns have immediate implications for diverse ecological questions (e.g., seasonal foraging ecology, Elbroch et al., 2013) and estimating animal abundance using models dependent upon animal speed (Moeller et al., 2018).
The MRME model can be extended to meet further practical needs. For some predators, an inactive period may have different purposes, such as resting and food handling. Pozdnyakov et al. (2020) introduced the moving-resting-handling (MRH) process that allows two different types of motionless states. The moving period may also represent different behaviors (Benhamou, 2011;Kranstauber et al., 2012), and this can be accommodated by allowing the volatility parameter to have multiple levels. More generally, animal behavior depends on a suite of intrinsic and extrinsic variables. Our MRME model provides an important benchmark for building these more realistic extensions.
Then we get that Similarly, we also have Let us mention that these formulas do not require numerical multiple integral evaluation and, as a consequence, are computationally efficient.