The likelihood of observations \(x_1, \dots, x_T\) is
$$\mathcal{L}(\boldsymbol{\theta}) = \boldsymbol{\delta}^{(1)} \mathbf{P}(x_1) \boldsymbol{\Gamma}^{(2)} \mathbf{P}(x_2) \boldsymbol{\Gamma}^{(3)} \cdots \boldsymbol{\Gamma}^{(T)} \mathbf{P}(x_T) \mathbf{1}^\top,$$
where \(\mathbf{P}(x_t)\) is a diagonal matrix with the \(t\)-th row of allprobs on its diagonal,
\(\boldsymbol{\Gamma}^{(t)}\) is the transition probability matrix into time \(t\),
and \(\mathbf{1}\) is a row vector of ones.
The forward algorithm evaluates this product recursively and on the log-scale.
If trackID is provided, the total log-likelihood is the sum of each track's likelihood
contribution. In this case, Gamma can be
a matrix (same transition probabilities for each track),
an array of dimension c(nStates, nStates, nTracks) (track-specific transition probability matrix), or
an array of dimension c(nStates, nStates, nObs) for time-varying transition probabilities, or
an array of dimension c(nStates, nStates, L) for periodic transition probabilities, if tod is provided.
Additionally, delta can be a vector (same initial distribution for each track) or a
matrix of dimension c(nTracks, nStates) (track-specific initial distributions).
With delta = "stationary", the initial distribution of each track is computed from Gamma:
the stationary distribution of Gamma (or of its slice for each track, if track-specific),
the periodically stationary distribution at the first time point of each track, if tod is provided (see stationary_p),
the stationary distribution of the first slice of each track, for time-varying transition probabilities, i.e. at the covariate values of the first observation.
This requires Gamma to contain transition probability matrices (not, e.g., sub-stochastic matrices), which is checked without automatic differentiation.
The computed delta is reported as a matrix of dimension c(nTracks, nStates).
Note: forward should only be called once in the likelihood function, with a trackID argument for multiple tracks.
Downstream functions like lama, viterbi, stateprobs or pseudo_res rely on the quantities it reports,
and only those of the last call are kept.