<figure><img src="image_1.png" /></figure>

HAL
open science

# Proper account of long-term correlations in the observations improves state-space models' performances

Nicolas Bez, Pierre Gloaguen, Marie-Pierre Etienne, Rocio Joo, Sophie Lanco, Etienne Rivot, Emily Walker, Mathieu Woillez, Stéphanie Mahévas

► To cite this version:

Nicolas Bez, Pierre Gloaguen, Marie-Pierre Etienne, Rocio Joo, Sophie Lanco, et al. Proper account of long-term correlations in the observations improves state-space models' performances. 2024. hal-04547315

HAL Id: hal-04547315

https://hal.science/hal-04547315

Preprint submitted on 15 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Proper account of auto-correlations improves decoding performances of state-space (semi) Markov models

Nicolas Bez, Pierre Gloaguen, Marie-Pierre Etienne, Sophie Lanco, Rocío Joo, Etienne Rivot, Emily Walker, Mathieu Woillez, Stéphanie Mahévas

<figure><img src="image_1.png" /></figure>

## To cite this version:

Nicolas Bez, Pierre Gloaguen, Marie-Pierre Etienne, Sophie Lanco, Rocío Joo, et al. Proper account of auto-correlations improves decoding performances of state-space (semi) Markov models. Peer Community Journal, 2025, 5, pp.e38. 10.24072/pcjournal.535. hal-04547315v4

HAL Id: hal-04547315

https://hal.science/hal-04547315v4

Submitted on 26 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Proper account of long-term correlations in the observations improves state-space models' performances

Nicolas Bez<fnref n="1" />, Pierre Gloaguen<fnref n="2" />, Marie-Pierre Etienne<fnref n="3" />, Rocio Joo<fnref n="4" />, Sophie Lanco<fnref n="1" />, Etienne Rivot<fnref n="5" />, Emily Walker<fnref n="6" />, Mathieu Woillez<fnref n="7" />, & Stéphanie Mahévas<fnref n="1" />

**Correspondence:** nicolas.bez@ird.fr

<footnote n="1">¹ Marbec, IRD, Ifremer, Univ Montpellier, CNRS - Sète, France</footnote>

<footnote n="2">² Université de Bretagne Sud - Vannes, France</footnote>

<footnote n="3">³ IRMAR, Rennes University - Rennes, France</footnote>

<footnote n="4">⁴ Global Fishing Watch - Washington, USA</footnote>

<footnote n="5">⁵ Decod (Ecosystem Dynamics and Sustainability), L'Institut Agro, Ifremer, INRAE - Rennes, France</footnote>

<footnote n="6">⁶ BioSP, INRAE - Avignon, France</footnote>

<footnote n="7">⁷ Decod (Ecosystem Dynamics and Sustainability), L'Institut Agro, Ifremer, INRAE - Brest, France</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Proper account of auto-correlations
improves decoding performances of
state-space (semi) Markov models

Nicolas Bez<fnref n="1" />, Pierre Gloaguen<fnref n="2" />, Marie-Pierre Etienne<fnref n="3" />, Rocio Joo<fnref n="4" />, Sophie Lanco<fnref n="1" />, Etienne Rivot<fnref n="5" />, Emily Walker<fnref n="6" />, Mathieu Woillez<fnref n="7" />, & Stéphanie Mahévas<fnref n="1" />

<footnote n="1">¹ Marbec, IRD, Ifremer, Univ Montpellier, CNRS - Sète, France</footnote>

<footnote n="2">² Université de Bretagne Sud - Vannes, France</footnote>

<footnote n="3">³ IRMAR, Rennes University - Rennes, France</footnote>

<footnote n="4">⁴ Global Fishing Watch - Washington, USA</footnote>

<footnote n="5">⁵ Decod (Ecosystem Dynamics and Sustainability), L'Institut Agro, Ifremer, INRAE - Rennes, France</footnote>

<footnote n="6">⁶ BioSP, INRAE - Avignon, France</footnote>

<footnote n="7">⁷ Decod (Ecosystem Dynamics and Sustainability), L'Institut Agro, Ifremer, INRAE - Brest, France</footnote>
Abstract

State-space models are widely used in ecology to infer hidden behaviors. This study develops an exten-
sive numerical simulation-estimation experiment to evaluate the state decoding accuracy of four simple
state-space models. These models are obtained by combining different Markovian specifications (Markov
and semi-Markov) for the hidden layer with the absence and presence of auto-correlation for the observa-
tion layer. Model parameters are issued from two sets of real annotated trajectories. Three metrics are
developed to help interpret model performance. The first is the Hellinger distance between Markov and
semi-Markov sojourn time probability distributions. The second is sensitive to the overlap between the
probability density functions of state-dependent variables (e.g., speed variables). The third quantifies the
deterioration of the inference conditions between AR0 and AR1 formulations. It emerges that the most sen-
sitive model choice concerns the auto-correlation of the random processes describing the state-dependent
variables. Opting for the absence of auto-correlation in the model while the state-dependent variables are ac-
tually auto-correlated, is detrimental to state decoding performance. Regarding the hidden layer, imposing a
Markov structure while the state process is semi-Markov does not impair the state decoding performances.
The real-life estimates are consistent with our experimental finding that performance deteriorates when
there are significant temporal correlations that are not accounted for in the model. In light of these findings,
we recommend that researchers carefully consider the structure of the statistical model they suggest and
confirm its alignment with the process being modeled, especially when considering the auto-correlation of
observed variables.

**Keywords:** Hidden Markov models; Hidden semi-Markov models; state decoding accuracy; simulation-estimation; auto-regressive process; Ground truth

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22
Abstract

State-space models are widely used in ecology to infer hidden behaviors. This study develops an extensive numerical simulation-estimation experiment to evaluate the state decoding accuracy of four simple state-space models. These models are obtained by combining different Markovian specifications (Markov and semi-Markov) for the hidden layer with the absence (model AR0) and presence (AR1) of auto-correlation for the observation layer. Model parameters are issued from two sets of real annotated trajectories. Three metrics are developed to help interpret model performance. The first is the Hellinger distance between Markov and semi-Markov sojourn time probability distributions. The second is sensitive to the overlap between the probability density functions of state-dependent variables (e.g., speed variables). The third quantifies the deterioration of the inference conditions between AR0 and AR1 formulations. It emerges that the most sensitive model choice concerns the auto-correlation of the random processes describing the state-dependent variables. Opting for the absence of auto-correlation in the model while the state-dependent variables are actually auto-correlated, is detrimental to state decoding performance. Regarding the hidden layer, imposing a Markov structure while the state process is semi-Markov (with negative Binomial sojourn times) does not impair the state decoding performances. The real-life estimates are consistent with our experimental finding that performance deteriorates when there are significant temporal correlations that are not accounted for in the model. In light of these findings, we recommend that researchers carefully consider the structure of the statistical model they suggest and confirm its alignment with the process being modeled, especially when considering the auto-correlation of observed variables.

**Keywords:** Hidden Markov models; Hidden semi-Markov models; state decoding accuracy; simulation-estimation; auto-regressive process; Ground truth

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22
Introduction

The analysis of tracking data reveals information about the behavior underlying the movement patterns. After establishing the principle of representing behavior with a few meaningful states, segmentation methods can be used to estimate the states (Edelhoff et al., 2016). The commonly used segmentation methods includes filtering or thresholding (Gerritsen and Lordan, 2011), mixed models (Owen-Smith et al., 2010), K-means (Schwager et al., 2007), and hidden Markov models (HMM) or hidden semi-Markov models (HSMM) (Jonsen et al., 2013; Joo et al., 2013; Morales et al., 2004). This paper focuses on the two latter that allow time dependency to be considered (however, it is generally limited to first order, whereas data can have long temporal memories), and that are hierarchical models. Typically, they are composed of two layers: one to model the dynamics of the unobserved state process (also called the hidden process) and one to model the observed variables along the track (see, Zucchini and MacDonald, 2009 or McClintock, Langrock, et al., 2020 for didactic descriptions). Each of these layers can be flexible in their specifications (Gloaguen et al., 2015; Langrock et al., 2012; McClintock, King, et al., 2012; Morales et al., 2004). For example, one could use a Markov or semi-Markov hidden process to model the sequence of hidden states, and/or a correlated or uncorrelated random walk for the observed state-dependent variables.

When the states are known by direct or indirect observations (supervised cases), it is possible to infer properly the structure and the parameters of a segmentation model on a labeled training set, and to use the calibrated model to predict the states for a new unlabeled dataset. Learning methods, e.g., random forests (Sur et al., 2017; Thums et al., 2008) can also be used in these situations. However, when it comes to the most common unsupervised frameworks, model specification flexibility comes at a cost. Histograms for state-dependent variables and sojourn times in different states cannot be access to. As a result, one cannot infer the probabilistic laws they follow. In such cases, the *type* of probability density functions (PDF) of the state-dependent variables and the *type* of probability mass functions (PMF) of the sojourn times in the various states are not empirically grounded. They cannot be formally verified (Auger-Méthé, Clair, et al., 2011; Avgar et al., 2013; Pohle et al., 2017). However, postulating that the hidden process is Markov (HMM), with no time dependence in the observed state-dependent variables (AR0), may have a negative impact on the state decoding accuracy if the data do not conform with these choices. This study aims to assess the implications of some standard model specification choices.

In the absence of ground truth (unsupervised situations), a model can be evaluated against competing models using penalized maximum likelihood criteria such as AIC or BIC (Auger-Méthé, Newman, et al., 2021; Joo et al., 2013; Langrock et al., 2012; Pohle et al., 2017). Based on likelihood, selection is made from the standpoint of the visible side of the model, i.e., the observation layer. The hidden layers are not explicitly considered, even though they are the estimation's target. Another challenge is that several hidden state sequences can result in the same likelihood, making it difficult to select models based on likelihood. Zucchini and MacDonald, 2009 suggest to compare pseudo-residuals and the conditional probabilities to observe one particular observation knowing all observations, except the current one, in a Normal-Normal scales. Similarly Auger-Méthé, Newman, et al., 2021 use probability integral transform (PIT). However, these approaches only evaluate model consistency regarding the observed state-dependent variables. They do not permit evaluating model performance in terms of state estimations, which is the primary goal of these models. Simulation-estimation experiments are widely used in unsupervised situations to understand and evaluate the strengths and weaknesses of models (e.g., for the HMM context, McClintock, 2021). They are especially suitable for evaluating state decoding accuracy. Simulation-estimation entails assuming a true model, i.e., a model type and its parameters, which are then used to simulate sequences of states as well as realizations of the dependent variables (simulation step). The simulated dependent variables are then used as in real life unsupervised situations to re-estimate the parameters of the model, the *real* one or another model, and the states (estimation step).
Introduction

The analysis of tracking data reveals information about the behavior underlying the movement patterns. This might for instance concern the detection of fishing activities (e.g. Bertrand et al., 2007; Bez et al., 2011; Gerritsen and Lordan, 2011; Gloaguen et al., 2015; Vermard et al., 2010), as in the present study, or the identi- fication of exploratory versus encamped activities of wildlife animals, e.g. els (Morales et al., 2004) or bisons (Langrock et al., 2012). After establishing the principle of representing behavior with a few meaningful states, segmentation methods can be used to estimate the states (Edelhoff et al., 2016). The commonly used segmen- tation methods includes filtering or thresholding (Gerritsen and Lordan, 2011), mixture models (Owen-Smith et al., 2010), K-means (Schwager et al., 2007), and hidden Markov models (HMM) or hidden semi-Markov mod- els (HSMM) (Jonsen et al., 2013; Joo et al., 2013; Morales et al., 2004). This paper focuses on the two latter that allow time dependency to be considered (however, it is generally limited to first order, whereas data can have long temporal memories), and that are hierarchical models. Typically, they are composed of two layers: one to model the dynamics of the unobserved state process (also called the hidden process) and one to model the observed variables along the track (see, Zucchini and MacDonald, 2009 or McClintock, Langrock, et al., 2020 for didactic descriptions). Each of these layers can be flexible in their specifications (Gloaguen et al., 2015; Lan- grock et al., 2012; McClintock, King, et al., 2012; Morales et al., 2004). For example, one could use a Markov or semi-Markov hidden process to model the sequence of hidden states, and/or a correlated or uncorrelated random walk for the observed state-dependent variables.

When the states are known by direct or indirect observations (supervised cases), it is possible to infer prop-
erly the structure and the parameters of a segmentation model on a labeled training set, and to use the
calibrated model to predict the states for a new unlabeled dataset. Learning methods, e.g., random forests
(Sur et al., 2017; Thums et al., 2008) can also be used in these situations. However, when it comes to the
most common unsupervised frameworks, model specification flexibility comes at a cost. Histograms for state-
dependent variables and sojourn times in different states cannot be access to. As a result, one cannot infer
the probabilistic laws they follow. In such cases, the *type* of probability density functions (PDF) of the state-
dependent variables and the *type* of probability mass functions (PMF) of the sojourn times in the various states
are not empirically grounded. They cannot be formally verified (Auger-Méthé, Clair, et al., 2011; Avgar et al.,
2013; Pohle et al., 2017). However, postulating that the hidden process is Markov (HMM), with no time depen-
dence in the observed state-dependent variables (AR0), may have a negative impact on the state decoding
accuracy if the data do not conform with these choices. This study aims to assess the implications of some
standard model specification choices.

In the absence of ground truth (unsupervised situations), a model can be evaluated against competing mod-
els using penalized maximum likelihood criteria such as AIC or BIC (Auger-Méthé, Newman, et al., 2021; Joo et
al., 2013; Langrock et al., 2012; Pohle et al., 2017). Based on likelihood, selection is made from the standpoint
of the visible side of the model, i.e., the observation layer. The hidden layers are not explicitly considered, even
though they are the estimation's target. Another challenge is that several hidden state sequences can result
in the same likelihood, making it difficult to select models based on likelihood. Zucchini and MacDonald, 2009
suggest to compare pseudo-residuals and the conditional probabilities to observe one particular observation
knowing all observations, except the current one, in a Normal-Normal scales. Similarly Auger-Méthé, Newman,
et al., 2021 use probability integral transform (PIT). However, these approaches only evaluate model consis-
tency regarding the observed state-dependent variables. They do not permit evaluating model performance
in terms of state estimations, which is the primary goal of these models. Simulation-estimation experiments
are widely used in unsupervised situations to understand and evaluate the strengths and weaknesses of mod-
els (e.g., for the HMM context, McClintock, 2021). They are especially suitable for evaluating state decoding
accuracy. Simulation-estimation entails assuming a true model, i.e., a model type and its parameters, which
The use of state-space models in fisheries is motivated by the fact a significant proportion of vessels are required by regulation to be equipped with GPS systems that allow for the observation of their trajectory at regular time steps. However, we don't know when they're fishing. The use of state-space models proved to be efficient to estimate where and when a vessel is fishing knowing its trajectories, and thus to monitor fishing pressure on certain fish stocks (Bez et al., 2011; Vermard et al., 2010). This results however in uncontrolled estimation precision and underlines the need for an evaluation of the state decoding performances of such models. We were fortunate to have access to a few fishing vessel tracks where the fishing operations were meticulously documented. In this paper, we conduct simulation-estimation experiments with four types of models, the parameter values of which are derived from these sets of real trajectories, and compare the estimation of the states with the real ones. In addition to the fact that behavioral state durations must be longer than the observation time step, understanding an individual's behavior from his or her trajectory depends on the model's time step. This often coincides with the observation time step of the trajectory, which influences the temporal scope of correlations. Consequently, the ability of models to infer behavior is time-dependent, so it is useful to investigate the sensitivity of model performance to their temporal resolutions.

This study aims to assess state-decoding accuracy or the models' ability to recover true hidden states. Following a general introduction to the main characteristics of hidden Markov and semi-Markov models and to the three metrics defined in this study to compare them, we focus on four types of models by combining two Markovian structures for the hidden layer (Markov or semi-Markov) with two auto-regressive orders for the observation layer (AR0 or AR1). In line with the literature dealing with model skill assessment, i.e., on the ability to assess the fidelity of model behavior to truth (Allen and Somerfield, 2009; Lynch et al., 2009), we assess different facets of state decoding accuracy simulation-estimations experiments. We perform a Model Skill Assessment (MSA), where the model used for the simulation and the estimation is the same. Then, we carry out a Model Robustness Assessment (MRA) where the model used for estimation is one of all the possible models except the model used for simulation. Thanks to the ground-truther fishing behavior associated with the vessel trajectories, we finally complete the analyses by checking the model assumptions and, more importantly, evaluating the ability of the four types of models to recover the known states.

# Material and methods

## Definition of the simulation models

Hidden Markov (HMM) and hidden semi-Markov models (HSMM) are two similar state-space models formulated and implemented here in discrete regular-time frameworks. Table 1 lists all notations.

The hidden layer hosts a random process $S_t$ taking $N_S$ discrete possible values that can be considered as a Markov (resp. semi-Markov) process for HMM (resp. HSMM) state-space models. Given that we want to estimate one specific behavioral state rather than all of them, we consider two states ($N_S = 2$). This is also consistent with the data that are used to define the models' settings. In general, HMM cannot be considered a particular case of HSMM. As explained below, the only case in which an HMM can be considered a subset of an HSMM is when the probability mass functions (PMFs) for the HSMM's sojourn time are chosen from the negative binomial distributions.

The eventually multiple observation layers host the state-dependent variables. Two observed state-dependent variables are considered here ($N_V = 2$) corresponding to the average persistent and rotation speeds associated with a pair of consecutive locations and located by convention on the first of these points (Figure 1, Gloaguen et al., 2015). Formally, the properties of the average velocity are strongly dependent on the time delay between observations, so it is recommended in practice to take regular time steps. It also leads us to
are then used to simulate sequences of states as well as realizations of the dependent variables (simulation step). The simulated dependent variables are then used as in real life unsupervised situations to re-estimate the parameters of the model, the *real* one or another model, and the states (estimation step).

The use of state-space models in fisheries is motivated by the fact that a significant proportion of vessels are required by regulation to be equipped with GPS systems that allow for the observation of their trajectory at regular time steps. However, we don't know when they're fishing. The use of state-space models proved to be efficient to estimate where and when a vessel is fishing knowing its trajectories, and thus to monitor fishing pressure on certain fish stocks (Bez et al., 2011; Vermard et al., 2010). This results however in uncontrolled estimation precision and underlines the need for an evaluation of the state decoding performances of such models. We were fortunate to have access to a few fishing vessel tracks where the fishing operations were meticulously documented. In this paper, we conduct simulation-estimation experiments with four types of models, the parameter values of which are derived from these sets of real trajectories, and compare the estimation of the states with the real ones. In addition to the fact that behavioral state durations must be longer than the observation time step, understanding an individual's behavior from his or her trajectory depends on the model's time step. This often coincides with the observation time step of the trajectory, which influences the temporal scope of correlations. Consequently, the ability of models to infer behavior is time-dependent, so it is useful to investigate the sensitivity of model performance to their temporal resolutions.

This study aims to assess state-decoding accuracy or the models' ability to recover true hidden states. Following a general introduction to the main characteristics of hidden Markov and semi-Markov models and to the three metrics defined in this study to compare them, we focus on four types of models by combining two Markovian structures for the hidden layer (Markov or semi-Markov) with two auto-regressive orders for the observation layer (AR0 or AR1). In line with the literature dealing with model skill assessment, i.e., on the ability to assess the fidelity of model behavior to truth (Allen and Somerfield, 2009; Lynch et al., 2009), we assess different facets of state decoding accuracy simulation-estimations experiments. We perform a Model Skill Assessment (MSA), where the model used for the simulation and the estimation is the same. Then, we carry out a Model Robustness Assessment (MRA) where the model used for estimation is one of all the possible models except the model used for simulation. Thanks to the ground-truther fishing behavior associated with the vessel trajectories, we finally complete the analyses by checking the model assumptions and, more importantly, evaluating the ability of the four types of models to recover the known states.

# Material and methods

## Definition of the simulation models

Hidden Markov (HMM) and hidden semi-Markov models (HSMM) are two similar state-space models formulated and implemented here in discrete regular-time frameworks. Table 2 lists all notations.

The hidden layer hosts a random process $S_t$, indexed on time $t$, taking $N_S$ discrete possible values that can be considered as a Markov (resp. semi-Markov) process for HMM (resp. HSMM) state-space models. Given that we want to estimate one specific behavioral state rather than all of them, we consider two states ($N_S = 2$), a targeted state and the others. This is also consistent with the data that are used to define the models' settings. In general, HMM cannot be considered a particular case of HSMM. As explained below, the only case in which an HMM can be considered a subset of an HSMM is when the probability mass functions (PMFs) for the HSMM's sojourn time are chosen from the negative binomial distributions. The sojourn times in the various states are constraint to last at least 1 time unit. Knowing that the process enters, say, state $s$, at
investigate model performance as a function of time resolution.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Directed acyclic graph for HMM and HSMM. The black arrows indicate conditional dependencies. Conditional independence is thus reflected by the absence of an arrow. The four tested models (m = 1,...,4) get one hidden layer with two states ($N_S = 2$) and one observation layer with two variables ($N_V = 2$) whose definition is indicated on the top-left panel. Circles represent the two states (white for State 1 and plain for State 2). The observation layer is made of two random variables depicted by squares and triangles. Conditionally on the states, the variables of the observation layer can follow an AR0 (no arrow between them) or an AR1 (arrows). The number of parameters needed by layer assumed a shifted-negative binomial PMF for the sojourn times in the HSMM formulation and a Gaussian PDF for the state-dependent variables. $N_H$ and $N_O$ represent the number of parameters required to parameterize the hidden and the observation layers, respectively.</figcaption></figure>

