mrf2d (version 0.2.0)

fit_sa: Stochastic Approximation algorithm for MRFs on 2d lattices

Description

Estimates the parameters of a MRF by successively sampling from a parameter configuration and updating it by comparing the sufficient statistics of the sampled field and the observed field.

This method aims to find the parameter value where the gradient of the likelihood function is equal to zero.

Usage

fit_sa(
  Z,
  mrfi,
  family = "onepar",
  gamma_seq,
  init = 0,
  cycles = 5,
  refresh_each = length(gamma_seq) + 1,
  refresh_cycles = 60
)

Arguments

Z

A matrix object containing the observed MRF. NA values can be used to create a subregion of the lattice for non-rectangular data.

mrfi

A mrfi object representing the interaction structure.

family

The family of parameter restrictions to potentials. Families are: 'onepar', 'oneeach', 'absdif', 'dif' or 'free'. See mrf2d-familiy.

gamma_seq

A numeric vector with the sequence of constants used in each step \(\gamma_t\).

init

The initial value to be used in the optimization. It can be:

• A valid array of parameter values according to family.

• 0. If set to 0 an array with `0`` in all entries is created.

cycles

The number of updates to be done (for each each pixel).

refresh_each

An integer with the number of iterations taken before a complete refresh (restart from a random state). This prevents the sample from being stuck in a mode for too long. Defaults to length(gamma_seq) + 1 (no refresh happens).

refresh_cycles

An iteger indicating how many Gibbs Sampler cycles are performed when a refresh happens. Larger is usually better, but slower.

Value

A list object with the following elements:

• theta: The estimated array of potentials.

• metrics: A data.frame containing the the euclidean distance between the sufficient statics computed for Z and the current sample.

Details

The stochastic approximation method consists of, given an observed field Z, and a starting parameters configuration \(\theta_0\), successively sample a field \(Z_t\) from the current parameter configuration and estimate the direction of the gradient of the likelihood function by comparing the sufficient statistics in the current sample and the observed field.

The solution is updated by moving in the estimated direction with a predefined step size \(\gamma_t\), a new field \(Z_{t+1}\) is sampled using the new parameter configuration and \(Z_t\) as an initial value, and the process is repeated.

$$\theta_{t+1} = \theta_t - \gamma_t(T(Z_t) - T(Z)),$$ where \(T(Z)\) is the sufficient statistics for the reference field, \(T(Z_t)\) is the sufficient statistics for a field sampled from \(\theta_t\).

gamma_seq is normalized internally by diving values by length(Z), so the choice of the sequence is invariant to the lattice dimensions. Typically, a sequence like seq(from = 1, to = 0, length.out = 1000) should be used for defining a sequence with 1000 steps. Some tuning of this sequence is required.

References

wiki_samrf2d

Examples

# NOT RUN {
set.seed(2)
fit1 <- fit_sa(Z_potts, mrfi(1), family = "onepar", gamma_seq = seq(1, 0, length.out = 50))
# Estimated parameters
fit1$theta
# A visualization of estimated gradient norm over iterations.
plot(fit1$metrics)

fit_sa(Z_potts, mrfi(1), family = "oneeach", gamma_seq = seq(1, 0, length.out = 50))
# }
# NOT RUN {
# }