Historically in HMMoce we used fixed migratory and resident "speeds" (i.e. allowable diffusion) and an expectation-maximization algorithm to do parameter estimation for the behavior state switching. The EM algorithm was implemented in a function called makePar() (link) which is no longer part of the package but is retained in git version control for posterity. For our purposes, the important part of this workflow is how movement "speeds" were treated.
makePar took a numeric migr.spd input (in m/s) and the working model grid and pushed these to a function called calc.param() (link) which converted the user defined migratory speed (migr.spd in m/s) to the appropriate grid units (lines 20-21) and assigned a residency speed that was some fraction of the migratory (default was 10%). The result was a set of movement parameters that were then used to populate the gaussian kernel (with gausskern()link) which was used to do the convolution in hmm.filter() and hmm.smoother().
More recently, @paulgatti worked to implement full parameter estimation for both diffusion/movement speed parameters and behavior switching probabilities. This functionality is now wrapped in opt.params() (link) which uses any number of optimization routines (there are several options, such as optim and nlminb) to do parameter estimation by minimizing the negative log-likelihood. For simplicity, this approach is all based on the parameters being in grid units, thus the output "speed" parameters need to be converted if they're going to be intuitive as a speed parameter (such as in m/s). @marosteg has been working on an implementation of HMMoce where these intuitive speed parameters are needed and is working to implement this conversion.
There are a number of complexities in dealing with these movement kernels, including, for example, the isotropic (or lack thereof) nature of the desired kernels. Thus, we have experimented with a number of different gausskern.xx functions that currently live in the -dev branch. All are working but we have yet to do systematic testing to better understand the pros/cons of each approach. Regardless of how the gaussian kernal is populated, it will be difficult to directly (and precisely) convert these grid units to a more intuitive metric; however, an approximate conversion can be implemented after the filter and smoother steps so that the parameters can be reported in a meaningful way.
Couple things to note:
the muadv input to the gausskern.xx functions is a placeholder for adding advection to the kernels. This has not yet been implemented and is thus fixed at 0.
the first input to the gausskern.xx functions is a size parameter that is a function of the max allowed diffusion (related to migr.spd) and the grid of interest. This will be the same for the kernels for each behavior state while the sigma value used to populate the 2D gaussian distribution of the kernel is where the speed parameters actually matter
If we assume a couple movement parameters (drawn from the albacore example and rounded for simplicity) of sigma1 = 3 and sigma2 = 0.3 (resid is 10% of migr) and a size parameter of 13 on a 0.25 degree grid:
The old approach
Historically in
HMMocewe used fixed migratory and resident "speeds" (i.e. allowable diffusion) and an expectation-maximization algorithm to do parameter estimation for the behavior state switching. The EM algorithm was implemented in a function calledmakePar()(link) which is no longer part of the package but is retained in git version control for posterity. For our purposes, the important part of this workflow is how movement "speeds" were treated.makePartook a numericmigr.spdinput (in m/s) and the working model grid and pushed these to a function calledcalc.param()(link) which converted the user defined migratory speed (migr.spdin m/s) to the appropriate grid units (lines 20-21) and assigned a residency speed that was some fraction of the migratory (default was 10%). The result was a set of movement parameters that were then used to populate the gaussian kernel (withgausskern()link) which was used to do the convolution inhmm.filter()andhmm.smoother().For an example implementation of this deprecated approach see an example run script for a thresher shark dataset from the Red Sea (line ~252).
The new approach
More recently, @paulgatti worked to implement full parameter estimation for both diffusion/movement speed parameters and behavior switching probabilities. This functionality is now wrapped in
opt.params()(link) which uses any number of optimization routines (there are several options, such asoptimandnlminb) to do parameter estimation by minimizing the negative log-likelihood. For simplicity, this approach is all based on the parameters being in grid units, thus the output "speed" parameters need to be converted if they're going to be intuitive as a speed parameter (such as in m/s). @marosteg has been working on an implementation ofHMMocewhere these intuitive speed parameters are needed and is working to implement this conversion.For an example implementation of the new approach see an example run script for an albacore dataset in the N Pacific (line ~252).
To do
There are a number of complexities in dealing with these movement kernels, including, for example, the isotropic (or lack thereof) nature of the desired kernels. Thus, we have experimented with a number of different
gausskern.xxfunctions that currently live in the-devbranch. All are working but we have yet to do systematic testing to better understand the pros/cons of each approach. Regardless of how the gaussian kernal is populated, it will be difficult to directly (and precisely) convert these grid units to a more intuitive metric; however, an approximate conversion can be implemented after the filter and smoother steps so that the parameters can be reported in a meaningful way.Couple things to note:
muadvinput to thegausskern.xxfunctions is a placeholder for adding advection to the kernels. This has not yet been implemented and is thus fixed at 0.gausskern.xxfunctions is a size parameter that is a function of the max allowed diffusion (related tomigr.spd) and the grid of interest. This will be the same for the kernels for each behavior state while thesigmavalue used to populate the 2D gaussian distribution of the kernel is where the speed parameters actually matterIf we assume a couple movement parameters (drawn from the albacore example and rounded for simplicity) of sigma1 = 3 and sigma2 = 0.3 (resid is 10% of migr) and a size parameter of 13 on a 0.25 degree grid:
We get (using
fields::image.plot):