# Retrieval of Similar Evolution Patterns from Satellite Image Time Series +
## 2 Materials and Methods
### Determining the Optimal Threshold Using Expectation-Maximization  
In the previous subsection, we described the process of determining the DTW distance image **DI** containing DTW similarity scores with respect to a query selected by an end-user. The following step in the automatic extraction of evolutions that have similar behavior to a given query is the optimal thresholding step. More precisely, we aim to separate the DTW distance scores between two classes, a class of _similar_ evolutions and a class of _non-similar_ spatio-temporal evolutions. Separating these two \begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline  & **0** & **1** & **0** & **2** & **1** & **3** & **0** \\ \hline
5 & 5 & 9 & 14 & 17 & 21 & 23 & 28 \\
4 & 9 & 8 & 12 & 14 & 17 & 18 & 22 \\
6 & 15 & 13 & 14 & 16 & 19 & 20 & 24 \\
3 & 18 & 15 & 16 & 15 & 16 & 16 & 19 \\
5 & 23 & 19 & 20 & 18 & 20 & 18 & 21 \\
4 & 27 & 22 & 23 & 20 & 21 & 19 & 22 \\
5 & 32 & 26 & 27 & 23 & 24 & 21 & 24 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Example of computation of the Dynamic Time Warping (DTW) distance matrix for two sequences, ‘5463545’ and ‘0102130’. classes of evolutions translates into finding the optimal threshold to be applied on the DTW similarity scores. In this sense, we formulate the problem of query-by-example SITS retrieval in terms of the Bayesian decision theory. The DTW similarity scores computed for a given query follow a multimodal distribution. This can be observed from Figure 3, where the query selected for retrieval is from an agricultural area in a long SITS (88 temporal samples). As expected, the first lobe corresponds to similar evolutions. This is argued by the fact that DTW similarity scores are smaller for similar evolutions than for non-similar evolutions. Furthermore, the histogram profile shown in Figure 3 confirms the idea of separating similar from non-similar evolutions using the Maximum Likelihood (ML) technique. In the following text, we denote the random variable associated to the computed DTW similarity scores by \(X\), whilst \(x\) represents a particular DTW score. Following the formalism initially presented in [18], let us assume that the scores are drawn from a mixture of two Gaussian distributions, \(\mathcal{N}(\mu_{s},\sigma_{s})\) and \(\mathcal{N}(\mu_{n},\sigma_{n})\). The first Gaussian distribution corresponds to evolutions that are _similar_ to the given query, whereas the second component refers to _non-similar_ evolutions. The Maximum Likelihood Estimation (MLE) framework determines the parameter values that make the observed data most likely, i.e., that maximize the likelihood function and fit the model to the data. Following the above considerations, the overall posterior probability can be decomposed as  
\[p_{\theta}(x)=\pi_{s}\mathcal{N}(x|\mu_{s},\sigma_{s})+\pi_{n}\mathcal{N}(x| \mu_{n},\sigma_{n}) \tag{6}\]  
where \(\theta=\{\pi_{s},\mu_{s},\sigma_{s},\pi_{n},\mu_{n},\sigma_{n}\}\) is the set of model parameters, whilst \(\pi_{s}\) and \(\pi_{n}\) are the mixture probabilities [19] such that \(\pi_{s}+\pi_{n}=1\). The mixture components are modeled as Gaussian densities:  
\[\mathcal{N}(x|\mu_{k},\sigma_{k})=\frac{1}{\sqrt{2\pi\sigma_{k}^{2}}}\exp \left[-\frac{(x-\mu_{k})^{2}}{2\sigma_{k}^{2}}\right] \tag{7}\]  
with \(k\in\{s,n\}\), whilst \((\mu_{s},\sigma_{s})\) and \((\mu_{n},\sigma_{n})\) are the mean and standard deviation pairs corresponding to similar and non-similar evolutions, respectively. The set of parameters \(\theta\) is determined iteratively using the Expectation-Maximization (EM) algorithm. Specifically, each iteration \(j\) is decomposed in two basic steps:  
Figure 3: Histogram over DTW similarity scores and a mixture of two Gaussian distributions fitted to these scores. 1. _Expectation step (\(E\)-step)_. Compute the log-likelihood (i.e., the logarithm of the posterior probability) with respect to the current values of the parameters \(\theta^{(t)}\): \[L(\theta^{(t)})=\ln p_{\theta^{(t)}}(x)\] (8)
2.