## Definition and characterization of the hidden layers

For HMMs, the Markov property (of order 1) implies that $S_{t+1}$ is independent of ($S_1, S_2, ..., S_{t-1}$) conditionally on $S_t$, which means that:

$$[S_{t+1}|S_t = s_t, ..., S_1 = s_1] = [S_{t+1}|S_t = s_t]$$

where, $[X|Y = y]$ stands for the probability density function (PDF) of the random variable X, knowing that the random variable Y equals y. The stationary PMF of the sojourn times in the different states is defined by the transition matrix whose elements $p_{i,j}$ are the probabilities to switch from state i to state j when going from one time step to the next. A consequence of the conditional independence is that the sojourn time $T_s$ in each possible state $s = 1, ..., N_S$ follows a geometric PMF:

$$T_s \sim \text{geom}(1 - p_{s,s})$$
time $t + 1$, the distribution of the sojourn time $T_s$ in this state is defined as

$$P(T_s = k) = P(S_{t+k+1} \neq s, S_{t+k} = s, \dots, S_{t+2} = s | S_{t+1} = s, S_t \neq s), \text{ for } k \geq 1, s = 1, \dots, N_s$$

The eventually multiple observation layers host the state-dependent variables. Two observed state-dependent variables are considered here ($N_V = 2$) corresponding to the average persistent and rotation speeds associated with a pair of consecutive locations and located by convention on the first of these points (Figure 1, Gloaguen et al., 2015). Formally, the properties of the average velocity are strongly dependent on the time delay between observations, so it is recommended in practice to take regular time steps. It also leads us to investigate model performance as a function of time resolution.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Directed acyclic graph for HMM and HSMM. The black arrows indicate conditional dependencies. In HSMM, the time duration is an outcome of a random variable. The four tested models ($m = 1,...,4$) get one hidden layer with a two states ($N_S = 2$) variable and one observation layer with two variables ($N_V = 2$) whose definition is indicated on the top-left panel. Circles represent the states. The observation layer is made of two random variables depicted by squares and triangles. Conditionally on the state, the variables of the observation layer can follow an AR0 (no arrow between them) or an AR1 (arrows). The number of parameters needed by layer assumed a shifted-negative binomial PMF for the sojourn times in the HSMM formulation and a Gaussian PDF for the state-dependent variables. $N_H$ and $N_O$ represent the number of parameters required to parameterize the hidden and the observation layers, respectively.</figcaption></figure>

## Definition and characterization of the hidden layers

For HMMs, the state process $S_t$ is Markov which implies that $S_{t+1}$ is independent of $(S_1, S_2, \dots, S_{t-1})$ conditionally on $S_t$:

$$P(S_{t+1} = s_{t+1} | S_t = s_t, \dots, S_1 = s_1) = P(S_{t+1} = s_{t+1} | S_t = s_t), \forall t \geq 1$$
where, $p_{s,s}$ is the probability of remaining in state $s$ when going from one time step to the next. The first two moments are fully determined by $p_{s,s}$ as:

$$ \mathbb{E}(T_s) = \frac{1}{(1 - p_{s,s})} \\ \mathbb{V}(T_s) = \frac{p_{s,s}}{(1 - p_{s,s})^2} $$

and the number of parameters required to specify the hidden layer of an HMM is $N_H = N_S(N_S - 1)$.

For HSMMs, the hidden layer itself can be decomposed into two sub-layers (R Core Team (2022), Figure 1). The first sub-layer corresponds to a Markov chain, not indexed on time, defined by a transition matrix that drives the realization of the states' sequence (the transition matrix gets a zero diagonal and lines summing to one so that it gets $N_S(N_S - 2)$ elements). The second sub-layer specifies the PMF for the sojourn times in each state, not necessarily geometric and not necessarily parametric even though non-parametric formulations are rare (Johnson and Willsky, 2013; Sonia Malefaki and Limnios, 2010). In parametric cases, the hidden layer of an HSMM is defined by $N_H = N_S(N_S - 2) + N_S N_{p,S}$ parameters, where $N_{p,S}$ represents the number of parameters defining the PMFs of the sojourn time. In this work, the PMFs for the sojourn time are chosen to be shifted negative binomial ($s \mathcal{N} B$) functions ($N_{p,S} = 3$):

$$ T_s \sim \text{shift}_s + \mathcal{N}B(n_s, p_s) $$

with

$$ \mathbb{E}(T_s) = \text{shift}_s + n_s \frac{1 - p_s}{p_s} \\ \mathbb{V}(T_s) = n_s \frac{1 - p_s}{p_s^2} $$

The rationale for using such PMF is twofold. Firstly, in practice, sojourn times last at least one time step and cannot be null (shift=1). It is also common practice to obtain tracking data at a finer resolution than the duration of the states so that the smallest stay in a state is likely to be equivalent to a certain number of time steps (shift > 1 or >> 1). Secondly, geometrical PMFs, i.e. PMFs of the sojourn times in HMM cases, are particular cases of negative binomial PMFs ($n_s = 1$). Under this formulation, the HMM models are thus nested in the HSMM models (HMM $\subset$ HSMM$_{s\mathcal{N}B}$), which is not true in general for other HSMM specifications, i.e., for PMFs different from the family of the negative binomial functions. The distance between a shifted negative binomial PMF and its equivalent geometric PMF is evaluated by the Hellinger distances between the two PMFs (Hellinger, 1909). For each state $s = 1, \dots, N_S$, the Hellinger distance $d_{T_s}$ is given by:

$$ d_{T_s} = \sqrt{1 - \sum_{t=0}^{\infty} \sqrt{\text{geom}_s(t) \cdot s \mathcal{N} B_s(t)}} $$

Thus, the overall measure of distance between an HMM and an HSMM$_{s\mathcal{N}B}$ is:

$$ d_T = \sum_s \pi_s d_{T_s} $$

where $\pi_s = \mathbb{E}(T_s) / \sum_{s'} \mathbb{E}(T_{s'})$ is the probability of being in state $s$.

## Definition and characterization of the observation layer

