Skip to contents

Find the Maximum A Posteriori Estimation

Source: R/find_MAP.R
find_MAP.Rd

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 used to find the maximum a posteriori estimate:

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

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

  • "HC" - Hill Climbing;

  • "continue" (the default for g optimized with MH or HC) - 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 when save_all_perms = TRUE. For:

  • optimizer = "MH" - normalize posterior weights over distinct visited cyclic groups, giving probabilities conditional on that set;

  • optimizer = "BF" - calculate the full posterior over all cyclic groups.

Each group contributes one weight, regardless of visit counts or the number of its generators encountered. MH continuation includes earlier MH and HC visits; keep save_all_perms = TRUE in every stage. Standalone HC does not support this option.

To examine probabilities after optimization, call get_probabilities_from_gips().

Value

An optimized object of a gips class.

Details

find_MAP() can produce a warning when:

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

  • the selected permutation has n0 > number_of_observations, or n0 > min(number_of_observations) for multiple samples. The threshold n0 already includes the additional observation required when was_mean_estimated = TRUE (for more information, see \(C_\sigma\) and n0 section in the vignette("Theory", package = "gips") or in its pkgdown page).

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 <= 10. 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.

Repeated optimization

MH and HC results can be optimized further, including with BF. If the last optimization used BF, any further call to find_MAP() is rejected, including after forget_perms(), because the whole permutation space has already been searched.

When all previous stages used MH, optimizer = "continue" resumes from the last visited state and records the shared boundary state only once. Repeated visits after rejected proposals remain in the history. An explicit optimizer = "MH" starts from the best model found so far and retains the starting state of the new stage.

Multi-sample

If g was created with gips(list(S1, S2, ...), c(n1, n2, ...)), the same optimizers search for one permutation shared by all groups. The optimized criterion is the combined unnormalized log-posterior, i.e. the sum of per-sample log weights. The default optimizer is still chosen from the common matrix size p. A sample-size warning is issued when min(n_g) < n0. Passing this check alone does not guarantee MLE availability: every projected covariance matrix must also be numerically positive definite. Use summary(g)$MLE_available on the returned object to check both conditions.

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

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
#> The permutation ():
#>  - was found after 5 posterior calculations;
#>  - is 1 times more likely than the () permutation.

g_map2 <- find_MAP(g_map, max_iter = 5, show_progress_bar = FALSE, optimizer = "continue")
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)
#> The optimized `gips` object.
#> 
#> Permutation:
#>  (1,2,3,4,5)
#> 
#> Log_posteriori:
#>  -15.1521
#> 
#> Times more likely than starting permutation:
#>  414.386
#> 
#> The p-value of Likelihood-Ratio test:
#>  0.4883
#> 
#> The number of observations:
#>  13
#> 
#> The mean in the `S` matrix was estimated.
#> Therefore, one degree of freedom was lost.
#> There are 12 degrees of freedom left.
#> 
#> n0 = number of cycles in a permutation + 1 (as the mean was estimated):
#>  2
#> 
#> The MLE estimator based on the found permutation does exist,
#> since the number of observations (13) is bigger than n0 (2),
#> and the projected covariance matrix is numerically positive definite.
#> 
#> The number of free parameters in the covariance matrix:
#>  3
#> 
#> BIC:
#>  129.0236
#> 
#> AIC:
#>  127.3288
#> 
#> --------------------------------------------------------------------------------
#> Optimization algorithm:
#>  brute_force
#> 
#> The number of log_posteriori calls:
#>  67
#> 
#> Optimization time:
#>  0.03756833 secs