find_MAP {gips}

Find the Maximum A Posteriori Estimation

Description

Use one of the optimization algorithms to find the permutation that maximizes a posteriori probability based on observed data. Not all optimization algorithms will always find the MAP, but they try to find a significant value. More information can be found in the "Possible algorithms to use as optimizers" section below.

Usage

find_MAP(
  g,
  max_iter = NA,
  optimizer = NA,
  show_progress_bar = TRUE,
  save_all_perms = FALSE,
  return_probabilities = FALSE
)

Arguments

g

Object of a gips class.

max_iter

The number of iterations for an algorithm to perform. At least 2. For optimizer = "BF", it is not used; for optimizer = "MH", it has to be finite; for optimizer = "HC", it can be infinite.

optimizer

The optimizer for the search of the maximum posteriori:

• "BF" (the default for unoptimized g with ⁠perm size <= 9⁠) - Brute Force;

• "MH" (the default for unoptimized g with ⁠perm size > 10⁠) - Metropolis-Hastings;

• "HC" - Hill Climbing;

• "continue" (the default for optimized g) - The same as the g was optimized by (see Examples).

See the Possible algorithms to use as optimizers section below for more details.

show_progress_bar

A boolean. Indicate whether or not to show the progress bar:

• When max_iter is infinite, show_progress_bar has to be FALSE;

• When return_probabilities = TRUE, then shows an additional progress bar for the time when the probabilities are calculated.

save_all_perms

A boolean. TRUE indicates saving a list of all permutations visited during optimization. This can be useful sometimes but need a lot more RAM.

return_probabilities

A boolean. TRUE can only be provided only when save_all_perms = TRUE. For:

• optimizer = "MH" - use Metropolis-Hastings results to estimate posterior probabilities;

• optimizer = "BF" - use brute force results to calculate exact posterior probabilities.

These additional calculations are costly, so a second and third progress bar is shown (when show_progress_bar = TRUE).

To examine probabilities after optimization, call get_probabilities_from_gips().

Details

find_MAP can produce a warning when:

• the optimizer "hill_climbing" gets to the end of its max_iter without converging.

• the optimizer will find the permutation with smaller n0 than number_of_observations (for more information on what it means, see C_\sigma and n0 section in vignette("Theory", package = "gips") or in its pkgdown page).

Value

Returns an optimized object of a gips class.

Possible algorithms to use as optimizers

For a in-depth explanation, see in vignette("Optimizers", package = "gips") or in its pkgdown page.

For every algorithm, there are some aliases available.

• "brute_force", "BF", "full" - use the Brute Force algorithm that checks the whole permutation space of a given size. This algorithm will find the actual Maximum A Posteriori Estimation, but it is very computationally expensive for bigger spaces. We recommend Brute Force only for p <= 9. For the time the Brute Force takes on our machines, see in vignette("Optimizers", package = "gips") or in its pkgdown page.

• "Metropolis_Hastings", "MH" - use the Metropolis-Hastings algorithm; see Wikipedia. The algorithm will draw a random transposition in every iteration and consider changing the current state (permutation). When the max_iter is reached, the algorithm will return the best permutation calculated as the MAP Estimator. This implements the Second approach from references, section 4.1.2. This algorithm used in this context is a special case of the Simulated Annealing the user may be more familiar with; see Wikipedia.

• "hill_climbing", "HC" - use the hill climbing algorithm; see Wikipedia. The algorithm will check all transpositions in every iteration and go to the one with the biggest a posteriori value. The optimization ends when all neighbors will have a smaller a posteriori value. If the max_iter is reached before the end, then the warning is shown, and it is recommended to continue the optimization on the output of the find_MAP() with optimizer = "continue"; see examples. Remember that p*(p-1)/2 transpositions will be checked in every iteration. For bigger p, this may be costly.

References

Piotr Graczyk, Hideyuki Ishi, Bartosz Kołodziejek, Hélène Massam. "Model selection in the space of Gaussian models invariant by symmetry." The Annals of Statistics, 50(3) 1747-1774 June 2022. arXiv link; doi:10.1214/22-AOS2174

See Also

• gips() - The constructor of a gips class. The gips object is used as the g parameter of find_MAP().

• plot.gips() - Practical plotting function for visualizing the optimization process.

• summary.gips() - Summarize the output of optimization.

• AIC.gips(), BIC.gips() - Get the Information Criterion of the found model.

• get_probabilities_from_gips() - When find_MAP(return_probabilities = TRUE) was called, probabilities can be extracted with this function.

• log_posteriori_of_gips() - The function that the optimizers of find_MAP() tries to find the argmax of.

• forget_perms() - When the gips object was optimized with find_MAP(save_all_perms = TRUE), it will be of considerable size in RAM. forget_perms() can make such an object lighter in memory by forgetting the permutations it considered.

• vignette("Optimizers", package = "gips") or its pkgdown page - A place to learn more about the available optimizers.

• vignette("Theory", package = "gips") or its pkgdown page - A place to learn more about the math behind the gips package.

Examples

require("MASS") # for mvrnorm()

perm_size <- 5
mu <- runif(perm_size, -10, 10) # Assume we don't know the mean
sigma_matrix <- matrix(
  data = c(
    1.0, 0.8, 0.6, 0.6, 0.8,
    0.8, 1.0, 0.8, 0.6, 0.6,
    0.6, 0.8, 1.0, 0.8, 0.6,
    0.6, 0.6, 0.8, 1.0, 0.8,
    0.8, 0.6, 0.6, 0.8, 1.0
  ),
  nrow = perm_size, byrow = TRUE
) # sigma_matrix is a matrix invariant under permutation (1,2,3,4,5)
number_of_observations <- 13
Z <- MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix)
S <- cov(Z) # Assume we have to estimate the mean

g <- gips(S, number_of_observations)

g_map <- find_MAP(g, max_iter = 5, show_progress_bar = FALSE, optimizer = "Metropolis_Hastings")
g_map

g_map2 <- find_MAP(g_map, max_iter = 5, show_progress_bar = FALSE, optimizer = "continue")

if (require("graphics")) {
  plot(g_map2, type = "both", logarithmic_x = TRUE)
}

g_map_BF <- find_MAP(g, show_progress_bar = FALSE, optimizer = "brute_force")
summary(g_map_BF)