find_MAP: Find the Maximum A Posteriori Estimation
In PrzeChoj/gips: Gaussian Model Invariant by Permutation Symmetry

find_MAP

R Documentation

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 needs 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 the 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 an in-depth explanation, see in the 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 the 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; \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1214/22-AOS2174")}

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)

PrzeChoj/gips documentation built on June 12, 2025, 12:23 a.m.