**Probability density functions** As in standard practice, we consider that the state-dependent random variables $V_i, i = 1, 2$, are mutually independent ($V_{1,t} \perp V_{2,t'}$, $\forall t, t'$) with PDFs defined conditionally to the hidden process. In the present work, the two state-dependent variables are chosen univariate Gaussian ($N_{p0} = 2$):

$$ [V_{i,t} | S_t = s] \sim \mathcal{N}(\mu_{i,s}, \sigma_{i,s}) = \mathcal{N}_{i,s} $$
where $s_t$ is the outcome of $S_t$. The temporal random process is thus characterized by the transition matrix whose elements $p_{ss'}$ are the probabilities to switch from state $s$ to state $s'$ when going from one time step to the next:

$$p_{ss'} = P(S_{t+1} = s'|S_t = s), \forall t$$

Following this equation, $p_{ss}$ is the probability of remaining in state $s$ when going from one time step to the next, while $q_s = 1 - p_{ss}$ is the probability to leave state $s$.

A consequence of the conditional independence is that the distribution of the sojourn time $T_s$ in each possible state $s = 1, \dots, N_S$ follows a geometric PMF:

$$T_s \sim G(q_s)$$

In the present work, the sojourn times last at least 1 time step so that:

$$P(T_s = k) = g_s(k) = \begin{cases} 0 & \text{if } k = 0 \\ (1 - q_s)^{k-1} q_s & \text{if } k \ge 1 \end{cases}$$

The first two moments are fully determined by $q_s$ as:

$$\mathbb{E}(T_s) = \frac{1}{q_s}$$

$$\mathbb{V}(T_s) = \frac{1 - q_s}{q_s^2}$$

and the number of free parameters required to specify the hidden layer of an HMM is $N_H = N_S(N_S - 1)$.

HSMM have a Markov property, with the difference that the memoryless property does not act on time but on the jumps between the state process (Barbu and Limnios, 2008). In HMM the state at time $t+1$ only depends on the state at time $t$, irrespective of all the previous time steps. As mentioned above, in this case, the memory acts on time $t$ and $t+1$. In HSMM, the current state, that lasts a random number of time steps according to some PMF, only depends on the prior state that also lasts several time steps. The memory refers to the sequence of visited state irrespective of the time steps. In HSMM, the states follow a Markov chain of order one $\tilde{S}_n, n \in \mathbb{N}$, not indexed on time (hence the term *chain* rather than *process*.), so that:

$$P(\tilde{S}_{n+1} = s_{n+1} | \tilde{S}_n = s_n, \dots, \tilde{S}_1 = s_1) = P(\tilde{S}_{n+1} = s_{n+1} | \tilde{S}_n = s_n), \forall n$$

The Markov chain is characterized by the transition probabilities

$$\tilde{p}_{ss'} = P(\tilde{S}_{n+1} = s'|\tilde{S}_n = s)$$

The transition matrix gets a zero diagonal $\tilde{p}_{ss} = 0$, and lines sum to one, so that it requires $N_S(N_S - 2)$ free parameters (note that in the case $N_S = 2$, the matrix has no free parameters). The hidden layer of an HSMM is naturally indexed on the jumping times defined by the successive time points $J_n, n \in \mathbb{N}$ of states' changes (Barbu and Limnios, 2008; Figure 1). A semi-Markov process is such that the future visited state and the sojourn time in the present state depend only on the present state. In this case, the semi-Markov kernel is given by:

$$\tilde{p}_{ss'}(k) = P(\tilde{S}_{n+1} = s', \tilde{T}_{n+1} = k | \tilde{S}_n = s) = \tilde{p}_{ss'} \cdot P(\tilde{T}_{n+1} = k | \tilde{S}_{n+1} = s', \tilde{S}_n = s)$$

where $\tilde{T}_{n+1}$ and $\tilde{S}_{n+1}$ are the sojourn time and the state after the $n^{th}$ jump, i.e. starting at the $n^{th}$ jumping time $J_n$. In the present study, the HSMM follows a simplified structure where the conditioning is reduced to the current state Barbu and Limnios, 2008; page 46):
**Discrimination between PDFs** The estimation of the hidden states is a priori all the easier when the speed distributions per state have no or low overlaps. As an alternative to a direct measure of overlap or distance between PDFs, the discrepancy between the PDFs is quantified here with regard to the quality of the state prediction they potentially permit (the larger, the better). It corresponds to the state decoding accuracy achieved when ignoring the hidden layer, i.e. when considering a mixed model (Figure 2). Indeed, in this simple case, state estimation entails selecting the most likely state based on speed values. The state decoding accuracy when considering both variables, considered independent in the current case, is defined as:

$$d_V = \sum_s \pi_s \iint P(\hat{S} = s | S = s, V_1 = v_1, V_2 = v_2) \mathcal{N}_1(v_1) \mathcal{N}_2(v_2) dv_1 dv_2 \\ = \sum_s \pi_s \iint 1_{\mathcal{N}_{1,s}(v_1) \mathcal{N}_{2,s}(v_2) > \mathcal{N}_{1,s'}(v_1) \mathcal{N}_{2,s'}(v_2)} \mathcal{N}_{1,s}(v_1) \mathcal{N}_{2,s}(v_2) dv_1 dv_2$$

Using the empirical speed frequency distribution for variable i in state s yields empirical equivalents. These metrics, ranging from 0 to 1, serve as summary statistics for the observation layers of the various models used in this study.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Overlap between PDFs and accuracy of state decoding assuming a simple mixed model with one single state-dependent variable. Each PDF is Gaussian. Knowing the variable value, the probability of getting the correct state corresponds to the PDF that is maximum weighted by the probability of being in the corresponding state.</figcaption></figure>

**Auto-regressive process of order 0 (AR0)** It is also standard to consider that, conditionally to the state, $V_{i,t}$ and $V_{i,t+1}$ are not temporally correlated (auto-regressive process of order 0; AR0):

$$\operatorname{cor}(V_{i,t}, V_{i,t+1} | S_t = S_{k+1} = s) = \tau_{i,s} = 0$$

The number of parameters characterizing the PDF of the state-dependent variables, the observation layer of an AR0, is defined by $N_O = N_S N_V N_{pO}$ parameters.

**Auto-regressive processes of order 1 (AR1)** A natural extension is to consider auto-regressive processes of order 1 (AR1) for $V_{i,t}$ (Gloaguen et al., 2015; McClintock, Langrock, et al., 2020; Poritz, 1982). Conditionally to a given state, say $S_t = S_{t+1} = s$, the auto-regressive processes can be formulated as follow: $V_{i,t+1} = \eta_{i,s} + \tau_{i,s} V_{i,t} + \sigma_{i,s} \epsilon$, where $\epsilon \sim \mathcal{N}(0, 1)$ so that

$$[V_{i,t} | S_t = s] \sim \mathcal{N}(\mu_{i,s} = \frac{\eta_{i,s}}{1 - \tau_{i,s}}, \sigma_{i,s})$$

For AR1, the number of parameters characterizing the PDF of the state-dependent variables is $N_O = N_S N_V (N_{p0} + 172)$. The extra parameter corresponds to the auto-regressive coefficient. AR0 being a particular case of AR1
$$ \tilde{p}_{ss'}(k) = \tilde{p}_{ss'} \cdot P(\tilde{T}_{n+1} = k | \tilde{S}_{n+1} = s') $$

The relation between the (non Markov) stochastic process $S_t$ and the Markov chain $\tilde{S}_n$ is given by

$$ S_t = \tilde{S}_{N(t)} $$

where $N(t) := \max(n \in \mathbb{N} | J_n \le t)$ is the discrete-time counting process of the number of jumps in $[1, t]$. The state sequence over time is no longer a Markov chain, which relaxes the constraint on the PMF of the sojourn time which no longer needs be geometric. Non-parametric formulations are rare (Johnson and Willsky, 2013; Sonia Malefaki and Limnios, 2010) and are costly in terms of number of parameters. In this work, the PMFs for the sojourn time are chosen to be shifted negative binomial

$$ T_s \sim \mathcal{SNB}(\text{shift}_s, n_s, q_s) $$

defined by $N_{pS} = 3$ state specific parameters, the shift ($\text{shift}_s$) corresponding to the minimum sojourn time, and the two standard parameters of the negative binomial PMF ($n_s, q_s$). For $n_s \in \mathbb{R}_+$, the PMF is given by:

$$ P(T_s = k) = b_s(k) = \begin{cases} 0 & \text{if } 0 \le k < \text{shift}_s \\ \frac{\Gamma(k'+n_s)}{k'\Gamma(n_s)} q_s^{n_s} p_s^{k'} & \text{if } k \ge \text{shift}_s \text{ and where } k' = k - \text{shift}_s, \end{cases} $$

where $\Gamma(\cdot)$ is the gamma function. Expectation and variance are given by

$$ \mathbb{E}(T_s) = \text{shift}_s + n_s \frac{1 - q_s}{q_s}, $$

$$ \mathbb{V}(T_s) = n_s \frac{1 - q_s}{q_s^2}. $$

The rationale for using such PMF is twofold. First, in practice, sojourn times last at least one time step and cannot be null (shift $\ge 1$). It is also common practice to obtain tracking data at a finer resolution than the duration of the states so that the smallest stay in a state is likely to be equivalent to a certain number of time steps (shift > 1 or $>> 1$). Secondly, geometrical PMFs, i.e. PMFs of the sojourn times in HMM cases, are particular cases of negative binomial PMFs ($n_s = 1$). The $q_s$ parameter of the negative binomial PMF is related, but, in general, not equal to the equivalent parameter of the geometrical PMF. They are equal when $n_s = 1$ in which case, if the shift is also set to 1, the two PMFs are similar ($\mathcal{SNB}(1, 1, q_s) = \mathcal{G}(q_s)$). The hidden layer of an HSMM is thus defined by $N_H = N_S(N_S - 2) + N_S N_{pS}$ parameters.

Under this formulation, the HMM models are thus nested in the HSMM models ($HMM \subset HSMM_{\mathcal{SNB}}$). This is not true in general for other HSMM specifications, i.e., for PMFs different from the family of the negative binomial functions. The distance between a shifted negative binomial PMF and its equivalent geometric PMF is evaluated by the Hellinger distances between the two PMFs (Hellinger, 1909). For each state $s = 1, \dots, N_S$, the Hellinger distance $d_{T_s}$ is given by:

$$ d_{T_s} = \sqrt{1 - \sum_{k=0}^{\infty} \sqrt{g_s(k) \cdot b_s(k)}} $$

Thus, the overall measure of distance between an HMM and an HSMM$_{\mathcal{SNB}}$ is:

$$ d_T = \sum_s \Pi_s d_{T_s} $$

where $\Pi_s = \mathbb{E}(T_s) / \sum_{s'} \mathbb{E}(T_{s'})$ is the probability of being in state $s$ given by the stationary (i.e. invariant) distribution of the Markov chain.
(AR0 ⊂ AR1), it is possible to estimate an AR1 on the basis of AR0 data expecting that the coefficient of correlation will be 0 ($\tau_{i,s} = 0$), while the reverse is unwise.

The discrete auto-correlation function of an AR1 writes $\rho(t) = \tau^t$. The correlation time (Yaglom, 1987) or the integral range (Lantuéjoul, 1991) is defined by $A = \sum_{t=0}^{\infty} \tau^t$. It is homogeneous to a time duration and characterizes the intrinsic temporal scale of a process ($A = 1$ unit of time for a nugget effect, $A = 2$ when $\tau = 0.5$ and $A = 10$ when $\tau = 0.9$). In terms of prediction, it contributes to define if a given observation time window can be regarded as large enough (resp. too low) for the mean of one realization to represent a precise prediction the mean of the process. This says that the time window of observation must be 10 times larger for an AR1 with $\tau = 0.9$ than for an AR0 to get the same level of precision when estimating the mean of the process. Reversely, when the time window is fixed, the precision deteriorates with increasing correlation coefficients.

The metric used to measure the difference between the AR1 and AR0 formulations is defined as the difference between the ratio of the integral range to the mean time duration of the various states. For an AR0 state-space model, the integral range is 1, and the ratio is simply $1/\mathbb{E}(T_s)$. The unit-less differences by variable and state ($d_{AR_{i,s}}$), by variable ($d_{AR_s}$) and for all variables and states ($d_{AR}$) between AR1 and AR0 formulations of a given state-space models are defined in this study by:

$$d_{AR_{i,s}} = \frac{A_{i,s} - 1}{\mathbb{E}(T_s)} \\ d_{AR_s} = \frac{1}{N_V} \sum_i d_{AR_{i,s}} \\ d_{AR} = \sum_s \pi_s d_{AR_s}$$

## Definition of the simulation models

**Four model types** Four types of models are investigated (Figure 1) by crossing the two types of Markov structures with the two types of observation processes:

$m_1 = \text{HMM-AR0}$

$m_3 = \text{HSMM-AR0}$

$m_2 = \text{HMM-AR1}$

$m_4 = \text{HSMM-AR1}$

The number of parameters ($N_\theta = N_H + N_O$) ranges from 10 to 18 depending on the model.

**Model parameters and time resolutions** Parameters are set using ground-truthed tracks from two different fishing vessels with distinct fishing behaviors. Vessel 1 is a trawler operating in the North-East Atlantic with onboard sensors that detect fishing operations using GPS acquisition at 0.001 Hz (every 15 minutes). Vessel 2 is a purse seiner fishing in the Indian Ocean, with observers on board who record when the vessel is fishing and GPS acquisition at 0.003 Hz (every 5 minutes). For both vessels, the reported states are non-fishing (State 1) and fishing (State 2). Persistent and rotation speeds are calculated at four different time resolutions by degrading observations at 30 minutes, 1 hour, and 2 hours for Vessel 1 and 15 minutes, 30 minutes, and 1 hour for Vessel 2, respectively. To facilitate model comparison and gain generality, the resolution is expressed as the average number of time steps in the fishing state (i.e., State 2). For example, for Vessel 2, if a fishing operation lasts an hour and a half on average, the four resolutions correspond to 18, 6, 3, and 1.5 observations per average fishing operation.

For simulation-estimation experiments, parameter coherence across the four model types is critical. For the AR1 formulations, the parameters $\eta_{i,s}$ are deduced from $\tau_{i,s}$ and $\mu_{i,s}$ ensuring that the expected values and variances of the respective PDFs are the same in AR0 and AR1 for a given resolution. For the sojourn time PMFs, the parameters of the shifted negative binomial PMFs (shift,n,p) are evaluated by maximum likelihood based on observed tracking data. The HMM formulation's parameters are then deduced to ensure that the
## Definition and characterization of the observation layer

**Probability density functions** As in standard practice, we consider that the state-dependent random variables $V_i, i = 1, 2$, are mutually independent conditionally to the state sequence with PDFs defined conditionally to the hidden process. In the present work, the two state-dependent variables are chosen univariate Gaussian ($N_{pO} = 2$):

$$[V_{i,t} | S_t = s] \sim \mathcal{N}(\mu_{i,s}, \sigma_{i,s}).$$

Their PDFs defined conditionally to the hidden process are thus denoted:

$$f_{V_i|S_t=s}(v) = \mathcal{N}_{i,s}(v)$$

**Discrimination between PDFs** The estimation of the hidden states is a priori all the easier when the speed distributions per state have no or low overlaps. As an alternative to a direct measure of overlap or distance between PDFs, the discrepancy between the PDFs is quantified here with regard to the quality of the state prediction they potentially permit (the larger, the better). It corresponds to the state decoding accuracy achieved when considering a mixture model (Figure 2). Indeed, in this simple case, the state is no longer considered as a stochastic process but as a random variable $S$ and the state estimator $\hat{S}$ entails selecting the most likely state based on speed values, following the classical Bayes rule. Then, when all parameters are known, the state estimator derived from a given velocity $V = (v_1, v_2)$ is given by:

$$\hat{S} = \underset{s}{\operatorname{argmax}} \left( \prod_s \mathcal{N}_{1,s}(v_1) \mathcal{N}_{2,s}(v_2) \right)$$

For each $1 \le s \le N_s$, we define $R_s = \{V = (v_1, v_2), \text{ such that } \hat{S} = s\}$, i.e., the region of $\mathbb{R}^2$ where the state is estimated as $s$. The state decoding accuracy derived from the overlap between the PDFs of the speed velocities under a mixture model is then defined as:

$$\begin{aligned}
d_V^{\text{mod}} &= \sum_s P(\hat{S} = s, S = s) \\
&= \sum_s \prod_s P(\hat{S} = s | S = s) \\
&= \sum_s \prod_s \int_{R_s} \mathcal{N}_{1,s}(v_1) \mathcal{N}_{2,s}(v_2) dv_1 dv_2.
\end{aligned}$$

When using labeled data, one can compute the empirical histograms of both speed variables by state and compute the following empirical equivalent:

$$d_V^{\text{emp}} = \sum_s \pi_s \sum_{(v_1, v_2) \in r_s} h_{1,s}(v_1) h_{2,s}(v_2),$$

where $\pi_s$ represents the empirical proportion of time effectively spent in state $s$, $\pi_s = \frac{\#\{s_i=s\}_{1,...,T}}{T}$, and $r_s$ is the empirical version of $R_s$.

These metrics, ranging from 0 to 1, will be used as summary statistics for the observation layers of the various real cases and models used in this study.

**Auto-regressive process of order 0 (AR0)** It is also standard to consider that, conditionally to the state, $V_{i,t}$ and $V_{i,t+1}$ are not temporally correlated (auto-regressive process of order 0; AR0):

$$\operatorname{cor}(V_{i,t}, V_{i,t+1} | S_t = S_{t+1} = s) = \tau_{i,s} = 0$$

The number of parameters characterizing the PDF of the state-dependent variables, the observation layer of an AR0, is defined by $N_O = N_S N_V N_{pO}$ parameters.
mean durations are the same in both frameworks:

$$p_{s,s} = 1 - \frac{1}{shift_s + n_s \frac{1-p_s}{p_s}}$$

The variance is automatically deduced. The variances in state durations between the HSMM and HMM formulations may thus differ. Settings 1 (resp. 2) denote the sets of model parameters estimated from the tracks of Vessel 1 (resp. 2). Tables 2 and 3 list the model parameters.

## Simulation-estimation experiments

### Dimensions of the numerical experiments

Simulation-estimation experiments (Figure 3) assess the accuracy and robustness of state decoding under controlled conditions for:

$$2 \text{ settings} \times 4 \text{ resolutions} \times 4 \text{ model types} = 32 \text{ different cases}$$

For each model, 100 simulation-estimation experiments are run. In each experiment, synthetic data (states and speeds) are simulated, and the model parameters ($\hat{\theta}$) and states are re-estimated using expectation-maximization algorithms (refer to below). Each simulation lasts long enough to produce 250 state changes. We also carried out experiments with 50 alternations.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Workflow illustrated for one case study (grey panels represent the workflow for one particular resolution). **Upper left panels:** observed time series in its original high-frequency resolution (r=1) and in three levels of up-scaling (r = 2,3,4). **Upper right panels:** computation of persistent and rotational speeds ($v_1, v_2$) for each resolution and estimation by maximum likelihood procedures of models’ parameters and settings ($\theta$). **Bottom Right panels:** the simulation-estimation experiments. For each tested model structure ($m_1, m_2, m_3$ and $m_4$), 100 data simulations are produced. For each simulation, the estimation is made using the EM algorithm with the same model (model skill assessments-MSA) or with a different model (model robustness assessment-MRA). Experiments get 250 or 50 alternations of states 1 and 2. **Bottom left panel:** estimation from the annotated real data with model structure $m_4$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Overlap between PDFs and accuracy of state decoding assuming a simple mixture model with one single state-dependent variable. Each PDF is Gaussian. Knowing the variable value, the probability of getting the correct state corresponds to the PDF that is maximum weighted by the probability of being in the corresponding state.</figcaption></figure>

**Auto-regressive processes of order 1 (AR1)** A natural extension is to consider auto-regressive processes of order 1 (AR1) for $V_{i,t}$ (Gloaguen et al., 2015; McClintock, Langrock, et al., 2020; Poritz, 1982). Conditionally to a given state, say $S_t = S_{t+1} = s$, the auto-regressive processes can be formulated as follow: $V_{i,t+1} = \eta_{i,s} + \tau_{i,s}V_{i,t} + \sigma_{i,s}\epsilon_t$, where $\epsilon_t \sim \mathcal{N}(0, 1)$ so that

$$[V_{i,t+1}|S_{t+1} = S_t = s] \sim \mathcal{N}(\mu_{i,s} = \frac{\eta_{i,s}}{1 - \tau_{i,s}}, \sigma_{i,s})$$

For AR1, the number of parameters characterizing the PDF of the state-dependent variables is $N_O = N_S N_V (N_{p0} + 1)$. The extra parameter corresponds to the auto-regressive coefficient. AR0 being a particular case of AR1 (AR0 $\subset$ AR1), it is possible to estimate an AR1 on the basis of AR0 data expecting that the coefficient of correlation will be 0 ($\tau_{i,s} = 0$), while the reverse is unwise.

The discrete auto-correlation function of an AR1 writes $\rho(t) = \tau^t$. The correlation time (Yaglom, 1987) or the integral range (Lantuéjoul, 1991) is defined by $A = \sum_{t=0}^{\infty} \tau^t$. It is homogeneous to a time duration and characterizes the intrinsic temporal scale of a process ($A = 1$ unit of time for a nugget effect, $A = 2$ when $\tau = 0.5$ and $A = 10$ when $\tau = 0.9$). In terms of prediction, it contributes to define if a given observation time window can be regarded as large enough (resp. too low) for the mean of one realization to represent a precise prediction the mean of the process. This says that the time window of observation must be 10 times larger for an AR1 with $\tau = 0.9$ than for an AR0 to get the same level of precision when estimating the mean of the process. Reversely, when the time window is fixed, the precision deteriorates with increasing correlation coefficients.

The metric used to measure the difference between the AR1 and AR0 formulations is defined as the difference between the ratio of the integral range to the mean time duration of the various states. For an AR0 state-space model, the integral range is 1, and the ratio is simply $1/\mathbb{E}(T_s)$. The unit-less differences by variable and state ($d_{AR_{i,s}}$), by state ($d_{AR_s}$) and for all variables and states ($d_{AR}$) between AR1 and AR0 formulations of a given state-space models are defined in this study by:

$$d_{AR_{i,s}} = \frac{A_{i,s} - 1}{\mathbb{E}(T_s)}$$
$$d_{AR_s} = \frac{1}{N_V} \sum_i d_{AR_{i,s}}$$
$$d_{AR} = \sum_s \pi_s d_{AR_s}$$
## Maximum likelihood inference

Inferences are made by maximizing the likelihood using Expectation-Maximization (EM) algorithms (see Dempster et al., 1977 for HMM and Guédon, 2003 in the context of HSMM). EM algorithms are iterative and converge to the (possibly local) maximum likelihood. When the shape of the likelihood is unknown, it is recommended to run the algorithm from a variety of starting points. Inferences were made using starting values equal to ¼, ½, 1, 2, or 4 times the input parameter values used in simulations. Convergence occurred at the same values regardless of starting values, indicating that, in our case, the likelihood was unimodal, with a tractable optimum free of local and detrimental holes. To avoid mixing states, label switching (assigning State 1 to slow speed instead of large speed, and vice versa) is resolved prior to analyzing the results. The state estimate represents the maximum a posteriori at each time step. Alternatively, the Viterbi algorithm (Rabiner, 1989) can reconstruct the most likely sequence of states. In the current study, both outputs are very similar, so we only present the former.

## State decoding accuracies

**General definition** The capacity to recover the model parameters ($\theta$) is investigated using the box-plots of the hundred estimations of each of the model parameters. However, the performance is mainly based on state decoding accuracies as the objectives of the recourse to these models are to estimate states. Decoding accuracy is defined by the probability of correctly estimating the states (proportions of true positives and true negatives):

$$ accuracy = P(\hat{S} = S) = \sum_{s} \pi_{s} P(\hat{S} = s | S = s) $$

Two facets of the accuracy are considered (Figure 3).

**Model skill assessment (MSA)** In sensitivity analysis (Allen and Somerfield, 2009), model skill assessment refers to scenarios where the simulation and estimation models are identical. This allows evaluating accuracy in the ideal conditions, that is when the data and model type are completely consistent (best case scenario):

$$ MSA(m_i, m_j) = P(\hat{S} = S | \text{simulation model} = m_i, \text{estimation model} = m_j) $$

**Model robustness assessment (MRA)** In reference to the above MSA, the model robustness assessment (MRA) considers accuracies when the estimation model differs from the model used for simulation:

$$ MRA(m_i, m_j) = P(\hat{S} = S | \text{simulation model} = m_i, \text{estimation model} = m_j) $$

For instance, if data are simulated according to $m_4$ = HSMM-AR1, estimations are performed with each of the three other possible model types, one by one ($m_1, m_2$, or $m_3$). An analysis of the impact of discrepancies between the Markov property of the hidden layer and/or the AR properties of the state-dependent layer allows testing which of the two discrepancies affects the most the percentage of loss of accuracy when using model $m'$ instead of model $m$:

$$ loss(m_i, m_j) = \frac{MRA(m_i, m_j) - MSA(m_i, m_j)}{MSA(m_i, m_j)} $$

## Real case state decoding performance

HSMM-AR1, the most flexible model type, is fitted to two sets of real annotated trajectories (Figure 3) at different time resolutions. Real states are used to calculate the associated accuracy. The most likely order of auto-regressive processes used to model the speed time series per state is an important piece of information for interpreting real-world decoding performance. To test if a (partial) coefficient of correlation of order $k$ is null, the appropriate Student test statistics is:

$$ \frac{R}{\sqrt{\frac{1-R^2}{n-k-2}}} \sim T_{n-k-2} $$
## Definition of the simulation models

**Four model types** Four types of models are investigated (Figure 1) by crossing the two types of Markov structures with the two types of observation processes:

$$m_1 = \text{HMM-AR0}$$

$$m_2 = \text{HMM-AR1}$$

$$m_3 = \text{HSMM-AR0}$$

$$m_4 = \text{HSMM-AR1}$$

The number of parameters ($N_\theta = N_H + N_O$) ranges from 10 to 18 depending on the model.

**Model parameters and time resolutions** Parameters are set using ground-truthed tracks from two different fishing vessels with distinct fishing behaviors. Vessel 1 is a trawler operating in the North-East Atlantic with onboard sensors that detect fishing operations using GPS acquisition at 0.001 Hz (every 15 minutes). Vessel 2 is a purse seiner fishing in the Indian Ocean, with observers on board who record when the vessel is fishing and GPS acquisition at 0.003 Hz (every 5 minutes). For both vessels, the reported states are non-fishing (State 1) and fishing (State 2). Persistent and rotation speeds are calculated at four different time resolutions by degrading observations at 30 minutes, 1 hour, and 2 hours for Vessel 1 and 15 minutes, 30 minutes, and 1 hour for Vessel 2, respectively. To facilitate model comparison and gain generality, the resolution is expressed as the average number of time steps in the fishing state (i.e., State 2). For example, for Vessel 2, if a fishing operation lasts an hour and a half on average, the four resolutions correspond to 18, 6, 3, and 1.5 observations per average fishing operation.

For simulation-estimation experiments, parameter coherence across the four model types is critical. For the AR1 formulations, the parameters $n_{i,s}$ are deduced from $\tau_{i,s}$ and $\mu_{i,s}$ ensuring that the expected values and variances of the respective PDFs are the same in AR0 and AR1 for a given resolution. For the sojourn time PMFs, the parameters of the shifted negative binomial PMFs ($shift,n,p$) are evaluated by maximum likelihood based on observed tracking data. The HMM formulation's parameters are then deduced to ensure that the mean durations are the same in both frameworks:

$$p_{ss} = 1 - \frac{1}{shift_s + n_s \frac{1-p_s}{p_s}}$$

The variance is automatically deduced. The variances in state durations between the HSMM and HMM formulations may thus differ. Settings 1 (resp. 2) denote the sets of model parameters estimated from the tracks of Vessel 1 (resp. 2). Tables 3 and 4 list the model parameters.

## Simulation-estimation experiments

### Dimensions of the numerical experiments

Simulation-estimation experiments (Figure 3) assess the accuracy and robustness of state decoding under controlled conditions for:

$$2 \text{ settings} \times 4 \text{ resolutions} \times 4 \text{ model types} = 32 \text{ different cases}$$

For each case, 100 simulation-estimation experiments are run. In each experiment, synthetic data (states and speeds) are simulated, and the model parameters ($\hat{\theta}$) and states are re-estimated using expectation-maximization algorithms (refer to below). Each simulation lasts long enough to produce 250 state changes. We also carried out experiments with 50 alternations.

### Maximum likelihood inference

Inferences are made by maximizing the likelihood using Expectation-Maximization (EM) algorithms (see Dempster et al., 1977 for HMM and Guédon, 2003 in the context of HSMM). EM algorithms are iterative and
where *n* represents the number of observations so that *n* − *k* represents the number of pairs of data *k* apart in time.

In this study, the time series are interrupted by each change of state, and the lengths of each section are variables (Figure 4). So the number of data available to test the value of partial coefficients of correlation of order *k* is rather *n*<sub>*k*</sub>, the number of pairs of observations available at time distance *k* accounting for the gaps in the time series. For each state-dependent variable and each state, the order of the auto-regressive process is estimated as the smallest order for which the empirical partial correlation falls within the 90% confidence intervals of an *T*(*n*<sub>*k*</sub> − 2).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Adaptation of the Student T-test to account for the effective number of available data *n*<sub>*k*</sub> for each tested order *k* for an auto-regressive process.</figcaption></figure>

Results

Characteristics of the simulation models

Comparing settings

Figure 5 depicts the features that distinguish the two settings. In both settings, *d*<sub>*T*</sub> and *d*<sub>AR</sub> are co-linear, and they increase when the resolutions decrease (up to the exception of *d*<sub>*T*</sub> for Settings 2, which is between the two highest resolutions). In other words, deviations from the HMM reference situation result in proportional deviations from the AR0 reference situation, which are significantly larger for Settings 1 than for Settings 2. Fitting an HMM-AR0 model to data simulated with an HSMM<sub>S,N,B</sub>-AR1 model is expected to be challenging, especially at high resolutions. The reverse is not a problem because AR0 and/or HMM models are special cases of AR1 and/or HSMM<sub>S,N,B</sub> models. The proportions of time spent in states 1 and 2 are strongly unbalanced, with opposite dominance between settings. The state-specific characteristics of the various settings are thus weighted in the opposite direction. The average sojourn time in the various states also differs significantly between the two settings. In State 1, the average duration is ten times shorter in Settings 1 than in Settings 2 (regardless of resolution). When compared to State 2, it is 1.5 times larger. Despite the similarity of the correlation coefficients, the observed difference in *d*<sub>AR</sub> can be attributed to this significant difference between the two settings. In State 1, the average correlation coefficients and integral ranges are very similar for the two settings (Figure 6, right column). We are using averages to explain the overall patterns; it is worth noting that the coefficients of correlation increase with increasing resolution. Despite similar coefficients of correlation, the distance *d*<sub>AR1</sub> for Settings 1 is 10 times larger due to their longer sojourn times. The picture for state two is different. Coefficients of correlation and integral ranges are slightly smaller in Settings 2. The slightly longer sojourn time results in a smaller *d*<sub>AR2</sub> for Settings 2. Therefore, *d*<sub>AR1</sub> and *d*<sub>AR2</sub> are larger for Settings 1. The unbalanced states' proportions between settings resulted in a larger *d*<sub>AR</sub> for Settings 1 but reduced the
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Workflow illustrated for one case study (grey panels represent the workflow for one particular resolution). **Upper left panels:** observed time series in its original high-frequency resolution (r=1) and in three levels of up-scaling (r = 2,3,4). **Upper right panels:** computation of persistent and rotational speeds ($v_1, v_2$) for each resolution and estimation by maximum likelihood procedures of models’ parameters and settings ($\theta$). **Bottom Right panels:** the simulation-estimation experiments. For each tested model structure ($m_1, m_2, m_3$ and $m_4$), 100 data simulations are produced. For each simulation, the estimation is made using the EM algorithm with the same model (model skill assessments-MSA) or with a different model (model robustness assessment-MRA). Experiments get 250 or 50 alternations of states 1 and 2. **Bottom left panel:** estimation from the annotated real data with model structure $m_4$.</figcaption></figure>

converge to the (possibly local) maximum likelihood. When the shape of the likelihood is unknown, it is recommended to run the algorithm from a variety of starting points. Inferences were made using starting values equal to ¼, ½, 1, 2, or 4 times the input parameter values used in simulations. Convergence occurred at the same values regardless of starting values, indicating that, in our case, the likelihood was unimodal, with a tractable optimum free of local and detrimental holes. To avoid mixing states, label switching (assigning State 1 to slow speed instead of large speed, and vice versa) is resolved prior to analyzing the results. The state estimate is the one that ensures the maximum likelihood of the state knowing the speeds at each time step:

$$ \hat{s}_t = \underset{s}{\operatorname{argmax}} P(S_t = s | V_{1,...,T} = v_{1,...,T}, V_{2,...,T} = v_{2,...,T}) $$

Alternatively, the Viterbi algorithm (Rabiner, 1989) can reconstruct the most likely sequence of states

$$ \{\hat{s}_t\}_{t=1...T} = \underset{s_1...s_T}{\operatorname{argmax}} P(S_1 = s_1, \dots, S_T = s_T | V_{1,...,T} = v_{1,...,T}, V_{2,...,T} = v_{2,...,T}) $$

Consistent with Hoek and Elliott, 2018 (p. 129), in the current study, the transition matrix gets no 0 values explaining why both outputs are very similar. So we only present the former.

## State decoding accuracies

**General definition** The capacity to recover the model parameters ($\theta$) is investigated using the box-plots of the hundred estimations of each of the model parameters. However, the performance is mainly based on
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Graphical representation of the model settings. Settings 1 and 2 are summarized by $d_T$ the Hellinger distances between the shifted-negative binomial (HSMM) and the geometric (HMM) PMF formulations for the sojourn time on the x-axis, by $d_{AR}$ the difference between AR1 and AR0 formulations on the y-axis, and by $d_V$ the state decoding accuracies in case of mixed model through the surface area of the circles. The proportions of time spent in each state are also represented by two half discs within the circle. The arrows in each setting indicate the path towards decreasing resolution.</figcaption></figure>

Finally, given the overlaps between the PDFs of the speed variables by states, and the expected values of the different sojourn times, the expected state decoding capacities when using mixed models ($d_V$) are larger for Settings 2 than for Settings 1 (Figure 5), with 95% of good predictions for the largest resolution for Settings 2 as opposed to 80% for Settings 1 (Figure 9). All of these statements point to disappointing decoding accuracy for Settings 1 and decreasing resolutions in general. This latter statement suffers one exception: for Settings 1; the accuracy appears slightly better for the lowest resolution than the penultimate resolution (Figure 9). The expected values of the PMFs are equal in HMM and HSMM<sub>$sN^B$</sub> formulations, but the standard deviations differ. For Settings 1 and 2, the standard deviations are co-linear between the two types of Markov formulations but are generally larger in the HMMs (Figure 7). It also appears that, for Settings 1, the standard deviations of state durations in States 1 and 2 are of the same order of magnitude, whereas for Settings 2, they are one order of magnitude larger in State 1 than in State 2.

MSA (estimation model = simulation model)

Parameter estimation

In general, parameter estimates under MSA conditions are unbiased and symmetrically distributed around the target value (Figure 8). One notable exception is the variance of the Gaussian PDF for $V_{1,2}$ for Settings 2 when using an HSMM<sub>$sN^B$</sub>-AR0 model. When the target value is 1.26, 50% of the estimations fall between $Q_{25\%} = 1.25$ and $Q_{75\%} = 1.44$, indicating a minor issue. The variance of the shifted negative binomial at
state decoding accuracies as the objectives of the recourse to these models are to estimate states. Decoding accuracy is defined by the proportion of correct estimations:

$$ accuracy = p(\hat{s} = s) = \frac{\#\{\hat{s}_t = s_t\}_{t=1...T}}{T} $$

Two facets of the accuracy are considered (Figure 3).

**Model skill assessment (MSA)** In sensitivity analysis (Allen and Somerfield, 2009), model skill assessment refers to scenarios where the simulation and estimation models are identical. This allows evaluating accuracy in the ideal conditions, that is when the data and model type are completely consistent (best case scenario):

$$ MSA(m_i, m_j) = p(\hat{s} = s | \text{simulation model} = m_i, \text{estimation model} = m_j) $$

**Model robustness assessment (MRA)** In reference to the above MSA, the model robustness assessment (MRA) considers accuracies when the estimation model differs from the model used for simulation:

$$ MRA(m_i, m_j) = p(\hat{s} = s | \text{simulation model} = m_i, \text{estimation model} = m_j) $$

For instance, if data are simulated according to $m_4$ = HSMM-AR1, estimations are performed with each of the three other possible model types, one by one ($m_1, m_2$, or $m_3$). An analysis of the impact of discrepancies between the Markov property of the hidden layer and/or the AR properties of the state-dependent layer allows testing which of the two discrepancies affects the most the percentage of loss of accuracy when using model $m'$ instead of model $m$:

$$ loss(m_i, m_j) = \frac{MRA(m_i, m_j) - MSA(m_i, m_j)}{MSA(m_i, m_j)} = \frac{MRA(m_i, m_j)}{MSA(m_i, m_j)} - 1 $$

## Real case state decoding performance

HSMM-AR1, the most flexible model type, is fitted to two sets of real annotated trajectories (Figure 3) at different time resolutions. Real states are used to calculate the associated accuracy. The most likely order of auto-regressive processes used to model the speed time series per state is an important piece of information for interpreting real-world decoding performance. To test if $\rho_k$, the various (partial) coefficient of correlation of order $k$ (Saporta, 1990, p. 137) are null, the appropriate Student test statistics $T$ is:

$$ \frac{\rho_k}{\sqrt{\frac{1-\rho_k^2}{n-k-2}}} \sim T_{n-k-2} $$

where $n$ represents the number of observations so that $n - k$ represents the number of pairs of data $k$ apart in time.

In this study, the time series are interrupted by each change of state, and the lengths of each section are variables (Figure 4). So the number of data available to test the value of partial coefficients of correlation of order $k$ is rather $n_k$, the number of pairs of observations available at time distance $k$ accounting for the gaps in the time series. For each state-dependent variable and each state, the order of the auto-regressive process is estimated as the smallest order for which the empirical partial correlation falls within the 90% confidence intervals of an $T(n_k - 2)$.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Auto-correlation coefficients for state-dependent variables in AR1 models, plotted against resolution. The average coefficients, regardless of resolution, appear in the right marginal column.</figcaption></figure>

high resolution appears to be slightly positively biased as well (the mean of the estimation is slightly larger than the target value). Finally, for a given parameter, setting, and degradation, the ranges of fluctuations of the estimations change with the model used. This is to the disadvantage of HSMM<sub>$S, N/B$</sub> frameworks (larger inter-quantile intervals).

291

292

293

294
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Adaptation of the Student *T*-test to account for the effective number of available data $n_k$ for each tested order *k* for an auto-regressive process. Squares represent a variable of the observation layer (e.g. speed variable $V_{1,t}$) with padding indicating the state.</figcaption></figure>

# Results

## Characteristics of the simulation models

### Comparing settings

Figure 5 depicts the features that distinguish the two settings. In both settings, $d_T$ and $d_{AR}$ are co-linear, and they increase when the resolutions decrease (up to the exception of $d_T$ for Settings 2, which is between the two highest resolutions). In other words, deviations from the HMM reference situation result in proportional deviations from the AR0 reference situation, which are significantly larger for Settings 1 than for Settings 2. Fitting an HMM-AR0 model to data simulated with an HSMM$_{S,N^B}$-AR1 model is expected to be challenging, especially at high resolutions. The reverse is not a problem because AR0 and/or HMM models are special cases of AR1 and/or HSMM$_{S,N^B}$ models. The proportions of time spent in states 1 and 2 are strongly unbalanced, with opposite dominance between settings. The state-specific characteristics of the various settings are thus weighted in the opposite direction. The average sojourn time in the various states also differs significantly between the two settings. In State 1, the average duration is ten times shorter in Settings 1 than in Settings 2 (regardless of resolution). When compared to State 2, it is 1.5 times larger. Despite the similarity of the correlation coefficients, the observed difference in $d_{AR}$ can be attributed to this significant difference between the two settings. In State 1, the average correlation coefficients and integral ranges are very similar for the two settings (Figure 6, right column). We are using averages to explain the overall patterns; it is worth noting that the coefficients of correlation increase with increasing resolution. Despite similar coefficients of correlation, the distance $d_{AR1}$ for Settings 1 is 10 times larger due to their longer sojourn times. The picture for state two is different. Coefficients of correlation and integral ranges are slightly smaller in Settings 2. The slightly longer sojourn time results in a smaller $d_{AR2}$ for Settings 2. Therefore, $d_{AR1}$ and $d_{AR2}$ are larger for Settings 1. The unbalanced states' proportions between settings resulted in a larger $d_{AR}$ for Settings 1 but reduced the difference.

Finally, given the overlaps between the PDFs of the speed variables by states, and the expected values of the different sojourn times, the expected state decoding capacities when using mixture models ($d_V$) are larger for Settings 2 than for Settings 1 (Figure 5), with 95% of good predictions for the largest resolution for Settings 2 as opposed to 80% for Settings 1 (Figure 9). All of these statements point to disappointing decoding accuracy for Settings 1 and decreasing resolutions in general. This latter statement suffers one exception: for Settings 1; the accuracy appears slightly better for the lowest resolution than the penultimate resolution (Figure 9). The expected values of the PMFs are equal in HMM and HSMM$_{S,N^B}$ formulations, but the standard deviations differ. For Settings 1 and 2, the standard deviations are co-linear between the two types of Markov
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The decimal logarithm of standard deviations of the sojourn time PMFs for the various settings. The x-axis represents the time duration of HSMM cases with shifted negative binomial PMFs. The y-axis refers to HMM cases with geometric PMFs.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Graphical representation of the model settings. Settings 1 and 2 are summarized by $d_T$ the Hellinger distances between the shifted-negative binomial (HSMM) and the geometric (HMM) PMF formulations for the sojourn time on the x-axis, by $d_{AR}$ the difference between AR1 and AR0 formulations on the y-axis, and by $d_V^{mod}$ the state decoding accuracies in case of mixture model through the surface area of the circles. The proportions of time spent in each state are also represented by two half discs within the circle. The arrows in each setting indicate the path towards decreasing resolution.</figcaption></figure>

formulations but are generally larger in the HMMs (Figure 7). It also appears that, for Settings 1, the standard deviations of state durations in States 1 and 2 are of the same order of magnitude, whereas for Settings 2, they are one order of magnitude larger in State 1 than in State 2.

MSA (estimation model = simulation model)

Parameter estimation

In general, parameter estimates under MSA conditions are unbiased and symmetrically distributed around the target value (Figure 8). One notable exception is the variance of the Gaussian PDF for $V_{1,2}$ for Settings 2 when using an HSMM$_{s,NB}$-AR0 model. When the target value is 1.26, 50% of the estimations fall between $Q_{25\%} = 1.25$ and $Q_{75\%} = 1.44$, indicating a minor issue. The variance of the shifted negative binomial at high resolution appears to be slightly positively biased as well (the mean of the estimation is slightly larger than the target value). Finally, for a given parameter, setting, and degradation, the ranges of fluctuations of the estimations change with the model used. This is to the disadvantage of HSMM$_{s,NB}$ frameworks (larger inter-quantile intervals).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Simulation-estimation experiments: box-plots of the 100 parameter estimation when the model used for the estimation is the same as the model used for the simulation (MSA experiments with 250 states alternations). Each panel represents a settings × a degradation. The sequence of the box-plots is always the same: HMM-AR0, HMM-AR1, HSMM-AR0, HSMM-AR1. For each parameter, the y-axis is set to be consistent across all four resolutions for given settings. A horizontal line represents the simulation model parameter's (reference) value. The shift for the sNB PMF is not represented.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Auto-correlation coefficients for state-dependent variables in AR1 models, plotted against resolution. The average coefficients, regardless of resolution, appear in the right marginal column.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** The decimal logarithm of standard deviations of the sojourn time PMFs for the various settings. The x-axis represents the time duration of HSMM cases with shifted negative binomial PMFs. The y-axis refers to HMM cases with geometric PMFs.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Accuracies (y-axis) of state estimations as a function of the resolution (x-axis expressed as the average number of time steps in State 2). **Simulation-estimation experiments (with 250 states alternations)**: in a row, the simulation models, and in the column, the estimation models. The light-colored envelopes are built on the basis of the 25% and 75% quantiles of the accuracies obtained over the 100 simulations performed for Settings 1 (blue) and 2 (red). The white envelopes correspond to the 10% and 90% quantiles. The white diagonal panels correspond to model skill assessment - MSA (estimation model equals the simulation model). The grey panels correspond to the model robustness assessment - MRA (estimation model equals all but the simulation model). **Model characteristics**: the continuous lines indicate the expected accuracies when considering a mixed model ($d_V$). **Real case accuracies**: the dashed lines in the bottom right panel correspond to the accuracies for real data estimations when using an HSMM<sub>sNB</sub>-AR1.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Simulation-estimation experiments: box-plots of the 100 parameter estimation when the model used for the estimation is the same as the model used for the simulation (MSA experiments with 250 states alternations). Each panel represents a settings × a degradation. The sequence of the box-plots is always the same: HMM-AR0, HMM-AR1, HSMM-AR0, HSMM-AR1. For each parameter, the y-axis is set to be consistent across all four resolutions for given settings. A horizontal line represents the simulation model parameter's (reference) value. The shift for the sNB PMF is not represented.</figcaption></figure>
## State decoding accuracies

Figure 9 (diagonal panels) shows that state decoding accuracy is higher in Settings 2 compared to Settings 1, regardless of the model. In both cases, increasing the resolution (or the number of observations per state) improves accuracy. There is no difference between the model types. The state decoding accuracy is similar for the most complex models (HSMM-AR1) to the simplest ones (HMM-AR0), which outperform by around 10% what would be expected under a more simple mixed model framework. The AR1 process, when present in the data and model structure, serves as a source of improvement. It's worth noting that the results are consistent with those obtained using mixed models. For Settings 1, performance improves when moving from the smallest to the penultimate smallest resolution. Across the 100 simulations, the accuracy values are very similar. This is obvious for Settings 2 at all resolutions, demonstrating no difference between simulations. For Settings 1, some variations occur at low resolutions but remain small, with an average accuracy equal to 80%. These results serve as the baseline for state decoding performance because they are obtained under the most favorable conditions, with the estimation model equal to the simulation model. It is thus expected that performance will suffer in less favorable situations.

## MRA (estimation model = all but the simulation model)

In MRA situations (non-diagonal panels of Figure 9 taken in a row), the overall picture remains the same: the accuracies are systematically larger for Settings 2 and get larger when the resolution increases. The largest changes in performances appear when the simulation model is AR1 and the estimation model is AR0 (Figure 9, panels of lines 2 and 4). The range of fluctuations of the accuracies remains similar to MSA conditions (very small for Settings 2 and small for Settings 1). As a result, average values in each numerical experiment provide a good summary of accuracy. The average relative accuracy losses (Figure 10) are shown below. The overall average loss is 1.2%. However, this average loss conceals contrasting situations, as follows.

Using the correct Markov structure (panels of the middle line), estimation with an AR1 has no impact (average loss equal to 0), which is consistent with the fact that an AR1 can capture an AR0 as a special case (AR0 ⊂ AR1). However, in the reverse situation (AR1 => AR0), accuracies deteriorate on average and significantly more for Settings 1. The loss of accuracy is linearly related to the distance between AR1 and AR0 formulations: the larger the distance $d_{AR}$, the greater the loss of accuracy (Figure 11).

Using the correct AR structure (panels of the middle column), an estimation with a HSMM$_{S,N^B}$ when the data is generated with an HMM has little impact. Accuracies are, however, slightly tampered with for Settings 2, whose $d_T$ values are nevertheless smaller than those for Settings 1 (Figure 5). These results support the idea that an HSMM$_{S,N^B}$ can capture an HMM as a special case (HMM ⊂ HSMM$_{S,N^B}$). The situation remains the same the other way round (HSMM$_{S,N^B}$ => HMM) contrary to expectations. One particular case of Settings 1 is particularly harmful, while the rest have no effect.

Estimating with an AR1 structure when the data is AR0 does not affect the state decoding accuracy (left column). The most significant impacts occur when estimating using an AR0 structure while the true process is an AR1 (right column panels), and the value case is when also forcing the estimation to be HMM while the data are HSMM$_{S,N^B}$ (average loss equals 3.6%). Changing from HSMM$_{S,N^B}$ to HMM is only a concern when there is a negative change in the AR structure of the estimating model (i.e. AR0 instead of AR1). In this worst case, a fine inspection of the parameters' re-estimations when simulating with an HSMM-AR1 and estimating with one of the three other models (Figure 12) shows that nearly all parameters are impacted by this bi-layers discrepancy between the simulating and estimating model. What is also clear is that the main impact occurs when the state-dependent variables are modeled using the AR0 structure. For example, when an HMM-AR0 or an HSMM-AR0 is used to estimate the variance parameters of the Gaussian PDFs of $V_1$, there is a significant
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Accuracies (y-axis) of state estimations as a function of the resolution (x-axis expressed as the average number of time steps in State 2). **Simulation-estimation experiments (with 250 states alternations)**: in a row, the simulation models, and in the column, the estimation models. The light-colored envelopes are built on the basis of the 25% and 75% quantiles of the accuracies obtained over the 100 simulations performed for Settings 1 (blue) and 2 (red). The extra transparent envelopes correspond to the 10% and 90% quantiles. The white diagonal panels correspond to model skill assessment - MSA (estimation model equals the simulation model). The grey panels correspond to the model robustness assessment - MRA (estimation model equals all but the simulation model). **Model characteristics**: the thick and coloured continuous lines indicate the expected accuracies when considering a mixture model ($d_V$).</figcaption></figure>
bias. Expected values of sojourn time are also significantly affected. However, the most significant effects are
observed when an HMM formulation is used for the estimation model.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Model robustness assessment-MRA (with 250 states alternations). The percentage of state estimation accuracy lost when using a different estimation model than the one used to simulate the data. The boxes show the ranges of fluctuations. Points are jittered for better readability. The black dots on the right of each panel represent the overall mean of each case (blue and red for Settings 1 and 2, respectively). The changes in model structures between simulation and estimation are shown in the left and top margins. Columns represent a possible change in the auto-regressive characteristics of the speed distributions, e.g., “AR0 => AR1” means simulating with AR0 and re-estimating with AR1. Rows correspond to a possible change of the Markov property of the model, e.g., “HMM => HSMM” means: simulating with HMM and re-estimating with HSMM<sub>sNB</sub>. The central panel corresponding to the MSA situation is not considered.</figcaption></figure>
## State decoding accuracies

Figure 9 (diagonal panels) shows that state decoding accuracy is higher in Settings 2 compared to Settings 1, regardless of the model. In both cases, increasing the resolution (or the number of observations per state) improves accuracy. There is no difference between the model types. The state decoding accuracy is similar for the most complex models (HSMM-AR1) to the simplest ones (HMM-AR0), which outperform by around 10% what would be expected under a more simple mixture model framework. The AR1 process, when present in the data and model structure, serves as a source of improvement. It's worth noting that the results are consistent with those obtained using mixture models. For Settings 1, performance improves when moving from the smallest to the penultimate smallest resolution. Across the 100 simulations, the accuracy values are very similar. This is obvious for Settings 2 at all resolutions, demonstrating no difference between simulations. For Settings 1, some variations occur at low resolutions but remain small, with an average accuracy equal to 80%. These results serve as the baseline for state decoding performance because they are obtained under the most favorable conditions, with the estimation model equal to the simulation model. It is thus expected that performance will suffer in less favorable situations.

## MRA (estimation model = all but the simulation model)

In MRA situations (non-diagonal panels of Figure 9 taken in a row), the overall picture remains the same: the accuracies are systematically larger for Settings 2 and get larger when the resolution increases. The largest changes in performances appear when the simulation model is AR1 and the estimation model is AR0 (Figure 9, panels of lines 2 and 4). The range of fluctuations of the accuracies remains similar to MSA conditions (very small for Settings 2 and small for Settings 1). As a result, average values in each numerical experiment provide a good summary of accuracy. The average relative accuracy losses (Figure 10) are shown below. The overall average loss is 1.2%. However, this average loss conceals contrasting situations, as follows.

Using the correct Markov structure (panels of the middle line), estimation with an AR1 has no impact (average loss equal to 0), which is consistent with the fact that an AR1 can capture an AR0 as a special case (AR0 ⊂ AR1). However, in the reverse situation (AR1 => AR0), accuracies deteriorate on average and significantly more for Settings 1. The loss of accuracy is linearly related to the distance between AR1 and AR0 formulations: the larger the distance $d_{AR}$, the greater the loss of accuracy (Figure 11).

Using the correct AR structure (panels of the middle column), an estimation with a HSMM$_{SVB}$ when the data is generated with an HMM has little impact. Accuracies are, however, slightly tampered for Settings 2, whose $d_T$ values are nevertheless smaller than those for Settings 1 (Figure 5). A geometric PMF is a particular case of a shifted negative binomial PMF ($G(q) = \text{SVB}(\text{shift} = 1, n = 1, q)$). However, the PMF for the (shifted) negative binomial PMF gets three parameters when the geometric PMF gets only one parameter. This certainly hampers inferences. In practice, a HSMM can tend to an HMM but might not reach it totally. Inferences might also be thwarted by the fact that the *n* parameter of the negative binomial PMF is not an integer, but a real. So it can tend to 1 without being equal 1. Looking at the *n* parameter values obtained when fitting a HSMM while the data are HMM (in supplementary material) allows evaluating the capacity of the EM algorithm to end up with *n* = 1 (Figure 15). It happens that *n* does fluctuate around 1 with alternate behavior (above/under) given the state. The state that lasts the shortest, i.e. state 1 for settings 1 and state 2 for settings 2, is the one with *n* smaller than one (and vice versa). The fluctuations in *n* increases when the data resolution increases, with larger fluctuations for settings 1 where $q_{75\%} = 2$ and where *n* sometimes equals 5. Finally, there is no difference in the output when using ARO or AR1 structures for the speed variables. These results support the idea that an HSMM$_{SVB}$ can capture an HMM as a special case (HMM ⊂ HSMM$_{SVB}$). The most interesting output is that the situation remains the same the other way round (HSMM$_{SVB}$ => HMM) contrary to expectations. One particular case of Settings 1 is particularly harmful, while the rest have no effect.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Simulation-estimation analyses-MRA (experiments with 250 states alternations). Loss of accuracy of state estimations when simulating with an AR1 model vs. estimating with an AR0 model as a function of the distance between AR1 and AR0 formulations for Settings 1 (blue) and 2 (red).</figcaption></figure>
Estimating with an AR1 structure when the data is AR0 does not affect the state decoding accuracy (left column). The most significant impacts occur when estimating using an AR0 structure while the true process is an AR1 (right column panels), and the value case is when also forcing the estimation to be HMM while the data are HSMM<sub>SNB</sub> (average loss equals 3.6%). Changing from HSMM<sub>SNB</sub> to HMM is only a concern when there is a negative change in the AR structure of the estimating model (i.e. AR0 instead of AR1). In this worst case, a fine inspection of the parameters’ re-estimations when simulating with an HSMM-AR1 and estimating with one of the three other models (Figure 12) shows that some but not all parameters are impacted by this bi-layers discrepancy between the simulating and estimating model. Parameters corresponding to the first moments of the speed PDFs (e.g. μ or τ) are not significantly impacted in the sense that their re-estimations are similar to those obtained in MSA conditions. What is also clear is that the main impact occurs when the state-dependent variables are modeled using the AR0 structure. For example, when an HMM-AR0 or an HSMM-AR0 is used to estimate the variance parameters of the Gaussian PDFs of V<sub>1</sub>, there is a significant bias. Expected values of sojourn time are also significantly affected. However, the most significant effects are observed when an HMM formulation is used for the estimation model.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Simulation-estimation experiments: box-plots of the 100 parameter estimation when the model used for the estimation all but the HSMM-AR1 model used for the simulation (MRA experiments with 250 states alternations). Each panel represents a settings × a degradation. The sequence of the box-plots is always the same: HMM-AR0, HMM-AR1, HSMM-AR0, HSMM-AR1. The (reference) value of the parameter of the simulation model is represented by a horizontal line. The shift for the sNB PMF is not represented.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Model robustness assessment-MRA (with 250 states alternations). The percentage of state estimation accuracy lost when using a different estimation model than the one used to simulate the data. The boxes show the ranges of fluctuations. Points are jittered for better readability. The black dots on the right of each panel represent the overall mean of each case (blue and red for Settings 1 and 2, respectively). The changes in model structures between simulation and estimation are shown in the left and top margins. Columns represent a possible change in the auto-regressive characteristics of the speed distributions, e.g., "AR0 => AR1" means simulating with AR0 and re-estimating with AR1. Rows correspond to a possible change of the Markov property of the model, e.g., "HMM => HSMM" means: simulating with HMM and re-estimating with HSMM<sub>SN/3</sub>. The central panel corresponding to the MSA situation is not considered.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Data characteristics for Vessel 1 (blue) and Vessel 2 (red). State 1 (non-fishing) is in open symbols, and State 2 (fishing) is in plain symbols.</figcaption></figure>

**Top panels:** real trajectories with centered coordinates for data resolution of 1 h. The histograms of the sojourn time are provided below.

**Middle panels:** correlation plots of the persistent ($V_p$) and the rotation ($V_r$) speeds. The marginal histograms represent the univariate distribution per state.

**Bottom panel:** summary metrics for the two vessels as a function of the data resolutions. The Hellinger distances between the sojourn time histograms and a geometric PMF ("G") or a shifted negative binomial PMF ("sNB"), and the accuracies expected under mixed model ($d_V$).
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Simulation-estimation analyses-MRA (experiments with 250 states alternations) where simulations are performed with AR1 models and estimations with AR0 models. Loss of accuracy of state estimations as a function of the distance between AR1 and AR0 formulations *d*<sub>AR</sub> for Settings 1 (blue) and 2 (red).</figcaption></figure>
# Real cases

343

## Choosing a model structure and describing its main characteristics

344

Figure 13 shows the general characteristics of two sets of real tracks. As the sequences of true behavioral states are known, it is possible to select the best model by evaluating parametric assumptions for both the hidden state variables and the state-dependent variables. However, given a large number of observations (from several hundred to a few thousand depending on resolution), any minor difference between model and data characteristics results in low p-values for null hypothesis testing, even at low-risk levels ("large N problem").

345

346

347

348

349

350

## Hidden layer

351

For Vessel 1, the Hellinger distances between the histograms of sojourn time duration in each state and either of the two possible PMFs to model the state process (i.e. geometric and shifted negative binomial) indicate that the most appropriate model is always an HSMM<sub>S,N̄B</sub> (Figure 13). For Vessel 2, this could be less relevant because the distances are greater. But choosing a geometric PMF would even be worth it.

352

353

354

355

## Observation layer

356

The state-dependent variables exhibit temporal auto-correlation at lag 1 (Figure 6). Among the four structures tested in this work, the most appropriate model structure is an HSMM<sub>S,N̄B</sub> - AR1, which is also the most flexible one with regards to inference objectives as it can converge down to an HMM - AR0 if necessary. However, in the framework of auto-regressive processes, which is an assumption itself, the orders of the AR processes of the speed variables are likely larger than 1 (Figure 14). This is very clear for Vessel 1, which has auto-regressive processes with an order greater than 1 in half of the cases (9/16), as opposed to 5/16 cases for Vessel 2, and has a larger partial auto-correlation than Vessel 2. The correlation plots show that the persistent and rotation speeds are not independent given the state (except for Vessel 2 in State 2, Figure 13) which contradicts the model's assumptions. Finally, Vessel 1 has strong mixing (respectively small $d_V$) between the empirical histograms of the speed in different states, particularly for the rotation speed (Figure 13), indicating poor state decoding performance.

357

358

359

360

361

362

363

364

365

366

367

## State decoding performances

368

HSMM<sub>S,N̄B</sub> - AR1 models and speed-based filters are applied to the speeds time series available for the two vessels at the four different time resolutions (Figures 3) to estimate the (known but removed) states. Given the number of fishing operations (corresponding to State 2 by convention here) observed along each trajectory (N = 276 for Vessel 1 and N = 31 for Vessel 2; Figure 13), real state decoding accuracy is compared to the simulation-estimation experiments based on 250 and 50 alternations for Settings 1 and 2 respectively. Real-case estimates (Figure 9) perform worse than simulation-estimation experiments. This is especially true for Vessel 1 at high resolution where the real-case state decoding performances clearly depart from the simulation-estimation experiment (70% of good state estimation for the HSMM<sub>S,N̄B</sub> - AR1 model, compared to 90% for the simulation-estimation experiment). Although not shown in this paper, these models surpass the conventional filtering approach based on a speed threshold to routinely identify behavior for fisheries management.

369

370

371

372

373

374

375

376

377

378

379
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Simulation-estimation experiments: box-plots of the 100 parameter estimation when the model used for the estimation equals all but the HSMM-AR1 model used for the simulation (MRA experiments with 250 states alternations). Each panel represents a settings × a degradation. The sequence of the box-plots is always the same: HMM-AR0, HMM-AR1, HSMM-AR0, HSMM-AR1 (the first three in color correspond to MRA conditions; the last one in white corresponds to MSA situations). The (reference) value of the parameter of the simulation model is represented by a horizontal line. The shift for the SNB PMF is not represented.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Partial auto-correlations of state-dependent variables. For each variable, vessel, resolution, and state, the partial auto-correlation coefficients are represented up to the point when they are no longer statistically different from 0. The inner panel summarizes the order obtained in each case with a reference being order 1.</figcaption></figure>

# Discussion

## Scope of the benchmarking

The current study's benchmark includes four model types (some with nesting properties) and two sets of contrasted settings. We chose to focus on these restricted configurations rather than following an experimental design crossing parameter values over a wider range. such a comprehensive experimental design would have required a significant amount of computing time. Instead, we defined the parameters of the four types of models using two sets of real fishing vessel trajectories that provided ground truth on the sequence of states. While this likely limits the overall generalizability of the results, it allows for useful recommendations for various real-world situations involving other natural foragers than fishers. According to Nettle et al., 2013, the growing interest in human behavioral ecology presents a significant opportunity to bridge the gap between natural and social sciences. Regarding fisheries science, Bertrand et al., 2007 show that fishermen follow similar patterns as natural top predators. In this context, specifying the models using annotated vessel trajectories places the study in a broader, more relevant ecological context. Without being too general, the scope of this work is quite broad. The model parameters are sufficiently contrasted to distinguish first-order

380

381

382

383

384

385

386

387

388

389

390

391

392

393
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Data characteristics for Vessel 1 (blue) and Vessel 2 (red). State 1 (non-fishing) is in open symbols, and State 2 (fishing) is in plain symbols.</figcaption></figure>

**Top panels:** real trajectories with centered coordinates for data resolution of 1 h. The histograms of the sojourn time are provided below.

**Middle panels:** correlation plots of the persistent ($V_p$) and the rotation ($V_r$) speeds. The marginal histograms represent the univariate distribution per state.

**Bottom left panel:** summary metrics for the two vessels as a function of the data resolutions:

* G: The Hellinger distances between the sojourn time histograms and a geometric PMF

* sNB: The Hellinger distances between the sojourn time histograms and a shifted negative binomial PMF

* d: the discrimination coefficient between speed histograms ($d_{V}^{emp}$)

**Bottom right panel:** state decoding accuracies when using an HSMM<sub>SNB</sub>-AR1 (dashed lines). The light-colored and white envelopes are alike those presented in figure 9. They correspond to the quantiles of the accuracies of the simulation-estimation experiment with an HSMM<sub>SNB</sub>-AR1.
auto-regressive processes from their AR0 counterparts and HSMMs from their HMM counterparts. So, the model parameters cause a significant contrast between the four model types.

The four model structures investigated in the experiments are reasonably simple (two hidden states with two independent observed variables). Most models in movement ecology (including fisheries) adopt simple modeling frameworks (e.g. Bez et al., 2011; Langrock et al., 2012; Vermard et al., 2010) that are still recommended in more recent studies (Auger-Méthé, Newman, et al., 2021; McClintock, 2021). The number of states and the identification of a specific state among them are sometimes the primary goals of using state space models, particularly in behavioral ecology. The goal of this analysis is not to infer the possibility of correctly estimating the number of states but rather to quantify the state decoding accuracy of hidden Markov models under controlled conditions. If you are interested in this issue, we recommend reading the pragmatical approach proposed by Pohle et al., 2017. Others sources of possible biases coming from temporal misalignment between model and observations are also not considered (Glennie et al., 2023). Although not shown, two intriguing side outputs are worth noting. First, the Estimation-Maximization (EM) algorithm produced reliable estimates of maximum likelihood, i.e., not impacted by changes in initialization values. Second, the most probable state sequences (Viterbi) were equal to marginal estimates, implying that the Viterbi algorithm is unnecessary in these cases.

## No best model in model skill assessment conditions

The best model performances are obtained from model skill assessments, which occur when the simulation and estimation models are identical. This is consistent with the fact that, given a fixed number of observations, one cannot expect to achieve higher accuracy when the data and model are no longer consistent. Given favorable inference situations (i.e., 250 state alternations), all four model types, including the most complex (e.g., HSMM-AR1), achieve similar state decoding accuracies and very good parameter estimations. This is not an exceptional situation; one of our two case studies shows the same amount of state alternation. This means that there is no optimal model type under MSA conditions. Model parameters and settings (including resolution) control performance rather than the model structure itself. The following examples illustrate this statement. In both settings, the variances of the sojourn times are greater in HMM formulation than in HSMM formulation (Figure 7).

The only distinction between the two is that, for Settings 2, the variance of the sojourn time is an order of magnitude greater in State 2 than in State 1, while the proportion of time spent in State 2 is quite small ($\pi_2$=0.15, $\pi_1$=0.85). This could have been detrimental to the fittings with Settings 2. However, we observe the opposite. The difference in performances is driven by the two metrics ($d_V$ and $d_{AR}$) that are clearly in favor of Settings 2 (Figure 5) and that explain strongly the shapes and levels of the different accuracies (Figure 9). Regardless of model type, the characteristics of the observation layer (more favorable for Settings 2 than for Settings 1) explain the performance better than those of the hidden layer. Under MSA conditions, state-space models consistently outperform mixed models (Figure 9). This finding suggests that the structures introduced by temporal state processes facilitate state estimation, at least as long as the estimation model agrees with the simulation model. Consider applying this result to real-world situations: in that case, the ability to select the true temporal process will have a greater impact on state inference quality than the change state assumption.

## Auto-correlation deteriorates the state decoding accuracy

The numerical experiments revealed that the deterioration in accuracy is primarily due to auto-correlation in the observation layers rather than the hidden layer's Markovian structure. Although $d_T$ (hidden layer) and $d_{AR}$ (observation layer) are co-linear in the settings with potential confounding effects, the current work demonstrates that the modeling choices related to the observation layer are of primary importance for state
# Real cases

## Choosing a model structure and describing its main characteristics

Figure 13 shows the general characteristics of two sets of real tracks. As the sequences of true behavioral states are known, it is possible to select the best model by evaluating parametric assumptions for both the hidden state variables and the state-dependent variables. However, given a large number of observations (from several hundred to a few thousand depending on resolution), any minor difference between model and data characteristics results in low p-values for null hypothesis testing, even at low-risk levels ("large N problem").

## Hidden layer : testing the Markov and semi-Markov hypotheses

For Vessel 1, the Hellinger distances between the histograms of sojourn time duration in each state and either of the two possible PMFs to model the state process (i.e. geometric and shifted negative binomial) indicate that the most appropriate model is always an HSMM<sub>SNB</sub> (Figure 13). For Vessel 2, this could be less relevant because the distances are greater. But choosing a geometric PMF would even be worth it.

## Observation layer : testing the AR hypothesis

The state-dependent variables exhibit temporal auto-correlation at lag 1 (Figure 6). Among the four structures tested in this work, the most appropriate model structure is an HSMM<sub>SNB</sub> - AR1, which is also the most flexible one with regards to inference objectives as it can converge down to an HMM - AR0 if necessary. However, in the framework of auto-regressive processes, which is an assumption itself, the orders of the AR processes of the speed variables are likely larger than 1 (Figure 14). This is very clear for Vessel 1, which has auto-regressive processes with an order greater than 1 in half of the cases (9/16), as opposed to 5/16 cases for Vessel 2, and has a larger partial auto-correlation than Vessel 2. The correlation plots show that the persistent and rotation speeds are not independent given the state (except for Vessel 2 in State 2, Figure 13) which contradicts the model's assumptions. Finally, Vessel 1 has strong mixing (respectively small $d_V$) between the empirical histograms of the speed in different states, particularly for the rotation speed (Figure 13), indicating poor state decoding performance.

## State decoding performances

HSMM<sub>SNB</sub> - AR1 models and speed-based filters are applied to the speeds time series available for the two vessels at the four different time resolutions (Figures 3) to estimate the (known but removed) states. Given the number of fishing operations (corresponding to State 2 by convention here) observed along each trajectory (N = 276 for Vessel 1 and N = 31 for Vessel 2; Figure 13), real state decoding accuracy is compared to the simulation-estimation experiments based on 250 and 50 alternations for Settings 1 and 2 respectively. Real-case estimates (Figure 13) perform worse than simulation-estimation experiments. This is especially true for Vessel 1 at high resolution where the real-case state decoding performances clearly depart from the simulation-estimation experiment (70% of good state estimation for the HSMM<sub>SNB</sub> - AR1 model, compared to 90% for the simulation-estimation experiment). Although not shown in this paper, these models surpass the conventional filtering approach based on a speed threshold to routinely identify behavior for fisheries management.
decoding performances. This finding is supported by both simulation-estimation experiments and real-data estimation.

Simulation-estimation experiments

The HSMM<sub>S<sup>NB</sup></sub> - AR1 model structure is the most robust to discrepancies between model hypotheses and data characteristics when retrieving hidden states, according to model robustness assessment (MRA). While the accuracy in state estimation deteriorates in 90% of cases when using a model structure different from the one used to simulate data, the relative loss is small (0.92% on average). Imposing no auto-correlation in the estimating model while allowing some for state-dependent variables has a negative impact on state decoding performance. The loss of accuracy is well explained by *d*<sub>AR</sub>, which is the distance between AR1 and AR0 formulations. This is not true for the hidden layer, where imposing an HMM structure while the state process is semi-Markov has no negative impact on the state decoding performance. This finding has direct practical implications. It demonstrates that it is preferable to concentrate on the hypothesis choices of the observation layer rather than the hidden layer. In some ways, this is reassuring, as only the observation layer contains empirical data. Even though, in theory, the time dependence of the observed processes or the distinction of the PDFs should be evaluated conditionally on the hidden states (which is impossible in practice because the states are, by definition, hidden), we recommend conducting exploratory analyses of the observed process data to define auto-correlation hypotheses.

Real cases

The analysis of model performance on real cases is consistent with the simulation results, indicating that discrepancies between the properties of the observations and the model's hidden layer assumptions do not affect performance. Both data sets demonstrate a visual deviation from the various model assumptions. However, even though the differences in residence time for Vessel 2 are greater, the model's performance improves. The most robust model (HSMM<sub>S<sup>NB</sup></sub> - AR1) results in lower state decoding performance for Vessel 1 and higher performance for Vessel 2. Interestingly, Vessel 2's performances do not deteriorate strongly with the resolution, whereas Vessel 1's do (Figure 9). The two primary differences in favor of Vessel 2 concern the speed variables:

* The discriminating power between the speed distributions measured through the integrated metric $d_V$ is low, respectively large, for Vessel 1 and Vessel 2 (Figure 13). This contributes to the low accuracies obtained for Vessel 1.

* In most cases, the partial autocorrelation of observed speed is greater than one (the model's assumed autocorrelation), particularly for Vessel 1. Not only are the auto-correlations of order 1 larger, but the most likely order of the auto-regressive process modeling the speed variable is greater than 1.

Beyond these considerations about the value of coefficients of auto-correlation and the orders of the auto-regressive processes, the observations may not be consistent with an auto-regressive process at all. In particular, auto-regressive processes assume that the relationship between successive observations is linear, which may not be the case. As an example, we provide the cross-plots between successive speed components (Additional Figure 15), which violate linear relationships over time. To form the correct hypotheses about the observation process, start by ensuring that the data are consistent with an auto-regressive process. The violation of this assumption corresponds to a structural error process that is rarely evaluated. Before discussing the order of the AR, consider the more fundamental choice of using an AR formulation for the observation layer (Gloaguen et al., 2015; Joo et al., 2013; Morales et al., 2004). Second, within the context of auto-regressive processes, the data may not correspond to a zero or first AR order (AR0 or AR1). This second type of violation is caused by mis-specification and can be more easily controlled by evaluating the model's goodness of fit.
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Partial auto-correlations of state-dependent variables. For each variable, vessel, resolution, and state, the partial auto-correlation coefficients are represented up to the point when they are no longer statistically different from 0. The inner panel summarizes the order obtained in each case with a reference being order 1.</figcaption></figure>

# Discussion

## Scope of the benchmarking

The current study's benchmark includes four model types (some with nesting properties) and two sets of contrasted settings. We chose to focus on these restricted configurations rather than following an experimental design crossing parameter values over a wider range. Such a comprehensive experimental design would have required a significant amount of computing time. Instead, we defined the parameters of the four types of models using two sets of real fishing vessel trajectories that provided ground truth on the sequence of states. While this likely limits the overall generalizability of the results, it allows for useful recommendations for various real-world situations involving other natural foragers than fishers. According to Nettle et al., 2013, the growing interest in human behavioral ecology presents a significant opportunity to bridge the gap between natural and social sciences. Regarding fisheries science, Bertrand et al., 2007 show that fishermen follow similar patterns as natural top predators. In this context, specifying the models using annotated vessel trajectories places the study in a broader, more relevant ecological context. Without being too general, the scope of this work is quite broad. The model parameters are sufficiently contrasted to distinguish first-order
Increasing the resolution: a tricky solution?

As expected, model characteristics are linked to data resolution. This affects not only the averages of the sojourn times, which are directly impacted because they are expressed in a series of time steps but also their distributions. This also affects the variability of the speeds (both persistent and rotational), their time dependencies, and the overlap of speed PDFs between states. The multidimensional nature of data resolution precludes statistical investigation of its effect on estimation performance, all other things being equal. Several authors (e.g., McClintock, King, et al., 2012; Vermard et al., 2010) have pointed out that time steps must be chosen to match the scale at which behavioral decisions are made using expert knowledge. Katara and Silva Katara and Silva, 2017 examined the mismatch between the temporal resolution of position records and the timescale of fishing activity. The authors investigated the effect of time resolution (from 10 minutes to 2 hours) on bias and errors in identifying fishing operations for vessels (Portuguese purse seiners) with characteristics similar to those of vessel 2. The results (7% of missed fishing sets with 2h time interval compared to 10 minutes) align with the estimated increase in accuracy with time scale resolution for Vessel 2.

Similar results were obtained for other vessels, emphasizing the importance of high-resolution time data Edelhoff et al., 2016; Mendo et al., 2019. Interestingly, when using data that were consistent with a given model's assumptions, such as in simulation-estimation experiments, increasing the resolution of simulated data resulted in higher accuracy. This was less obvious for real data, where the results differed between the two case studies and were not linear with data resolution. There is a contradiction between the need for high-frequency data to track short-term states or to get as close as possible to the tracked individual's actual behavior and the ability to model this behavior using a random decision. If we suppose that memory and/or inertia play a role in behavioral changes at high frequencies, in this case, it is clear that a state is dependent not only on previous states but also on multiple previous ones. As a result, a compromise must be made between the pursuit of better statistical conditions and a realism threshold to achieve an optimal level of data resolution.

In conclusion, in silico, when data characteristics are controlled by and consistent with model characteristics, simulation-estimation experiments show that a change in Markov property for the estimation model has no negative impact on state decoding performance unless one forces the state-dependent variable to be AR0 when they are indeed AR1. This aligns with the fact that $HSMM_{s,N^B}$ and AR1 formulations are generalizations of the usual HMM and AR0 formulations. This study also shows that model performance is more sensitive to the auto-regressive characteristics of the state-dependent variables than to the Markov properties of the hidden states. This is somewhat satisfactory because, in practice, one may have access to some aspects of the observation layer but not to the hidden, unobserved layer. This allows for empirical controls on what appear to be key model assumptions in terms of state decoding performance. However, in real cases, the robustness of $HSMM_{s,N^B}$ - AR1 appears low, with instances where its state decoding performances are worse than simple filtering. The confusion between the distributions of state-dependent variables in different states, whether observable (usually not possible) or more realistically speculated, is the most important factor influencing state decoding performance in real-world scenarios. Blind applications of state-space models in the absence of prior knowledge about the processes lead to radical errors and poor state-decoding performance. Given that HMMs are a popular choice in ecology due to their ease of use, we believe there is a need for serious consideration of the consequences of the numerous potential deviations between data and model assumptions.

Appendices

Supplementary figure
auto-regressive processes from their AR0 counterparts and HSMMs from their HMM counterparts. So, the model parameters cause a significant contrast between the four model types.

The four model structures investigated in the experiments are reasonably simple (two hidden states with two independent observed variables). Most models in movement ecology (including fisheries) adopt simple modeling frameworks (e.g. Bez et al., 2011; Langrock et al., 2012; Vermard et al., 2010) that are still recom- mended in more recent studies (Auger-Méthé, Newman, et al., 2021; McClintock, 2021). The number of states and the identification of a specific state among them are sometimes the primary goals of using state space models, particularly in behavioral ecology. The goal of this analysis is not to infer the possibility of correctly estimating the number of states but rather to quantify the state decoding accuracy of hidden Markov mod- els under controlled conditions. If you are interested in this issue, we recommend reading the pragmatical approach proposed by Pohle et al., 2017. Others sources of possible biases coming from temporal misalign- ment between model and observations are also not considered (Glennie et al., 2023). Although not shown, two intriguing side outputs are worth noting. First, the Estimation-Maximization (EM) algorithm produced reli- able estimates of maximum likelihood, i.e., not impacted by changes in initialization values. Second, the most probable state sequences (Viterbi) were equal to marginal estimates, implying that the Viterbi algorithm is unnecessary in these cases.

No best model in model skill assessment conditions

The best model performances are obtained from model skill assessments, which occur when the simulation and estimation models are identical. This is consistent with the fact that, given a fixed number of observations, one cannot expect to achieve higher accuracy when the data and model are no longer consistent. Given fa- vorable inference situations (i.e., 250 state alternations), all four model types, including the most complex (e.g., HSMM-AR1), achieve similar state decoding accuracies and very good parameter estimations. This is not an exceptional situation; one of our two case studies shows the same amount of state alternation. This means that there is no optimal model type under MSA conditions. Model parameters and settings (including resolution) control performance rather than the model structure itself. The following examples illustrate this statement. In both settings, the variances of the sojourn times are greater in HMM formulation than in HSMM formulation (Figure 7).

The only distinction between the two is that, for Settings 2, the variance of the sojourn time is an order of magnitude greater in State 2 than in State 1, while the proportion of time spent in State 2 is quite small ($\pi_2$=0.15, $\pi_1$=0.85). This could have been detrimental to the fittings with Settings 2. However, we observe the opposite. The difference in performances is driven by the two metrics ($d_V$ and $d_{AR}$) that are clearly in favor of Settings 2 (Figure 5) and that explain strongly the shapes and levels of the different accuracies (Figure 9). Regardless of model type, the characteristics of the observation layer (more favorable for Settings 2 than for Settings 1) explain the performance better than those of the hidden layer. Under MSA conditions, state-space models consistently outperform mixture models (Figure 9). This finding suggests that the structures introduced by temporal state processes facilitate state estimation, at least as long as the estimation model agrees with the simulation model. Consider applying this result to real-world situations: in that case, the ability to select the true temporal process will have a greater impact on state inference quality than the change state assumption.

Auto-correlation deteriorates the state decoding accuracy

The numerical experiments revealed that the deterioration in accuracy is primarily due to auto-correlation
in the observation layers rather than the hidden layer's Markovian structure. Although $d_T$ (hidden layer)
and $d_{AR}$ (observation layer) are co-linear in the settings with potential confounding effects, the current work
demonstrates that the modeling choices related to the observation layer are of primary importance for state
<figure><img src="image_18.png" /><figcaption>**Figure 15.** SUPPLEMENTARY. Delayed correlation plots for Vessel 1 (left, blue panels) and Vessel 2 (right, red panels). Persistent speed is in squares, and rotation speed is in triangles. State 1 is in empty symbols, and State 2 is in plain. Linear regressions for which p-value < 5% are represented.</figcaption></figure>
decoding performances. This finding is supported by both simulation-estimation experiments and real-data estimation.

## Simulation-estimation experiments

The HSMM<sub>SNB</sub> - AR1 model structure is the most robust to discrepancies between model hypotheses and data characteristics when retrieving hidden states, according to model robustness assessment (MRA). While the accuracy in state estimation deteriorates in 90% of cases when using a model structure different from the one used to simulate data, the relative loss is small (0.92% on average). Imposing no auto-correlation in the estimating model while allowing some for state-dependent variables has a negative impact on state decoding performance. The loss of accuracy is well explained by $d_{AR}$, which is the distance between AR1 and AR0 formulations. This is not true for the hidden layer, where imposing an HMM structure while the state process is semi-Markov has no negative impact on the state decoding performance. This finding has direct practical implications. It demonstrates that it is preferable to concentrate on the hypothesis choices of the observation layer rather than the hidden layer. In some ways, this is reassuring, as only the observation layer contains empirical data. Even though, in theory, the time dependence of the observed processes or the distinction of the PDFs should be evaluated conditionally on the hidden states (which is impossible in practice because the states are, by definition, hidden), we recommend conducting exploratory analyses of the observed process data to help setting as much as possible relevant auto-correlation hypotheses as performed on our real cases. However, this step requires access to (a sample) of trajectories with observations of the (normally hidden) states along the trajectories. If modeling is envisaged at the time of data collection, we recommend building up such a sample during a trajectory collection phase.

## Real cases

The analysis of model performance on real cases is consistent with the simulation results, indicating that discrepancies between the properties of the observations and the model's hidden layer assumptions do not affect performance. Both data sets demonstrate a visual deviation from the various model assumptions. However, even though the differences in residence time for Vessel 2 are greater, the model's performance improves. The use of the most flexible, and thus the most robust, model structure (i.e., HSMM<sub>SNB</sub> - AR1) produces state decoding performances that under-perform the simulation experiments. The losses in accuracies are the consequences of the discrepancies between the data characteristics and the model structure. These discrepancies and losses are low (resp. large) for vessel 2 (resp. vessel 1). Interestingly, Vessel 2's performances do not deteriorate strongly with the resolution, whereas Vessel 1's do (Figure 9). The two primary differences in favor of Vessel 2 concern the speed variables:

* The discriminating power between the speed distributions measured through the integrated metric $d_V$ is low, respectively large, for Vessel 1 and Vessel 2 (Figure 13). This contributes to the low accuracies obtained for Vessel 1.

* In most cases, the partial autocorrelation of observed speed is greater than one (the model's assumed autocorrelation), particularly for Vessel 1. Not only are the auto-correlations of order 1 larger, but the most likely order of the auto-regressive process modeling the speed variable is greater than 1.

Beyond these considerations about the value of coefficients of auto-correlation and the orders of the auto-regressive processes, the observations may not be consistent with an auto-regressive process at all. In particular, auto-regressive processes assume that the relationship between successive observations is linear, which may not be the case. As an example, we provide the cross-plots between successive speed components (Additional Figure 16), which violate linear relationships over time. To form the correct hypotheses about the observation process, start by ensuring that the data are consistent with an auto-regressive process. The violation of this assumption corresponds to a structural error process that is rarely evaluated. Before discussing
<figure><img src="image_19.png" /><figcaption>**Figure 16.** SUPPLEMENTARY. Settings 1. Evolution of the log-likelihood of each EM under MSA conditions (estimation model = simulation model. Rows correspond to the different model types. Columns correspond to the different resolutions. On each panel, the evolution of the log-likelihood during each of the 100 estimations is represented. The log-likelihood is normalized between the value for the initialization and the value obtained after convergence. The x-axis represents the number of iterations (different for each experiment). Tick marks represent the number of iterations reached at each estimation. When the EM reaches the maximum possible number of iterations, the output is represented by a black curve (which never happens as the numerical settings have been chosen so that this does not happen).</figcaption></figure>
the order of the AR, consider the more fundamental choice of using an AR formulation for the observation layer (Gloaguen et al., 2015; Joo et al., 2013; Morales et al., 2004). Second, within the context of auto-regressive processes, the data may not correspond to a zero or first AR order (AR0 or AR1). This second type of violation is caused by mis-specification and can be more easily controlled by evaluating the model's goodness of fit.

## Increasing the resolution: a tricky solution?

As expected, model characteristics are linked to data resolution. This affects not only the averages of the sojourn times, which are directly impacted because they are expressed in a series of time steps but also their distributions. This also affects the variability of the speeds (both persistent and rotational), their time dependencies, and the overlap of speed PDFs between states. The multidimensional nature of data resolution precludes statistical investigation of its effect on estimation performance, all other things being equal. Several authors (e.g., McClintock, King, et al., 2012; Vermard et al., 2010) have pointed out that time steps must be chosen to match the scale at which behavioral decisions are made using expert knowledge. Katara and Silva Katara and Silva, 2017 examined the mismatch between the temporal resolution of position records and the timescale of fishing activity. The authors investigated the effect of time resolution (from 10 minutes to 2 hours) on bias and errors in identifying fishing operations for vessels (Portuguese purse seiners) with characteristics similar to those of vessel 2. The results (7% of missed fishing sets with 2h time interval compared to 10 minutes) align with the estimated increase in accuracy with time scale resolution for Vessel 2.

Similar results were obtained for other vessels, emphasizing the importance of high-resolution time data Edelhoff et al., 2016; Mendo et al., 2019. Interestingly, when using data that were consistent with a given model's assumptions, such as in simulation-estimation experiments, increasing the resolution of simulated data resulted in higher accuracy. This was less obvious for real data, where the results differed between the two case studies and were not linear with data resolution. There is a contradiction between the need for high-frequency data to track short-term states or to get as close as possible to the tracked individual's actual behavior and the ability to model this behavior using a random decision. If we suppose that memory and/or inertia play a role in behavioral changes at high frequencies, in this case, it is clear that a state is dependent not only on previous states but also on multiple previous ones. As a result, a compromise must be made between the pursuit of better statistical conditions and a realism threshold to achieve an optimal level of data resolution.

In conclusion, in silico, when data characteristics are controlled by and consistent with model characteristics, simulation-estimation experiments show that a change in Markov property for the estimation model has no negative impact on state decoding performance unless one forces the state-dependent variable to be AR0 when they are indeed AR1. This aligns with the fact that $HSMM_{SNB}$ and AR1 formulations are generalizations of the usual HMM and AR0 formulations. This study also shows that model performance is more sensitive to the auto-regressive characteristics of the state-dependent variables than to the Markov properties of the hidden states. This is somewhat satisfactory because, in practice, one may have access to some aspects of the observation layer but not to the hidden, unobserved layer. This allows for empirical controls on what appear to be key model assumptions in terms of state decoding performance. However, in real cases, the robustness of $HSMM_{SNB}$ - AR1 appears low, with instances where its state decoding performances are worse than simple filtering. The confusion between the distributions of state-dependent variables in different states, whether observable (usually not possible) or more realistically speculated, is the most important factor influencing state decoding performance in real-world scenarios. Blind applications of state-space models in the absence of prior knowledge about the processes lead to radical errors and poor state-decoding performance. Given that HMMs are a popular choice in ecology due to their ease of use, we believe there is a need for serious consideration of the consequences of the numerous potential deviations between data and model assumptions.
<figure><img src="image_20.png" /><figcaption>**Figure 17.** SUPPLEMENTARY. Settings 2. Same caption as before.</figcaption></figure>

Additional information: On two occasions, the estimation of the shift for the shifted negative binomial PMFs causes a decreasing evolution of the likelihood. In these cases, we changed the seed to recover the normal behavior of the convergence. These cases correspond to experiment n° 4 - simulation n°72, and to experiment n° 12 - simulation n°73. See the code to get details (https://github.com/nicobez/H-S-MM_MSA).
# Appendices

## Supplementary figures

<figure><img src="image_19.png" /><figcaption>**Figure 15.** SUPPLEMENTARY. Boxplots of *n* parameter values of the *SNB* PMFs obtained when the simulation model is HMM (with AR0 or AR1 structures for the speed variables) and the estimation model is HSMM (with the same AR structures as that used for simulation ). Settings 1 (blue) and settings 2 (red). The white boxes correspond to state 1 while the coloured boxes are for state 2. Data resolution, is represented along the x-axis.</figcaption></figure>
Acknowledgements

527

The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM) for providing DATARMOR storage, data access, computational resources, visualization, web services, consultation, and support services. URL: http://www.ifremer.fr/pcdm.

528

529

530

A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions (https://creativecommons.org/licenses/by/4.0/legalcode.en).

531

532

533

534

Several colleagues have contributed to discussions and meetings without being formally co-authors. Their inputs are worth to be acknowledged: Pascal MONESTIEZ, David NERINI, Hélène De PONTUAL, Youen VER-MARD.

535

536

537

Fundings

538

The participation in the various workshops within which this work was discussed was financed by the INRAE network Path Tools and analysis - PathTIS.

539

540

Conflict of interest disclosure

541

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

542

543

Data, script, code, and supplementary information availability

544

Model simulations and states' estimations are done within R version 4.3.1 "Beagle Scouts" (R Core Team, 2022) using the package mhsmm (O'Connell and Højsgaard, 2011), which allows performing EM inferences with several realizations of the same underlying model. Modifications of the mhsmm package were, however, required. We extended the package:

545

546

547

548

* to the case of AR1 for the state-dependent variables that were not in the package,

549

* to the case of shifted negative binomial PMF for the sojourn time that was not fully implemented,

550

* to force independence between the state-dependent variables during the inference process (by default, a correlation is re-estimated during the EM procedure).

551

552

When considering 250 alternations of State 1 and 2, the 32 MSA experiments corresponding to the 2 settings × 4 degradations × 4 model types are fully reproducible from https://doi.org/10.5281/zenodo.10678877. Tolerances for stopping the iterative search of the maximum likelihood in the EM algorithm are adapted to avoid very long searches associated with teeny marginal improvements of the likelihood (the iterative procedure stops when the gain in log-likelihood is smaller than the tolerance). The tolerance is reduced to 10⁻⁰⁴ for HMM, and to 10⁻⁰³ for HSMM when the default values were 10⁻⁰⁸ for HMM, and 10⁻⁰⁴ for HSMM. Meanwhile, the maximum number of iterations is also adapted to ensure that the convergence is reached, at least with regard to the chosen tolerance. The maximum number of iterations is reduced to 200 for HMM (instead of 1000 by default) but increased to 150 for HSMM (instead of 100 by default). With these numerical settings, convergence is systematically reached before the maximum possible number of iterations (Figures 16 and 17).

553

554

555

556

557

558

559

560

561

562

563

With these numerical settings, the full MSA experiments last 10 hours and 20 minutes on an i7 PC which probably ruins their fluid and rapid full reproduction. This certainly does for the MRA procedures that are

564

565
<figure><img src="image_20.png" /><figcaption>**Figure 16.** SUPPLEMENTARY. Delayed correlation plots for Vessel 1 (left, blue panels) and Vessel 2 (right, red panels). The axes correspond to the speeds at time *t* (x-axis) and *t+1* hour (y-axis). This corresponds to the third and the fourth degradation for vessel 1 and 2 respectively. Persistent speed is in squares, and rotation speed is in triangles. State 1 is in empty symbols, and State 2 is in plain. Linear regressions for which p-value < 5% are represented.</figcaption></figure>
3 times longer. Thus, the MRA simulation-estimation experiments have been performed remotely using the computational resources of a super-computer (DATARMOR, http://www.ifremer.fr/pcdm). All the simulation-estimation outputs have been stored and the paper is reproducible based on this stored and open-access information. The code to re-run all the MRA simulations is, however, available on demand. After the upload of all the outputs of the simulation-estimation experiments, all the analyses and figures are reproducible from https://doi.org/10.5281/zenodo.10679448. For anonymity concerns, the coordinates of the trajectories are centered and standardized.

References

Allen J and P Somerfield (2009). A multivariate approach to model skill assessment. *Journal of Marine Systems* 76, 83-94. ISSN: 0924-7963. https://doi.org/10.1016/j.jmarsys.2008.05.009.

Auger-Méthé M, CCS Clair, MA Lewis, and AE Derocher (Aug. 2011). Sampling rate and misidentification of Lévy and non-Lévy movement paths: comment. en. *Ecology* 92, 1699-1701. ISSN: 0012-9658. https://doi.org/10.1890/10-1704.1.

Auger-Méthé M, K Newman, D Cole, F Empacher, R Gryba, AA King, V Leos-Barajas, J Mills Flemming, A Nielsen, G Petris, and L Thomas (2021). A guide to state-space modeling of ecological time series. en. *Ecological Monographs* 91, e01470. ISSN: 1557-7015. https://doi.org/10.1002/ecm.1470.

Avgar T, R Deardon, and JM Fryxell (Feb. 2013). An empirically parameterized individual based model of animal movement, perception, and memory. *Ecological Modelling* 251, 158-172. ISSN: 03043800. https://doi.org/10.1016/j.ecolmodel.2012.12.002.

Bertrand S, A Bertrand, R Guevara-Carrasco, and F Gerlotto (2007). Scale-Invariant Movements of Fishermen: The Same Foraging Strategy as Natural Predators. en. *Ecological Applications* 17, 331-337. ISSN: 1939-5582. https://doi.org/10.1890/06-0303.

Bez N, E Walker, D Gaertner, J Rivoirard, and P Gaspar (Nov. 2011). Fishing activity of tuna purse seiners estimated from vessel monitoring system (VMS) data. en. *Canadian Journal of Fisheries and Aquatic Sciences* 68. Ed. by Walters C, 1998-2010. ISSN: 0706-652X, 1205-7533. https://doi.org/10.1139/f2011-114.

Dempster AP, NM Laird, and DB Rubin (Sept. 1977). Maximum Likelihood from Incomplete Data Via the EM Algorithm. en. *Journal of the Royal Statistical Society: Series B (Methodological)* 39, 1-22. ISSN: 00359246. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x.

Edelhoff H, J Signer, and N Balkenhol (Dec. 2016). Path segmentation for beginners: an overview of current methods for detecting changes in animal movement patterns. en. *Movement Ecology* 4, 21. ISSN: 2051-3933. https://doi.org/10.1186/s40462-016-0086-5.

Gerritsen H and C Lordan (Jan. 2011). Integrating vessel monitoring systems (VMS) data with daily catch data from logbooks to explore the spatial distribution of catch and effort at high resolution. en. *ICES Journal of Marine Science* 68, 245-252. ISSN: 1095-9289, 1054-3139. https://doi.org/10.1093/icesjms/fsq137.

Glennie R, T Adam, V Leos-Barajas, T Michelot, T Photopoulou, and BT McClintock (2023). Hidden Markov models: Pitfalls and opportunities in ecology. en. *Methods in Ecology and Evolution* 14, 43-56. ISSN: 2041-210X. https://doi.org/10.1111/2041-210X.13801.

Gloaguen P, S Mahevas, E Rivot, M Woillez, J Guitton, Y Vermard, and MP Etienne (Feb. 2015). An autoregressive model to describe fishing vessel movement and activity. *Environmetrics* 26. WOS:000347896200003, 17-28. ISSN: 1180-4009. https://doi.org/10.1002/env.2319.

Guédon Y (2003). Estimating Hidden Semi-Markov Chains from Discrete Sequences. en. *Journal of Computational and Graphical Statistics* 12, 604-639.

Hellinger E (July 1909). Neue Begründung der Theorie quadratischer Formen von unendlichvielen Veränderlichen. en. *Journal für die reine und angewandte Mathematik* 1909, 210-271. ISSN: 1435-5345, 0075-4102. https://doi.org/10.1515/crll.1909.136.210.
<figure><img src="image_21.png" /><figcaption>**Figure 17.** SUPPLEMENTARY. Settings 1. Evolution of the log-likelihood of each EM under MSA conditions (estimation model = simulation model). Rows correspond to the different model types. Columns correspond to the different resolutions. On each panel, the evolution of the log-likelihood during each of the 100 estimations is represented. The log-likelihood is normalized between the value for the initialization and the value obtained after convergence. The x-axis represents the number of iterations (different for each experiment). Tick marks represent the number of iterations reached at each estimation. The EM never reached the maximum possible number of iterations set to 200 for HMM and 150 for HSMM.</figcaption></figure>
Johnson MJ and AS Willsky (2013). Bayesian Nonparametric Hidden Semi-Markov Models. *Journal of Machine Learning Research* 14, 673–701.

Jonsen I, M Basson, S Bestley, M Bravington, T Patterson, M Pedersen, R Thomson, U Thygesen, and S Wother- spoon (Apr. 2013). State-space models for bio-loggers: A methodological road map. en. *Deep Sea Research Part II: Topical Studies in Oceanography* 88-89, 34–46. ISSN: 09670645. https://doi.org/10.1016/j.dsr2.2012. 07.008.

Joo R, S Bertrand, J Tam, and R Fablet (2013). Hidden markov models: the best models for forager movements? *PloS one* 8, e71246.

Katara I and A Silva (Apr. 2017). Mismatch between VMS data temporal resolution and fishing activity time scales. en. *Fisheries Research* 188, 1–5. ISSN: 01657836. https://doi.org/10.1016/j.fishres.2016.11.023.

Langrock R, R King, J Matthiopoulos, L Thomas, D Fortin, and JM Morales (2012). Flexible and practical modeling of animal telemetry data: hidden Markov models and extensions. *Ecology* 93, 2336–2342. ISSN: 0012-9658. https://doi.org/10.1890/11-2241.1.

Lantuéjoul C (1991). Ergodicity and integral range. *Journal of Microscopy* 161, 387–403.

Lynch DR, DJJ McGillicuddy, and FE Werner (Feb. 2009). Skill assessment for coupled biological/physical models of marine systems. en. *Journal of Marine Systems* 76, 1–3. ISSN: 0924-7963. https://doi.org/10.1016/j. jmarsys.2008.05.002.

McClintock BT (2021). Worth the effort? A practical examination of random effects in hidden Markov models for animal telemetry data. en. *Methods in Ecology and Evolution* 12, 1475–1497. ISSN: 2041-210X. https://doi.org/10.1111/2041-210X.13619.

McClintock BT, R King, L Thomas, J Matthiopoulos, BJ McConnell, and JM Morales (2012). A general discrete- time modeling framework for animal movement using multistate random walks. *Ecological Monographs* 82, 335–349.

McClintock BT, R Langrock, O Gimenez, E Cam, DL Borchers, R Glennie, and TA Patterson (Dec. 2020). Uncov- ering ecological state dynamics with hidden Markov models. en. *Ecology Letters* 23. Ed. by Coulson T, 1878– 1903. ISSN: 1461-023X, 1461-0248. https://doi.org/10.1111/ele.13610.

Mendo T, S Smout, T Russo, L D'Andrea, and M James (Dec. 2019). Effect of temporal and spatial resolution on identification of fishing activities in small-scale fisheries using pots and traps. en. *ICES Journal of Marine Science* 76. Ed. by Maravelias C, 1601–1609. ISSN: 1054-3139, 1095-9289. https://doi.org/10.1093/icesjms/ fsz073.

Morales JM, DT Haydon, J Frair, KE Holsinger, and JM Fryxell (Sept. 2004). Extracting more out of relocation data: building movement models as mixtures of random walks. en. *Ecology* 85, 2436–2445. ISSN: 0012- 9658. https://doi.org/10.1890/03-0269.

Nettle D, MA Gibson, DW Lawson, and R Sear (2013). Human behavioral ecology: current research and future prospects. en. *Behavioral Ecology* 24, 1031–1040. ISSN: 1465-7279, 1045-2249. https://doi.org/10.1093/ beheco/ars222.

O'Connell J and S Højsgaard (2011). Hidden Semi Markov Models for Multiple Observation Sequences: The **mhsmm** Package for R. en. *Journal of Statistical Software* 39. ISSN: 1548-7660. https://doi.org/10.18637/jss. v039.i04.

Owen-Smith N, JM Fryxell, and EH Merrill (2010). Foraging theory upscaled: the behavioural ecology of herbi- vore movement. *Philosophical Transactions of the Royal Society B: Biological Sciences* 365, 2267–2278. https://doi.org/10.1098/rstb.2010.0095.

Pohle J, R Langrock, FM van Beest, and NM Schmidt (Sept. 2017). Selecting the Number of States in Hidden Markov Models: Pragmatic Solutions Illustrated Using Animal Movement. en. *Journal of Agricultural, Biological and Environmental Statistics* 22, 270–293. ISSN: 1085-7117, 1537-2693. https://doi.org/10.1007/s13253- 017-0283-8.

Poritz A (1982). Linear predictive hidden Markov models and the speech signal. *Proc.ICASSP '82. pp.1291-1294, Paris, France, May 1982.
<figure><img src="image_22.png" /><figcaption>**Figure 18.** SUPPLEMENTARY. Settings 2. Same caption as before.</figcaption></figure>

Additional information: On two occasions, the estimation of the shift for the shifted negative binomial PMFs causes a decreasing evolution of the likelihood. In these cases, we changed the seed to recover the normal behavior of the convergence. These cases correspond to experiment n° 4 - simulation n°72, and to experiment n° 12 - simulation n°73. See the code to get details (https://github.com/nicobez/H-S-MM_MSA).

**Acknowledgements**

The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM) for providing DATARMOR storage, data access, computational resources, visualization, web services, consultation, and support services.
R Core Team (2022). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing. Vienna, Austria.

Rabiner L (Feb. 1989). A tutorial on hidden Markov models and selected applications in speech recognition. *Proceedings of the IEEE* 77, 257-286. ISSN: 00189219. https://doi.org/10.1109/5.18626.

Schwager M, DM Anderson, Z Butler, and D Rus (Mar. 2007). Robust classification of animal tracking data. en. *Computers and Electronics in Agriculture* 56, 46-59. ISSN: 0168-1699. https://doi.org/10.1016/j.compag.2007.01.002.

Sonia Malefaki ST and N Limnios (2010). An EM and a Stochastic Version of the EM Algorithm for Nonparametric Hidden Semi-Markov Models. *Communications in Statistics - Simulation and Computation* 39, 240-261. https://doi.org/10.1080/03610910903411185.

Sur M, T Suffredini, SM Wessells, PH Bloom, M Lanzone, S Blackshire, S Sridhar, and T Katzner (2017). Improved supervised classification of accelerometry data to distinguish behaviors of soaring birds. en. *PLOS ONE* 12, e0174785. https://doi.org/10.1371/journal.pone.0174785.

Thums M, CJ Bradshaw, and MA Hindell (Dec. 2008). A validated approach for supervised dive classification in diving vertebrates. en. *Journal of Experimental Marine Biology and Ecology* 363, 75-83. ISSN: 0022-0981. https://doi.org/10.1016/j.jembe.2008.06.024.

Vermard Y, E Rivot, S Mahévas, P Marchal, and D Gascuel (2010). Identifying fishing trip behaviour and estimating fishing effort from VMS data using Bayesian Hidden Markov Models. *Ecological Modelling* 221, 1757-1769. ISSN: 0304-3800.

Yaglom A (1987). *Correlation Theory of Stationary and Related Random Functions*. en.

Zucchini W and IL MacDonald (2009). *Hidden Markov models for time series: an introduction using R*. Monographs on statistics and applied probability 110. OCLC: ocn144565537. Boca Raton: CRC Press. ISBN: 9781584885733.
URL: http://www.ifremer.fr/pcdm.

A CC-BY public copyright license has been applied by the authors to the present document and will be
applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in
accordance with the grant's open access conditions (https://creativecommons.org/licenses/by/4.0/legalcode.
en).

Several colleagues have contributed to discussions and meetings without being formally co-authors. Their inputs are worth to be acknowledged: Pascal MONESTIEZ, David NERINI, Hélène De PONTUAL, Youen VER-MARD.

The authors want to underline the very positive kindness in reviewers comments and requests with per-
sonal suggestions and their detailed analysis of the manuscript. Reading their feed backs and their comments
has been a real improvement and a pleasure. Would all the reviews be done with the same good-will and
rigour.

Fundings

The participation in the various workshops within which this work was discussed was financed by the INRAE network Path Tools and analysis - PathTIS.

**Conflict of interest disclosure**

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

Authors’ contributions

All the authors took part in the various workshops that led to this manuscript.

• Nicolas BEZ, contributed to the conception and the design of the analysis, to data analysis, and to write the paper. He also wrote all the R codes including some contributions of some co-authors.

* Pierre GLOAGUEN, contributed to the R code specific to the AR1 case, and to the proofreading of the paper.

* Marie-Pierre ETIENNE contributed to the conception and the design of part of the analysis. She also provided feedback on the analysis and writing.

* Rocio JOO led to some of the ideas of this paper, contributed to the first data analysis; she also provided critical feedbacks on the manuscript.

* Sophie LANCO provided critical feedbacks on the ideas and on the manuscript and took part in the various workshops that led to this manuscript.

* Etienne RIVOT provided critical feedback on the ideas and on the manuscript.

* Emily WALKER contributed to the initial data sampling and to the first data analysis, she also provided critical feedback on the manuscript.

* Mathieu WOILLEZ provided feedbacks on the manuscript.

* Stéphanie MAHEVAS, contributed to the conception and the design of the analysis, to the computation on the super-calculator, to data analysis, and to write the paper.
<table><caption>**Table 1.** Definition of the notations</caption><thead><tr><th></th><th>Notation</th><th>Meaning</th></tr></thead><tbody><tr><td></td><td>t</td><td>Discrete time step</td></tr><tr><td>Hidden layer</td><td>s<sub>t</sub>, S<sub>t</sub></td><td>State at time step t (lower case: observed; upper case:random variable)</td></tr><tr><td></td><td>N<sub>S</sub></td><td>Number of possible states</td></tr><tr><td></td><td>t<sub>s</sub>, T<sub>s</sub></td><td>Sojourn time in state s (lower case: observed; upper case:random variable)</td></tr><tr><td></td><td>p<sub>s,s</sub></td><td>For HMM, probability to stay in state s when moving from t to t + 1 so that T<sub>s</sub> ~ geom(1 - p<sub>s,s</sub>)</td></tr><tr><td></td><td>N<sub>pS</sub></td><td>Number of parameters of the PMFs of the sojourn time</td></tr><tr><td></td><td>shift<sub>s</sub></td><td>Shift for the PMF of T<sub>s</sub> for HSMM</td></tr><tr><td></td><td>(n<sub>s</sub>, p<sub>s</sub>)</td><td>Parameters of the negative binomial of the PMF of T<sub>s</sub> for HSMM so that T<sub>s</sub> ~ shift<sub>s</sub> + N<sub>B</sub>(n<sub>s</sub>, p<sub>s</sub>)</td></tr><tr><td></td><td>N<sub>H</sub></td><td>Number of parameters for the hidden layer</td></tr><tr><td></td><td>d<sub>T<sub>s</sub></sub></td><td>Hellinger distance between the geometric and the shifted negative binomial PMFs of T<sub>s</sub></td></tr><tr><td></td><td>d<sub>T</sub></td><td>Mean Hellinger distance between the geometric and the shifted negative binomial PMFs over all possible states</td></tr><tr><td></td><td>π<sub>s</sub></td><td>Probability to be in state s</td></tr><tr><td>Observation layer</td><td>N<sub>V</sub></td><td>Number of state-dependent variables</td></tr><tr><td></td><td>ν<sub>i</sub>, V<sub>i</sub></td><td>State-dependent variable number i (lower case: observed; upper case:random variable)</td></tr><tr><td></td><td>(V<sub>p</sub>, V<sub>r</sub>)</td><td>Persistent and rotation speeds (V<sub>1</sub>, V<sub>2</sub>)</td></tr><tr><td></td><td>V<sub>i,t</sub></td><td>Temporal process for state-dependent variables</td></tr><tr><td></td><td>N<sub>(μ<sub>i,s</sub>, σ<sub>i,s</sub>)</sub> = N<sub>i,s</sub></td><td>Uni-variate Gaussian PDFs of the state-dependent variable i for state s</td></tr><tr><td></td><td>μ<sub>i,s</sub></td><td>Expected value of V<sub>i,s</sub></td></tr><tr><td></td><td>σ<sub>i,s</sub></td><td>Standard deviation of V<sub>i,s</sub></td></tr><tr><td></td><td>N<sub>p0</sub></td><td>Number of parameters of the PDFs of the state-dependent variables</td></tr><tr><td></td><td>N<sub>O</sub></td><td>Number of parameters for the observation layer</td></tr><tr><td></td><td>d<sub>V</sub></td><td>Discrimination between PDFs</td></tr><tr><td></td><td>τ<sub>i,s</sub></td><td>Coefficient of correlation for auto-regressive processes of order 1 for state-dependent variable i in state s</td></tr><tr><td></td><td>η<sub>i,s</sub></td><td>μ = <sup>eta</sup>/<sub>1 - tau</sub></td></tr><tr><td></td><td>A</td><td>Integral range of AR1 process</td></tr><tr><td></td><td>d<sub>AR<sub>i,s</sub></sub></td><td>Difference between AR0 and AR1 formulations for state-dependent variable i and state s</td></tr><tr><td></td><td>d<sub>AR<sub>s</sub></sub></td><td>Difference between AR0 and AR1 formulations for state-dependent variables in state s</td></tr><tr><td></td><td>d<sub>AR</sub></td><td>Difference between AR0 and AR1 formulations</td></tr></tbody></table>
## Data, script, code, and supplementary information availability

Model simulations and states' estimations are done within R version 4.3.1 "Beagle Scouts" (R Core Team, 2022) using the package mhsmm (O'Connell and Højsgaard, 2011), which allows performing EM inferences with several realizations of the same underlying model. Modifications of the mhsmm package were, however, required. We extended the package:

* to the case of AR1 for the state-dependent variables that were not in the package,

* to the case of shifted negative binomial PMF for the sojourn time that was not fully implemented,

* to force independence between the state-dependent variables during the inference process (by default, a correlation is re-estimated during the EM procedure).

When considering 250 alternations of State 1 and 2, the 32 MSA experiments corresponding to the 2 settings × 4 degradations × 4 model types are fully reproducible from https://doi.org/10.5281/zenodo.10678877. Tolerances for stopping the iterative search of the maximum likelihood in the EM algorithm are adapted to avoid very long searches associated with tiny marginal improvements of the likelihood (the iterative procedure stops when the gain in log-likelihood is smaller than the tolerance). The tolerance is reduced to $10^{-04}$ for HMM, and to $10^{-03}$ for HSMM when the default values were $10^{-08}$ for HMM, and $10^{-04}$ for HSMM. Meanwhile, the maximum number of iterations is also adapted to ensure that the convergence is reached, at least with regard to the chosen tolerance. The maximum number of iterations is reduced to 200 for HMM (instead of 1000 by default) but increased to 150 for HSMM (instead of 100 by default). With these numerical settings, convergence is systematically reached before the maximum possible number of iterations (Figures 17 and 18).

With these numerical settings, the full MSA experiments last 10 hours and 20 minutes on an i7 PC which probably ruins their fluid and rapid full reproduction. This certainly does for the MRA procedures that are 3 times longer. Thus, the MRA simulation-estimation experiments have been performed remotely using the computational resources of a super-computer (DATARMOR, http://www.ifremer.fr/pcdm). All the simulation-estimation outputs have been stored and the paper is reproducible based on this stored and open-access information. The code to re-run all the MRA simulations is, however, available on demand. After the upload of all the outputs of the simulation-estimation experiments, all the analyses and figures are reproducible from https://doi.org/10.5281/zenodo.10679448. For anonymity concerns, the coordinates of the trajectories are centered and standardized.

## References

Allen J and P Somerfield (2009). A multivariate approach to model skill assessment. *Journal of Marine Systems* 76, 83-94. ISSN: 0924-7963. https://doi.org/10.1016/j.jmarsys.2008.05.009.

Auger-Méthé M, CCS Clair, MA Lewis, and AE Derocher (Aug. 2011). Sampling rate and misidentification of Lévy and non-Lévy movement paths: comment. en. *Ecology* 92, 1699-1701. ISSN: 0012-9658. https://doi.org/10.1890/10-1704.1.

Auger-Méthé M, K Newman, D Cole, F Empacher, R Gryba, AA King, V Leos-Barajas, J Mills Flemming, A Nielsen, G Petris, and L Thomas (2021). A guide to state-space modeling of ecological time series. *Ecological Monographs* 91, e01470. ISSN: 1557-7015. https://doi.org/10.1002/ecm.1470.

Avgar T, R Deardon, and JM Fryxell (Feb. 2013). An empirically parameterized individual based model of animal movement, perception, and memory. *Ecological Modelling* 251, 158-172. ISSN: 03043800. https://doi.org/10.1016/j.ecolmodel.2012.12.002.

Barbu VS and N Limnios (2008). *Semi-Markov Chains and Hidden Semi-Markov Models toward Applications*. Springer.
<table><caption>**Table 2.** Statistics for Vessel 1. For the sojourn time in the various states, the parameters of the PMFs in the model settings (e.g. $p_{s,s}, n_s, p_s$) are deduced from them according to the equations provided in the method section. The statistics for the sojourn time are provided in the number of time steps. For the speed variables, they are expressed in knots.</caption><thead><tr><th></th><th colspan="2">Resolution 15 min (N=9506)</th><th colspan="2">Resolution 30 min (N=4748)</th><th colspan="2">Resolution 60 min (N=2370)</th><th colspan="2">Resolution 120 min (N=1171)</th></tr></thead><tbody><tr><td></td><td colspan="8" style="text-align: center;">T</td></tr><tr><td></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td></tr><tr><td>E[T<sub>s</sub>]</td><td>11</td><td>24</td><td>5.7</td><td>12</td><td>2.9</td><td>6.1</td><td>1.6</td><td>3.5</td></tr><tr><td>σ[T<sub>s</sub>]</td><td>7.4</td><td>11</td><td>3.9</td><td>5.7</td><td>2</td><td>3</td><td>1.1</td><td>2.1</td></tr><tr><td>shift<sub>s</sub></td><td>2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>1</sub> = V<sub>p</sub></td></tr><tr><td></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td></tr><tr><td>μ<sub>1,s</sub></td><td>4.5</td><td>3</td><td>3.8</td><td>2.4</td><td>2.7</td><td>1.4</td><td>1.8</td><td>0.3</td></tr><tr><td>σ<sub>1,s</sub></td><td>4.6</td><td>1.6</td><td>4.9</td><td>2.0</td><td>5.1</td><td>2.4</td><td>4.7</td><td>2.2</td></tr><tr><td>τ<sub>1,s</sub></td><td>0.60</td><td>0.34</td><td>0.46</td><td>0.15</td><td>0.43</td><td>0.17</td><td>0.41</td><td>0.30</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>2</sub> = V<sub>r</sub></td></tr><tr><td></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td></tr><tr><td>μ<sub>2,s</sub></td><td>0.07</td><td>0.13</td><td>0.11</td><td>0.24</td><td>0.08</td><td>0.22</td><td>0.12</td><td>0.009</td></tr><tr><td>σ<sub>2,s</sub></td><td>1.7</td><td>1.4</td><td>2</td><td>1.5</td><td>2.2</td><td>1.4</td><td>2.3</td><td>1.3</td></tr><tr><td>τ<sub>2,s</sub></td><td>0.15</td><td>0.52</td><td>0.11</td><td>0.50</td><td>0.07</td><td>0.43</td><td>-0.008</td><td>0.27</td></tr></tbody></table>

<table><caption>**Table 3.** Statistics for Vessel 2. For the sojourn time in the various states, the parameters of the PMFs in the model settings (e.g. $p_{s,s}, n_s, p_s$) are deduced from them according to the equations provided in the method section. The statistics for the sojourn time are provided in the number of time steps. For the speed variables, they are expressed in knots.</caption><thead><tr><th></th><th colspan="2">Resolution 5 min (N=6181)</th><th colspan="2">Resolution 15 min (N=2011)</th><th colspan="2">Resolution 30 min (N=968)</th><th colspan="2">Resolution 60 min (N=446)</th></tr></thead><tbody><tr><td></td><td colspan="8" style="text-align: center;">T</td></tr><tr><td></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td></tr><tr><td>E[T<sub>s</sub>]</td><td>104</td><td>22</td><td>41</td><td>7.4</td><td>20</td><td>4.4</td><td>11</td><td>2.3</td></tr><tr><td>σ[T<sub>s</sub>]</td><td>112</td><td>4.4</td><td>45</td><td>1.9</td><td>22</td><td>1.8</td><td>12</td><td>1.2</td></tr><tr><td>shift<sub>s</sub></td><td>4</td><td>15</td><td>1</td><td>4</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>1</sub> = V<sub>p</sub></td></tr><tr><td></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td></tr><tr><td>μ<sub>1,s</sub></td><td>9.3</td><td>0.46</td><td>8.3</td><td>0.43</td><td>7.3</td><td>0.41</td><td>6.5</td><td>0.32</td></tr><tr><td>σ<sub>1,s</sub></td><td>4.7</td><td>1.2</td><td>5.2</td><td>0.9</td><td>5.5</td><td>1.3</td><td>5.5</td><td>1.6</td></tr><tr><td>τ<sub>1,s</sub></td><td>0.70</td><td>0.56</td><td>0.49</td><td>0.27</td><td>0.44</td><td>0.12</td><td>0.35</td><td>0.20</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>2</sub> = V<sub>r</sub></td></tr><tr><td></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td></tr><tr><td>μ<sub>2,s</sub></td><td>0.13</td><td>-0.08</td><td>-0.03</td><td>0.11</td><td>-0.24</td><td>0.20</td><td>-0.31</td><td>0.20</td></tr><tr><td>σ<sub>2,s</sub></td><td>2.9</td><td>0.80</td><td>3.8</td><td>0.80</td><td>4.3</td><td>1.3</td><td>4.6</td><td>1.6</td></tr><tr><td>τ<sub>2,s</sub></td><td>0.17</td><td>-0.10</td><td>0.12</td><td>-0.04</td><td>0.07</td><td>-0.09</td><td>0.04</td><td>-0.07</td></tr></tbody></table>
Bertrand S, A Bertrand, R Guevara-Carrasco, and F Gerlotto (2007). Scale-Invariant Movements of Fishermen: The Same Foraging Strategy as Natural Predators. en. *Ecological Applications* 17, 331-337. ISSN: 1939-5582. https://doi.org/10.1890/06-0303.

Bez N, E Walker, D Gaertner, J Rivoirard, and P Gaspar (Nov. 2011). Fishing activity of tuna purse seiners estimated from vessel monitoring system (VMS) data. en. *Canadian Journal of Fisheries and Aquatic Sciences* 68. Ed. by Walters C, 1998-2010. ISSN: 0706-652X, 1205-7533. https://doi.org/10.1139/f2011-114.

Dempster AP, NM Laird, and DB Rubin (Sept. 1977). Maximum Likelihood from Incomplete Data Via the EM Algorithm. en. *Journal of the Royal Statistical Society: Series B (Methodological)* 39, 1-22. ISSN: 00359246. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x.

Edelhoff H, J Signer, and N Balkenhol (Dec. 2016). Path segmentation for beginners: an overview of current methods for detecting changes in animal movement patterns. en. *Movement Ecology* 4, 21. ISSN: 2051-3933. https://doi.org/10.1186/s40462-016-0086-5.

Gerritsen H and C Lordan (Jan. 2011). Integrating vessel monitoring systems (VMS) data with daily catch data from logbooks to explore the spatial distribution of catch and effort at high resolution. en. *ICES Journal of Marine Science* 68, 245-252. ISSN: 1095-9289, 1054-3139. https://doi.org/10.1093/icesjms/fsq137.

Glennie R, T Adam, V Leos-Barajas, T Michelot, T Photopoulou, and BT McClintock (2023). Hidden Markov models: Pitfalls and opportunities in ecology. en. *Methods in Ecology and Evolution* 14, 43-56. ISSN: 2041-210X. https://doi.org/10.1111/2041-210X.13801.

Gloaguen P, S Mahevas, E Rivot, M Woillez, J Guitton, Y Vermard, and MP Etienne (Feb. 2015). An autoregressive model to describe fishing vessel movement and activity. *Environmetrics* 26. WOS:000347896200003, 17-28. ISSN: 1180-4009. https://doi.org/10.1002/env.2319.

Guédon Y (2003). Estimating Hidden Semi-Markov Chains from Discrete Sequences. en. *Journal of Computational and Graphical Statistics* 12, 604-639.

Hellinger E (July 1909). Neue Begründung der Theorie quadratischer Formen von unendlichvielen Veränderlichen. en. *Journal für die reine und angewandte Mathematik* 1909, 210-271. ISSN: 1435-5345, 0075-4102. https://doi.org/10.1515/crll.1909.136.210.

Hoek J van der and RJ Elliott (2018). *Introduction to Hidden Semi-Markov Models*. Cambridge University Press.

Johnson MJ and AS Willsky (2013). Bayesian Nonparametric Hidden Semi-Markov Models. *Journal of Machine Learning Research* 14, 673-701.

Jonsen I, M Basson, S Bestley, M Bravington, T Patterson, M Pedersen, R Thomson, U Thygesen, and S Wotherspoon (Apr. 2013). State-space models for bio-loggers: A methodological road map. en. *Deep Sea Research Part II: Topical Studies in Oceanography* 88-89, 34-46. ISSN: 09670645. https://doi.org/10.1016/j.dsr2.2012.07.008.

Joo R, S Bertrand, J Tam, and R Fablet (2013). Hidden markov models: the best models for forager movements? *PloS one* 8, e71246.

Katara I and A Silva (Apr. 2017). Mismatch between VMS data temporal resolution and fishing activity time scales. en. *Fisheries Research* 188, 1-5. ISSN: 01657836. https://doi.org/10.1016/j.fishres.2016.11.023.

Langrock R, R King, J Matthiopoulos, L Thomas, D Fortin, and JM Morales (2012). Flexible and practical modeling of animal telemetry data: hidden Markov models and extensions. *Ecology* 93, 2336-2342. ISSN: 0012-9658. https://doi.org/10.1890/11-2241.1.

Lantuéjoul C (1991). Ergodicity and integral range. *Journal of Microscopy* 161, 387-403.

Lynch DR, DJJ McGillicuddy, and FE Werner (Feb. 2009). Skill assessment for coupled biological/physical models of marine systems. en. *Journal of Marine Systems* 76, 1-3. ISSN: 0924-7963. https://doi.org/10.1016/j.jmarsys.2008.05.002.

McClintock BT (2021). Worth the effort? A practical examination of random effects in hidden Markov models for animal telemetry data. en. *Methods in Ecology and Evolution* 12, 1475-1497. ISSN: 2041-210X. https://doi.org/10.1111/2041-210X.13619.
McClintock BT, R King, L Thomas, J Matthiopoulos, BJ McConnell, and JM Morales (2012). A general discrete-time modeling framework for animal movement using multistate random walks. *Ecological Monographs* 82, 335–349.

McClintock BT, R Langrock, O Gimenez, E Cam, DL Borchers, R Glennie, and TA Patterson (Dec. 2020). Uncov-ering ecological state dynamics with hidden Markov models. en. *Ecology Letters* 23. Ed. by Coulson T, 1878–1903. ISSN: 1461-023X, 1461-0248. https://doi.org/10.1111/ele.13610.

Mendo T, S Smout, T Russo, L D'Andrea, and M James (Dec. 2019). Effect of temporal and spatial resolution on identification of fishing activities in small-scale fisheries using pots and traps. en. *ICES Journal of Marine Science* 76. Ed. by Maravelias C, 1601–1609. ISSN: 1054-3139, 1095-9289. https://doi.org/10.1093/icesjms/ fsz073.

Morales JM, DT Haydon, J Frair, KE Holsinger, and JM Fryxell (Sept. 2004). Extracting more out of relocation data: building movement models as mixtures of random walks. en. *Ecology* 85, 2436–2445. ISSN: 0012-9658. https://doi.org/10.1890/03-0269.

Nettle D, MA Gibson, DW Lawson, and R Sear (2013). Human behavioral ecology: current research and future prospects. en. *Behavioral Ecology* 24, 1031–1040. ISSN: 1465-7279, 1045-2249. https://doi.org/10.1093/beheco/ars222.

O'Connell J and S Højsgaard (2011). Hidden Semi Markov Models for Multiple Observation Sequences: The **mhsmm** Package for R. en. *Journal of Statistical Software* 39. ISSN: 1548-7660. https://doi.org/10.18637/jss.v039.i04.

Owen-Smith N, JM Fryxell, and EH Merrill (2010). Foraging theory upscaled: the behavioural ecology of herbi-vore movement. *Philosophical Transactions of the Royal Society B: Biological Sciences* 365, 2267–2278. https://doi.org/10.1098/rstb.2010.0095.

Pohle J, R Langrock, FM van Beest, and NM Schmidt (Sept. 2017). Selecting the Number of States in Hidden Markov Models: Pragmatic Solutions Illustrated Using Animal Movement. en. *Journal of Agricultural, Biological and Environmental Statistics* 22, 270–293. ISSN: 1085-7117, 1537-2693. https://doi.org/10.1007/s13253-017-0283-8.

Poritz A (1982). Linear predictive hidden Markov models and the speech signal. *Proc.ICASSP '82*. pp.1291-1294, Paris, France, May 1982.

R Core Team (2022). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing. Vienna, Austria.

Rabiner L (Feb. 1989). A tutorial on hidden Markov models and selected applications in speech recognition. *Proceedings of the IEEE* 77, 257–286. ISSN: 00189219. https://doi.org/10.1109/5.18626.

Saporta G (1990). *Probabilités, analyses des données et statistique*. Technip.

Schwager M, DM Anderson, Z Butler, and D Rus (Mar. 2007). Robust classification of animal tracking data. en. *Computers and Electronics in Agriculture* 56, 46–59. ISSN: 0168-1699. https://doi.org/10.1016/j.compag.2007.01.002.

Sonia Malefaki ST and N Limnios (2010). An EM and a Stochastic Version of the EM Algorithm for Nonparametric Hidden Semi-Markov Models. *Communications in Statistics - Simulation and Computation* 39, 240–261. https://doi.org/10.1080/03610910903411185.

Sur M, T Suffredini, SM Wessells, PH Bloom, M Lanzone, S Blackshire, S Sridhar, and T Katzner (2017). Improved supervised classification of accelerometry data to distinguish behaviors of soaring birds. en. *PLOS ONE* 12, e0174785. https://doi.org/10.1371/journal.pone.0174785.

Thums M, CJ Bradshaw, and MA Hindell (Dec. 2008). A validated approach for supervised dive classification in diving vertebrates. en. *Journal of Experimental Marine Biology and Ecology* 363, 75–83. ISSN: 0022-0981. https://doi.org/10.1016/j.jembe.2008.06.024.

Vermard Y, E Rivot, S Mahévas, P Marchal, and D Gascuel (2010). Identifying fishing trip behaviour and estimat-ing fishing effort from VMS data using Bayesian Hidden Markov Models. *Ecological Modelling* 221, 1757–1769. ISSN: 0304-3800.
Yaglom A (1987). *Correlation Theory of Stationary and Related Random Functions*. Springer.

Zucchini W and IL MacDonald (2009). *Hidden Markov models for time series: an introduction using R*. Monographs on statistics and applied probability 110. OCLC: ocn144565537. Boca Raton: CRC Press. ISBN: 9781584885733.
<table><caption>**Table 1.** Definition of the notations</caption><thead><tr><th></th><th>t</th><th>Discrete time step</th></tr></thead><tbody><tr><td>Hidden layer</td><td>N<sub>S</sub></td><td>Number of possible states</td></tr><tr><td></td><td>s</td><td>State code number; s = 1, ..., N<sub>S</sub></td></tr><tr><td></td><td>s<sub>t</sub>, S<sub>t</sub></td><td>(HMM) State at time step t (lower case: observed; upper case: random variable); S<sub>t</sub> ~ Markov process</td></tr><tr><td></td><td>˜s<sub>n</sub>, ˜S<sub>n</sub></td><td>(HSMM) n<sup>th</sup> visited state (lower case: observed; upper case: random variable); ˜S<sub>t</sub> ~ Markov chain</td></tr><tr><td></td><td>t<sub>s</sub>, T<sub>s</sub></td><td>Sojourn time in state s (lower case: observed; upper case: random variable)</td></tr><tr><td></td><td>p<sub>ss'</sub></td><td>(HMM) probability to switch from state s to state s' when moving from t to t + 1</td></tr><tr><td></td><td>p̃<sub>ss'</sub></td><td>(HSMM) probability to switch from state s to state s' when leaving state s; p̃<sub>ss</sub> = 0</td></tr><tr><td></td><td>q<sub>s</sub></td><td>(HMM) probability to leave state s when moving from t to t + 1</td></tr><tr><td></td><td>G, SNB</td><td>(HSMM) Third parameter of the SNB PMF for state s. Corresponds to the q<sub>s</sub> parameter of the HMM case when shift<sub>s</sub> = 1 and when n<sub>s</sub> = 1.</td></tr><tr><td></td><td>(shift<sub>s</sub>, n<sub>s</sub>, q<sub>s</sub>)</td><td>Geometric (HMM) and shifted negative binomial (HSMM) PMFs of T<sub>s</sub></td></tr><tr><td></td><td>d<sub>T</sub></td><td>Hellinger distance between the geometric and the shifted negative binomial PMFs of T<sub>s</sub></td></tr><tr><td></td><td>d<sub>T</sub></td><td>Mean Hellinger distance between the geometric and the shifted negative binomial PMFs over all possible states</td></tr><tr><td></td><td>π<sub>s</sub>, Π<sub>s</sub></td><td>Probability to be in state s (lower case: empirical; upper case: model)</td></tr><tr><td></td><td>N<sub>pS</sub></td><td>Number of parameters of the PMFs of the sojourn time</td></tr><tr><td></td><td>N<sub>H</sub></td><td>Number of parameters defining the hidden layer</td></tr><tr><td>Observation layer</td><td>N<sub>V</sub></td><td>Number of state-dependent variables</td></tr><tr><td></td><td>ν<sub>i</sub>, V<sub>i</sub></td><td>State-dependent variable number i (lower case: observed; upper case: random variable); typically, (V<sub>p</sub>, V<sub>r</sub>), the persistent and rotation speeds</td></tr><tr><td></td><td>V<sub>i,t</sub></td><td>Random process for state-dependent variables</td></tr><tr><td></td><td>N<sup>g</sup>(μ<sub>i,s</sub>, σ<sub>i,s</sub>)</td><td>Uni-variate Gaussian PDFs of the state-dependent variable i for state s: [V<sub>i,t</sub>|S<sub>t</sub> = s] ~ N<sup>g</sup>(μ<sub>i,s</sub>, σ<sub>i,s</sub>)</td></tr><tr><td></td><td>N<sub>i,s</sub>(ν)</td><td>PDF of the state-dependent variable given the state s: f<sub>V<sub>i</sub>|S<sub>i</sub>=s</sub>(ν) = N<sub>i,s</sub>(ν)</td></tr><tr><td></td><td>μ<sub>i,s</sub>, σ<sub>i,s</sub><sup>2</sup></td><td>Expected value and variance of V<sub>i,t</sub> knowing the state: μ<sub>i,s</sub> = E(V<sub>i,t</sub>|S<sub>t</sub> = s), σ<sub>i,s</sub><sup>2</sup> = V(V<sub>i,t</sub>|S<sub>t</sub> = s)</td></tr><tr><td></td><td>N<sub>p0</sub></td><td>Number of parameters of the PDFs of the state-dependent variables</td></tr><tr><td></td><td>N<sub>O</sub></td><td>Number of parameters for the observation layer</td></tr><tr><td></td><td>R<sub>S</sub>, r<sub>s</sub></td><td>Region of ℝ<sup>2</sup> where the state is estimated as s in a mixture model: R<sub>s</sub> = {V = (ν<sub>1</sub>, ν<sub>2</sub>), such that Š = s} (r<sub>s</sub> its empirical version)</td></tr><tr><td></td><td>d<sup>emp</sup><sub>V</sub>, d<sup>mod</sup><sub>V</sub></td><td>Discrimination coefficient between PDFs (<sup>emp</sup>: based on empirical speed histograms; <sup>mod</sup>: based on PDFs)</td></tr><tr><td></td><td>τ<sub>i,s</sub></td><td>Coefficient of correlation for auto-regressive processes of order 1 for state-dependent variable i in state s</td></tr><tr><td></td><td>η<sub>i,s</sub></td><td>μ = <sup>eta</sup>&frasl;<sub>1-tau</sub></td></tr><tr><td></td><td>A</td><td>Integral range of AR1 process</td></tr><tr><td></td><td>d<sub>AR<sub>i,s</sub></sub></td><td>Difference between AR0 and AR1 formulations for state-dependent variable i and state s</td></tr><tr><td></td><td>d<sub>AR<sub>s</sub></sub></td><td>Difference between AR0 and AR1 formulations for state-dependent variables in state s</td></tr><tr><td></td><td>d<sub>AR</sub></td><td>Difference between AR0 and AR1 formulations</td></tr></tbody></table>
<table><caption>**Table 2.** Definition of the notations (continue)</caption><tr><td>Models</td><td>m<sub>1</sub></td><td>HMM-AR0</td></tr><tr><td></td><td>m<sub>2</sub></td><td>HMM-AR1</td></tr><tr><td></td><td>m<sub>3</sub></td><td>HSMM-AR0</td></tr><tr><td></td><td>m<sub>4</sub></td><td>HSMM-AR1</td></tr><tr><td></td><td>θ, ̂θ</td><td>Random vector of parameters and its estimator</td></tr><tr><td></td><td>N<sub>θ</sub></td><td>Number of model parameters N<sub>θ</sub> = N<sub>H</sub> + N<sub>0</sub></td></tr></table>
<table><caption>**Table 3.** Statistics for Vessel 1. For the sojourn time in the various states, the parameters of the PMFs in the model settings (e.g. $p_{s,s}, n_s, p_s$) are deduced from them according to the equations provided in the method section. The statistics for the sojourn time are provided in the number of time steps. For the speed variables, they are expressed in knots.</caption><thead><tr><th></th><th colspan="2">Resolution 15 min<br>(N=9506)</th><th colspan="2">Resolution 30 min<br>(N=4748)</th><th colspan="2">Resolution 60 min<br>(N=2370)</th><th colspan="2">Resolution 120 min<br>(N=1171)</th></tr></thead><tbody><tr><td></td><td colspan="8" style="text-align: center;">T</td></tr><tr><td></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td></tr><tr><td>E[T<sub>s</sub>]</td><td>11</td><td>24</td><td>5.7</td><td>12</td><td>2.9</td><td>6.1</td><td>1.6</td><td>3.5</td></tr><tr><td>σ[T<sub>s</sub>]</td><td>7.4</td><td>11</td><td>3.9</td><td>5.7</td><td>2</td><td>3</td><td>1.1</td><td>2.1</td></tr><tr><td>shift<sub>s</sub></td><td>2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>1</sub> = V<sub>p</sub></td></tr><tr><td></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td></tr><tr><td>μ<sub>1,s</sub></td><td>4.5</td><td>3</td><td>3.8</td><td>2.4</td><td>2.7</td><td>1.4</td><td>1.8</td><td>0.3</td></tr><tr><td>σ<sub>1,s</sub></td><td>4.6</td><td>1.6</td><td>4.9</td><td>2.0</td><td>5.1</td><td>2.4</td><td>4.7</td><td>2.2</td></tr><tr><td>τ<sub>1,s</sub></td><td>0.60</td><td>0.34</td><td>0.46</td><td>0.15</td><td>0.43</td><td>0.17</td><td>0.41</td><td>0.30</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>2</sub> = V<sub>r</sub></td></tr><tr><td></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td></tr><tr><td>μ<sub>2,s</sub></td><td>0.07</td><td>0.13</td><td>0.11</td><td>0.24</td><td>0.08</td><td>0.22</td><td>0.12</td><td>0.009</td></tr><tr><td>σ<sub>2,s</sub></td><td>1.7</td><td>1.4</td><td>2</td><td>1.5</td><td>2.2</td><td>1.4</td><td>2.3</td><td>1.3</td></tr><tr><td>τ<sub>2,s</sub></td><td>0.15</td><td>0.52</td><td>0.11</td><td>0.50</td><td>0.07</td><td>0.43</td><td>-0.008</td><td>0.27</td></tr></tbody></table>

<table><caption>**Table 4.** Statistics for Vessel 2. For the sojourn time in the various states, the parameters of the PMFs in the model settings (e.g. $p_{s,s}, n_s, p_s$) are deduced from them according to the equations provided in the method section. The statistics for the sojourn time are provided in the number of time steps. For the speed variables, they are expressed in knots.</caption><thead><tr><th></th><th colspan="2">Resolution 5 min<br>(N=6181)</th><th colspan="2">Resolution 15 min<br>(N=2011)</th><th colspan="2">Resolution 30 min<br>(N=968)</th><th colspan="2">Resolution 60 min<br>(N=446)</th></tr></thead><tbody><tr><td></td><td colspan="8" style="text-align: center;">T</td></tr><tr><td></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td><td>T<sub>1</sub></td><td>T<sub>2</sub></td></tr><tr><td>E[T<sub>s</sub>]</td><td>104</td><td>22</td><td>41</td><td>7.4</td><td>20</td><td>4.4</td><td>11</td><td>2.3</td></tr><tr><td>σ[T<sub>s</sub>]</td><td>112</td><td>4.4</td><td>45</td><td>1.9</td><td>22</td><td>1.8</td><td>12</td><td>1.2</td></tr><tr><td>shift<sub>s</sub></td><td>4</td><td>15</td><td>1</td><td>4</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>1</sub> = V<sub>p</sub></td></tr><tr><td></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td><td>V<sub>1,1</sub></td><td>V<sub>1,2</sub></td></tr><tr><td>μ<sub>1,s</sub></td><td>9.3</td><td>0.46</td><td>8.3</td><td>0.43</td><td>7.3</td><td>0.41</td><td>6.5</td><td>0.32</td></tr><tr><td>σ<sub>1,s</sub></td><td>4.7</td><td>1.2</td><td>5.2</td><td>0.9</td><td>5.5</td><td>1.3</td><td>5.5</td><td>1.6</td></tr><tr><td>τ<sub>1,s</sub></td><td>0.70</td><td>0.56</td><td>0.49</td><td>0.27</td><td>0.44</td><td>0.12</td><td>0.35</td><td>0.20</td></tr><tr><td></td><td colspan="8" style="text-align: center;">V<sub>2</sub> = V<sub>r</sub></td></tr><tr><td></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td><td>V<sub>2,1</sub></td><td>V<sub>2,2</sub></td></tr><tr><td>μ<sub>2,s</sub></td><td>0.13</td><td>-0.08</td><td>-0.03</td><td>0.11</td><td>-0.24</td><td>0.20</td><td>-0.31</td><td>0.20</td></tr><tr><td>σ<sub>2,s</sub></td><td>2.9</td><td>0.80</td><td>3.8</td><td>0.80</td><td>4.3</td><td>1.3</td><td>4.6</td><td>1.6</td></tr><tr><td>τ<sub>2,s</sub></td><td>0.17</td><td>-0.10</td><td>0.12</td><td>-0.04</td><td>0.07</td><td>-0.09</td><td>0.04</td><td>-0.07</td></tr></tbody></table>