Package {gips}

Type:	Package
Title:	Gaussian Model Invariant by Permutation Symmetry
Version:	1.3.1
Description:	Find the permutation symmetry group such that the covariance matrix of the given data is approximately invariant under it. Discovering such a permutation decreases the number of observations needed to fit a Gaussian model, which is of great use when it is smaller than the number of variables. Even if that is not the case, the covariance matrix found with 'gips' approximates the actual covariance with less statistical error. The methods implemented in this package are described in Graczyk et al. (2022) <doi:10.1214/22-AOS2174>. Documentation about 'gips' is provided via its website at https://przechoj.github.io/gips/ and the paper by Chojecki, Morgen, Kołodziejek (2025, <doi:10.18637/jss.v112.i07>).
License:	GPL (≥ 3)
URL:	https://github.com/PrzeChoj/gips, https://przechoj.github.io/gips/
BugReports:	https://github.com/PrzeChoj/gips/issues
Depends:	R (≥ 4.0.0)
Imports:	numbers, permutations, Rcpp, rlang (≥ 0.4.10), stringi
LinkingTo:	Rcpp
Suggests:	DAAG, dplyr, ggplot2 (≥ 3.4.0), HSAUR2, knitr, MASS (≥ 7.3-39), mvtnorm, rmarkdown, spelling, testthat (≥ 3.1.5), tibble, tidyr (≥ 1.0.0), withr
VignetteBuilder:	knitr
Config/testthat/edition:	3
Encoding:	UTF-8
Language:	en-US
Config/roxygen2/version:	8.1.0
NeedsCompilation:	yes
Packaged:	2026-11-08 20:28:03 UTC; adam
Author:	Adam Przemysław Chojecki [aut, cre], Paweł Morgen [aut], Bartosz Kołodziejek [aut]
Maintainer:	Adam Przemysław Chojecki <adam.prze.choj@gmail.com>
Repository:	CRAN
Date/Publication:	2026-11-09 00:10:12 UTC

Akaike's An Information Criterion for gips class

Description

Akaike's An Information Criterion for gips class

Usage

## S3 method for class 'gips'
AIC(object, ..., k = 2)

## S3 method for class 'gips'
BIC(object, ...)

Arguments

Value

AIC.gips() returns calculated Akaike's An Information Criterion

If the sample-size requirement for the MLE is not met, it returns NULL. If a projected covariance matrix fails the numerical positive-definiteness check, it returns Inf.

In both failure situations, shows a warning. More information can be found in the MLE availability section of logLik.gips().

BIC.gips() returns calculated Schwarz's Bayesian Information Criterion.

Functions

• BIC(gips): Schwarz's Bayesian Information Criterion

Calculation details

For more details and used formulas, see the Information Criterion - AIC and BIC section in vignette("Theory", package = "gips") or its pkgdown page.

Multi-sample

For multi-sample gips objects, AIC.gips() and BIC.gips() are based on the combined log-likelihood (sum across G groups). The penalty term uses G * dim_omega(sigma) parameters, and the BIC sample size is sum(n_g). See logLik.gips() for details.

See Also

• AIC(), BIC() - Generic functions this AIC.gips() and BIC.gips() extend.

• find_MAP() - Usually, the AIC.gips() and BIC.gips() are called on the output of find_MAP().

• logLik.gips() - Calculates the log-likelihood for the gips object. An important part of the Information Criteria.

Examples

S <- matrix(c(
  5.15, 2.05, 3.10, 1.99,
  2.05, 5.09, 2.03, 3.07,
  3.10, 2.03, 5.21, 1.97,
  1.99, 3.07, 1.97, 5.13
), nrow = 4)
g <- gips(S, 14)
g_map <- find_MAP(g, optimizer = "brute_force")

AIC(g) # 238
AIC(g_map) # 224 < 238, so g_map is better than g according to AIC
# ================================================================================
BIC(g) # 244
BIC(g_map) # 226 < 244, so g_map is better than g according to BIC

Transform the gips object to a character vector

Description

Implementation of the S3 method.

Usage

## S3 method for class 'gips'
as.character(x, ...)

Arguments

Value

An object of a character type.

See Also

• as.character.gips_perm() - The underlying gips_perm of the gips object is passed to as.character.gips_perm().

• permutations::as.character.cycle() - The underlying permutation of the gips object is passed to permutations::as.character.cycle().

Examples

A <- matrix(rnorm(4 * 4), nrow = 4)
S <- crossprod(A)
g <- gips(S, 14, perm = "(123)")
as.character(g)

Transform the gips_perm object to a character vector

Description

Implementation of the S3 method.

Usage

## S3 method for class 'gips_perm'
as.character(x, ...)

Arguments

Value

An object of a character type.

See Also

• as.character.gips() - The underlying gips_perm of the gips object is passed to as.character.gips_perm().

• permutations::as.character.cycle() - The underlying permutation of the gips object is passed to permutations::as.character.cycle().

Examples

g_perm <- gips_perm("(5,4)", 5)
as.character(g_perm)

Calculate Gamma function

Description

Calculates the integral defined in Definition 11 from references. Implements Theorem 8 from references and uses the formula (19) from references.

Usage

calculate_gamma_function(perm, lambda)

Arguments

Value

The value of the Gamma function of the colored cone (for the definition of the colored cone, see the Basic definitions section in vignette("Theory", package = "gips") or in its pkgdown page).

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

• get_structure_constants() - The function useful inside the calculate_gamma_function().

• log_posteriori_of_gips() - The function that uses the values of the gamma function.

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

Examples

id_perm <- gips_perm("()", 2)
calculate_gamma_function(id_perm, 0.5001) # 44417.54
calculate_gamma_function(id_perm, 0.50000001) # 444288280
calculate_gamma_function(id_perm, 0.500000000001) # 4.442981e+12

oldw <- getOption("warn")
options(warn = -1)
calculate_gamma_function(id_perm, 0.5) # Inf
# Integral diverges; returns Inf and warning
options(warn = oldw)

Compare the posterior probabilities of 2 permutations

Description

Check which permutation is more likely and how much more likely.

Usage

compare_posteriories_of_perms(
  perm1,
  perm2 = "()",
  S = NULL,
  number_of_observations = NULL,
  delta = 3,
  D_matrix = NULL,
  was_mean_estimated = TRUE,
  print_output = TRUE,
  digits = 3
)

compare_log_posteriories_of_perms(
  perm1,
  perm2 = "()",
  S = NULL,
  number_of_observations = NULL,
  delta = 3,
  D_matrix = NULL,
  was_mean_estimated = TRUE,
  print_output = TRUE,
  digits = 3
)

Arguments

Value

The function compare_posteriories_of_perms() returns the value of how many times the perm1 is more likely than perm2.

The function compare_log_posteriories_of_perms() returns the logarithm of how many times the perm1 is more likely than perm2.

Functions

• compare_log_posteriories_of_perms(): More stable, logarithmic version of compare_posteriories_of_perms(). The natural logarithm is used.

Multi-sample

When comparing permutations for a multi-sample gips object, or when passing multi-sample arguments directly, the comparison uses the combined unnormalized log-posterior (the sum of per-sample log weights). Two gips objects can be compared only when they are both single-sample or both multi-sample with the same data and hyperparameters.

See Also

• print.gips() - The function that prints the posterior of the optimized gips object compared to the starting permutation.

• summary.gips() - The function that calculates the posterior of the optimized gips object compared to the starting permutation.

• find_MAP() - The function that finds the permutation that maximizes log_posteriori_of_gips().

• log_posteriori_of_gips() - The function this compare_posteriories_of_perms() calls underneath.

Examples

require("MASS") # for mvrnorm()

perm_size <- 6
mu <- runif(6, -10, 10) # Assume we don't know the mean
sigma_matrix <- matrix(
  data = c(
    1.05, 0.8, 0.6, 0.4, 0.6, 0.8,
    0.8, 1.05, 0.8, 0.6, 0.4, 0.6,
    0.6, 0.8, 1.05, 0.8, 0.6, 0.4,
    0.4, 0.6, 0.8, 1.05, 0.8, 0.6,
    0.6, 0.4, 0.6, 0.8, 1.05, 0.8,
    0.8, 0.6, 0.4, 0.6, 0.8, 1.05
  ),
  nrow = perm_size, byrow = TRUE
) # sigma_matrix is a matrix invariant under permutation (1,2,3,4,5,6)
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 = 10, show_progress_bar = FALSE, optimizer = "Metropolis_Hastings")

compare_posteriories_of_perms(g_map, g, print_output = TRUE)
exp(compare_log_posteriories_of_perms(g_map, g, print_output = FALSE))

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

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).

Value

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 <= 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 the last stage used MH, optimizer = "continue" resumes from the last visited state and records the shared boundary state only once, even when earlier stages used HC. 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

• 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 visited.

• 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")
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)

Forget the permutations for gips object optimized with save_all_perms = TRUE

Description

Slim the gips object by forgetting the visited permutations from find_MAP(save_all_perms = TRUE).

Usage

forget_perms(g)

Arguments

Details

For example, perm_size = 150 and max_iter = 150000 we checked forget_perms() saves ~350 MB of RAM.

Value

The same object g as given, but without the visited permutation list.

See Also

• find_MAP() - The forget_perms() is called on the output of find_MAP(save_all_perms = TRUE).

Examples

A <- matrix(rnorm(10 * 10), nrow = 10)
S <- crossprod(A)
g <- gips(S, 13, was_mean_estimated = FALSE)
g_map <- find_MAP(g,
  max_iter = 10, optimizer = "Metropolis_Hastings",
  show_progress_bar = FALSE, save_all_perms = TRUE
)

object.size(g_map) # ~18 KB
g_map_slim <- forget_perms(g_map)
object.size(g_map_slim) # ~8 KB

Extract probabilities for gips object optimized with return_probabilities = TRUE

Description

Read probabilities previously calculated by find_MAP(return_probabilities = TRUE, save_all_perms = TRUE). This function does not calculate missing probabilities. BF gives the full posterior over cyclic groups; MH gives the posterior conditional on the distinct visited groups. See find_MAP().

Usage

get_probabilities_from_gips(g)

Arguments

Details

Saved visit-frequency estimates from older MH objects are returned unchanged with a warning; run a new optimization for weight-based results.

Value

A numeric vector of probabilities sorted in decreasing order, named by representative generators of cyclic groups. Returns NULL with a message if the optimized object has no saved probabilities.

Multi-sample

For multi-sample gips objects, these are probabilities of the shared permutation symmetry, calculated from the combined posterior used by find_MAP().

See Also

• find_MAP() - The get_probabilities_from_gips() is called on the output of find_MAP(return_probabilities = TRUE, save_all_perms = TRUE).

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

Examples

g <- gips(matrix(c(1, 0.5, 0.5, 1.3), nrow = 2), 13, was_mean_estimated = FALSE)
g_map <- find_MAP(g,
  optimizer = "BF", show_progress_bar = FALSE,
  return_probabilities = TRUE, save_all_perms = TRUE
)

get_probabilities_from_gips(g_map)

Get Structure Constants

Description

Finds constants necessary for internal calculations of integrals and eventually the posterior probability in log_posteriori_of_gips().

Usage

get_structure_constants(perm)

Arguments

Details

Uses Theorem 5 from references to calculate the constants.

Value

A list with 5 items: r, d, k, L, dim_omega - vectors of constants from Theorem 1 from references and the beginning of section 3.1. from references.

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

• calculate_gamma_function(), log_posteriori_of_gips() - The functions that rely heavily on get_structure_constants().

Examples

perm <- gips_perm("(1)(2)(3)(4,5)", 5)
get_structure_constants(perm)

The constructor of a gips class.

Description

Create a gips object. This object will contain initial data and all other information needed to find the most likely invariant permutation. It will not perform optimization; one must call the find_MAP() function to do it. See the examples below.

Usage

gips(
  S,
  number_of_observations,
  delta = 3,
  D_matrix = NULL,
  was_mean_estimated = TRUE,
  perm = ""
)

new_gips(
  list_of_gips_perm,
  S,
  number_of_observations,
  delta,
  D_matrix,
  was_mean_estimated,
  optimization_info
)

validate_gips(g)

Arguments

When `Z` is the observed data:

Value

gips() returns an object of a gips class after the safety checks.

new_gips() returns an object of a gips class without safety checks.

validate_gips() returns its argument unchanged. If the argument is not a proper element of a gips class, it produces an error.

Functions

• new_gips(): Constructor. It is only intended for low-level use.

• validate_gips(): Validator. It is only intended for low-level use.

Methods for a gips class

• summary.gips()

• plot.gips()

• print.gips()

• logLik.gips()

• AIC.gips()

• BIC.gips()

• as.character.gips()

Hyperparameters

We encourage you to try D_matrix = d * I, where I is a p \times p identity matrix and d > 0 for some different d. Smaller d favor models with more equality constraints and fewer free covariance parameters. This tendency depends on the data and delta. For fixed dimension, sample size and delta, as d tends to infinity, the posterior distribution over symmetry groups converges to a distribution independent of S.

Using S with D_matrix = d * I gives the same results as using S / d with D_matrix = I, with the other parameters unchanged.

The default for D_matrix is D_matrix = d * I, where d = mean(diag(S)). For group selection, this is equivalent to using S / d with D_matrix = I, so that the mean diagonal entry of S / d is 1.

The Bayesian model uses a Wishart-type prior on the precision matrix (the inverse covariance matrix).

The hyperparameters affect symmetry selection. For a fixed permutation, project_matrix() does not use delta or D_matrix.

For a brief introduction, see the Bayesian model selection section in vignette("Theory", package = "gips") or in its pkgdown page.

For analysis of the Hyperparameters influence, see Section 3.2. of "Learning Permutation Symmetry of a Gaussian Vector with gips in R" by gips developers Adam Chojecki, Paweł Morgen, and Bartosz Kołodziejek, Journal of Statistical Software; doi:10.18637/jss.v112.i07.

Multi-sample

The usual and recommended starting point is the single-sample call gips(S, number_of_observations).

The optional multi-sample interface is for G independent groups that are assumed to share the same permutation symmetry, while each group keeps its own covariance matrix and sample size. Pass a list of G covariance matrices as S and a numeric vector of length G as number_of_observations:

S_by_group <- list(control = S1, treatment = S2, follow_up = S3)
g <- gips(S_by_group, c(n1, n2, n3))

All matrices in S must be square covariance matrices of the same size, and their rows and columns must refer to the same variables in the same order. gips() applies one shared permutation to matching matrix indices; it does not match or reorder variables using matrix names. Names of the list elements can identify groups, while matching row and column names on the matrices can identify variables. These names are used as labels in multi-sample heatmaps in plot.gips(). D_matrix should then be a list of G positive-definite matrices (defaulting to diag(mean(diag(S_g)), p) for each group). delta can be a scalar (broadcast to all groups) or a vector of length G (one value per group); it defaults to rep(3, G). was_mean_estimated is intentionally one boolean value applied to all groups, not a per-group setting.

See Also

• stats::cov() - The S parameter, as an empirical covariance matrix, is most of the time a result of the cov() function. For more information, see Wikipedia - Estimation of covariance matrices.

• find_MAP() - The function that finds the Maximum A Posteriori (MAP) Estimator for a given gips object.

• gips_perm() - The constructor of a gips_perm class. The gips_perm object is used as the base object for the gips object. To be more precise, the base object for gips is a one-element list of a gips_perm object.

Examples

require("MASS") # for mvrnorm()

perm_size <- 5
mu <- runif(5, -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, show_progress_bar = FALSE, optimizer = "brute_force")
g_map

summary(g_map)
plot(g_map, type = "both", logarithmic_x = TRUE)

Permutation object

Description

Create permutation objects to be passed to other functions of the gips package.

Usage

gips_perm(x, size)

new_gips_perm(rearranged_cycles, size)

validate_gips_perm(g)

Arguments

Value

gips_perm() returns an object of a gips_perm class after the safety checks.

new_gips_perm() returns an object of a gips_perm class without the safety checks.

validate_gips_perm() returns its argument unchanged. If the argument is not a proper element of a gips_perm class, it produces an error.

Functions

• new_gips_perm(): Constructor. Only intended for low-level use.

• validate_gips_perm(): Validator. Only intended for low-level use.

Methods for a gips_perm class

• as.character.gips_perm()

• print.gips_perm()

See Also

• project_matrix() - gips_perm is the perm parameter of project_matrix().

• permutations::permutation() - The constructor for the x parameter.

• gips() - The constructor for the gips class uses the gips_perm object as the base object.

Examples

# All 7 following lines give the same output:
gperm <- gips_perm("(12)(45)", 5)
gperm <- gips_perm("(1,2)(4,5)", 5)
gperm <- gips_perm(as.matrix(c(2, 1, 3, 5, 4)), 5)
gperm <- gips_perm(t(as.matrix(c(2, 1, 3, 5, 4))), 5) # both matrix orientations work
gperm <- gips_perm(list(list(c(2, 1), c(4, 5))), 5)
gperm <- gips_perm(permutations::as.word(c(2, 1, 3, 5, 4)), 5)
gperm <- gips_perm(permutations::as.cycle("(1,2)(4,5)"), 5)
gperm

# note the necessity of the `size` parameter:
gperm <- gips_perm("(12)(45)", 5)
gperm <- gips_perm("(12)(45)", 7) # this one is a different permutation

try(gperm <- gips_perm("(12)(45)", 4))
# Error, `size` was set to 4, while the permutation has the element 5.

Extract the Log-Likelihood for gips class

Description

Calculates Log-Likelihood of the sample based on the gips object.

Usage

## S3 method for class 'gips'
logLik(object, ...)

Arguments

Details

The value is largest for perm = "()" when the MLEs in both the unrestricted and the selected model are available.

Value

Log-Likelihood of the sample. Object of class logLik.

Possible failure situations:

• If the sample-size requirement is not met, it returns NULL.

• If a projected covariance matrix fails the numerical positive-definiteness check, it returns -Inf.

In both failure situations, it shows a warning. Here, -Inf signals an unavailable result, not a calculated log-likelihood. More information can be found in the MLE availability section below.

MLE availability

The sample-size requirement is number_of_observations >= n0. Here, n0 is the number of permutation cycles, increased by one when was_mean_estimated = TRUE; see summary(g)$n0. Meeting this threshold does not guarantee availability for a particular S: its projection must also be numerically positive definite. See summary(g)$MLE_available for the result of these checks.

Multi-sample

For multi-sample gips objects the log-likelihood is the sum of the per-group log-likelihoods. The sample-size requirement is min(n_g) >= n0, and every projected covariance matrix must be numerically positive definite. The df attribute equals G * dim_omega(sigma) (each group has its own free covariance parameters), and nobs equals sum(n_g).

For more information, refer to C_\sigma and n0 section in vignette("Theory", package = "gips") or its pkgdown page.

Calculation details

For more details and used formulas, see the Information Criterion - AIC and BIC section in vignette("Theory", package = "gips") or its pkgdown page.

See Also

• logLik() - Generic function this logLik.gips() extends.

• find_MAP() - Usually, the logLik.gips() is called on the output of find_MAP().

• AIC.gips(), BIC.gips() - Often, one is more interested in an Information Criterion AIC or BIC.

• summary.gips() - One can get n0 by calling summary(g)$n0. To see why one may be interested in n0, see the MLE availability section above.

• project_matrix() - Project the known matrix onto the found permutations space. It is mentioned in the Calculation details section above.

Examples

S <- matrix(c(
  5.15, 2.05, 3.60, 1.99,
  2.05, 5.09, 2.03, 3.57,
  3.60, 2.03, 5.21, 1.97,
  1.99, 3.57, 1.97, 5.13
), nrow = 4)
g <- gips(S, 5)
logLik(g) # -32.67048
# For perm = "()", which is default, there is p + choose(p, 2) degrees of freedom

g_map <- find_MAP(g, optimizer = "brute_force")
logLik(g_map) # -32.6722 # this will always be smaller than `logLik(gips(S, n, perm = ""))`

g_n_too_small <- gips(S, number_of_observations = 4)
logLik(g_n_too_small) # NULL # the sample size is below n0
summary(g_n_too_small)$n0 # 5, but we set number_of_observations = 4, which is smaller

g_MAP <- find_MAP(g_n_too_small)
logLik(g_MAP) # -24.94048, this is no longer NULL
summary(g_MAP)$n0 # 2

Log of the A Posteriori of the permutation that the covariance matrix is invariant under

Description

More precisely, this is the logarithm of an unnormalized posterior probability. It is the goal function for optimization algorithms in the find_MAP() function. The perm_proposal that maximizes this function is the Maximum A Posteriori (MAP) Estimator.

Usage

log_posteriori_of_gips(g)

Arguments

Details

It is calculated using formulas (33) and (27) from references.

If Inf or NaN is reached, it produces a warning.

Value

The logarithm of an unnormalized A Posteriori.

Multi-sample

When g is a multi-sample gips object (created with a list of matrices, e.g. gips(list(S1, S2), c(n1, n2))), each sample has its own covariance matrix, while all samples share the same cyclic symmetry group. Under the package's assumptions of conditional independence of the samples and their precision-matrix priors given that group, and a uniform prior over cyclic groups,

P(\Gamma \mid S_1, \ldots, S_G) \propto \exp\left\{\sum_{j=1}^{G} \ell_j(\Gamma)\right\}.

Here, \ell_j(\Gamma) is the single-sample unnormalized log-posterior weight, computed using that sample's number_of_observations, delta, and D_matrix. The function returns the sum of these log weights, without normalizing over symmetry groups.

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

• calculate_gamma_function() - The function that calculates the value needed for log_posteriori_of_gips().

• get_structure_constants() - The function that calculates the structure constants needed for log_posteriori_of_gips().

• find_MAP() - The function that optimizes the log_posteriori_of_gips function.

• compare_posteriories_of_perms() - Uses log_posteriori_of_gips() to compare A Posteriori of two permutations.

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

Examples

# In the space with p = 2, there are only 2 permutations:
perm1 <- permutations::as.cycle("(1)(2)")
perm2 <- permutations::as.cycle("(1,2)")
S1 <- matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE)
g1 <- gips(S1, 100, perm = perm1)
g2 <- gips(S1, 100, perm = perm2)
log_posteriori_of_gips(g1) # -134.1615, this is the MAP Estimator
log_posteriori_of_gips(g2) # -138.1695

exp(log_posteriori_of_gips(g1) - log_posteriori_of_gips(g2)) # 55.0
# g1 is 55 times more likely than g2.
# This is the expected outcome because S[1,1] significantly differs from S[2,2].

compare_posteriories_of_perms(g1, g2)
# The same result, but presented in a more pleasant way

# ========================================================================

S2 <- matrix(c(1, 0.5, 0.5, 1.1), nrow = 2, byrow = TRUE)
g1 <- gips(S2, 100, perm = perm1)
g2 <- gips(S2, 100, perm = perm2)
log_posteriori_of_gips(g1) # -98.40984
log_posteriori_of_gips(g2) # -95.92039, this is the MAP Estimator

exp(log_posteriori_of_gips(g2) - log_posteriori_of_gips(g1)) # 12.05
# g2 is 12 times more likely than g1.
# This is the expected outcome because S[1,1] is very close to S[2,2].

compare_posteriories_of_perms(g2, g1)
# The same result, but presented in a more pleasant way

Plot optimized matrix or optimization gips object

Description

Plot the heatmap of the projected covariance matrix or the convergence of the optimization method. The plot depends on the type argument.

Usage

## S3 method for class 'gips'
plot(
  x,
  type = NA,
  logarithmic_y = (type != "n0"),
  logarithmic_x = FALSE,
  color = NULL,
  title_text = "Convergence plot",
  xlabel = NULL,
  ylabel = NULL,
  show_legend = TRUE,
  ylim = NULL,
  xlim = NULL,
  ...
)

Arguments

Value

An object of class ggplot.

See Also

• find_MAP() - Usually, the plot.gips() is called on the output of find_MAP().

• project_matrix() - The function used with type = "MLE".

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

Examples

require("MASS") # for mvrnorm()

perm_size <- 6
mu <- runif(6, -10, 10) # Assume we don't know the mean
sigma_matrix <- matrix(
  data = c(
    1.0, 0.8, 0.6, 0.4, 0.6, 0.8,
    0.8, 1.0, 0.8, 0.6, 0.4, 0.6,
    0.6, 0.8, 1.0, 0.8, 0.6, 0.4,
    0.4, 0.6, 0.8, 1.0, 0.8, 0.6,
    0.6, 0.4, 0.6, 0.8, 1.0, 0.8,
    0.8, 0.6, 0.4, 0.6, 0.8, 1.0
  ),
  nrow = perm_size, byrow = TRUE
) # sigma_matrix is a matrix invariant under permutation (1,2,3,4,5,6)
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)
plot(g, type = "MLE")

g_map <- find_MAP(g, max_iter = 30, show_progress_bar = FALSE, optimizer = "hill_climbing")
plot(g_map, type = "both", logarithmic_x = TRUE)

plot(g_map, type = "MLE")
# Now, the output is (most likely) different because the permutation
  # `g_map[[1]]` is (most likely) not an identity permutation.

g_map_MH <- find_MAP(g, max_iter = 30, show_progress_bar = FALSE, optimizer = "MH")
plot(g_map_MH, type = "n0", logarithmic_y = FALSE)

Prepare orthogonal matrix

Description

Calculate the orthogonal matrix U_Gamma for decomposition in Theorem 1 from references.

Usage

prepare_orthogonal_matrix(perm, perm_size = NULL, basis = NULL)

Arguments

Details

Let X be a matrix invariant under the permutation perm. Let the permutations cyclic group for it be called Gamma: \Gamma = <perm> = \{perm, perm^2, ...\}.

Then U_\Gamma is an orthogonal matrix that block-diagonalizes X.

To be more precise, the matrix t(U_Gamma) %*% X %*% U_Gamma has a block-diagonal structure, which is ensured by Theorem 1 from references.

The formula for U_Gamma can be found in Theorem 6 from references.

A nice example is demonstrated in the Block Decomposition - [1], Theorem 1 section of vignette("Theory", package="gips") or its pkgdown page.

Value

A square matrix of size perm_size by perm_size with columns from vector elements v_k^{(c)} according to Theorem 6 from references.

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

• project_matrix() - A function used in examples to show the properties of prepare_orthogonal_matrix().

• Block Decomposition - [1], Theorem 1 section of vignette("Theory", package = "gips") or its pkgdown page - A place to learn more about the math behind the gips package and see more examples of prepare_orthogonal_matrix().

Examples

gperm <- gips_perm("(1,2,3)(4,5)", 5)
U_Gamma <- prepare_orthogonal_matrix(gperm)

number_of_observations <- 10
X <- matrix(rnorm(5 * number_of_observations), number_of_observations, 5)
S <- cov(X)
X <- project_matrix(S, perm = gperm) # this matrix is invariant under gperm

block_decomposition <- t(U_Gamma) %*% X %*% U_Gamma
round(block_decomposition, 5) # the non-zeros only on diagonal and [1,2] and [2,1]

Printing gips object

Description

Printing function for a gips class.

Usage

## S3 method for class 'gips'
print(
  x,
  digits = 3,
  compare_to_original = TRUE,
  log_value = FALSE,
  oneline = FALSE,
  ...
)

Arguments

Value

An invisible NULL.

Multi-sample

For multi-sample gips objects, posterior comparisons use the combined unnormalized log-posterior (sum of per-sample log weights). Use summary(x) to display the sample sizes; see summary.gips().

See Also

• find_MAP() - The function that makes an optimized gips object out of the unoptimized one.

• compare_posteriories_of_perms() - The function that prints the compared posteriors between any two permutations, not only compared to the starting one or id.

Examples

S <- matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE)
g <- gips(S, 10, perm = "(12)")
print(g, digits = 4, oneline = TRUE)

Printing gips_perm object

Description

Printing function for a gips_perm class.

Usage

## S3 method for class 'gips_perm'
print(x, ...)

Arguments

Value

An invisible NULL.

Examples

gperm <- gips_perm("(5,4)", 5)
print(gperm)

Project matrix onto a permutation-invariant space

Description

Project S onto the space of symmetric matrices invariant under the cyclic group generated by perm. Entries that must be equal under this symmetry are replaced by their average.

Usage

project_matrix(S, perm, precomputed_equal_indices = NULL)

Arguments

Details

The permutation can be specified in advance or selected using find_MAP(). Model selection by MAP concerns the cyclic group; the subsequent projection estimates the covariance matrix within the selected model.

The projection implements Definition 3 from references.

Value

The matrix S projected on the space of symmetric matrices invariant by a cyclic group generated by perm. See Covariance estimation for its statistical interpretation.

When S is a list, returns a list of projected matrices.

Covariance estimation

Let Z contain n observations in rows. The function returns the projection of the supplied S without changing its normalization:

• If the mean is known to be zero and S = crossprod(Z) / n, the projection is the covariance MLE in the specified symmetry model, when that MLE exists.

• If S = cov(Z), the projection retains the denominator n - 1 used by stats::cov(). For the full Gaussian likelihood maximized over both an unrestricted mean and a covariance matrix with the specified symmetry, the covariance MLE would be (n - 1) / n * project_matrix(S, perm), when it exists.

For S = cov(Z) with was_mean_estimated = TRUE, gips uses the projection without this rescaling, and logLik.gips() uses n - 1 in its likelihood calculation. This is the package's convention based on centered data, distinct from the full Gaussian likelihood described above. See the Estimated mean subsection in vignette("Theory", package = "gips").

MLE availability

project_matrix() does not receive the sample size and does not check whether an MLE is available. It can return a singular projected matrix. For a corresponding gips object g, summary.gips() reports the sample-size threshold in summary(g)$n0. The field summary(g)$MLE_available combines the sample-size check with a numerical positive-definiteness check of the projection, using the package's convention described above.

Multi-sample

When S is a list of matrices (as stored in a multi-sample gips object), project_matrix() applies the same permutation projection to each element and returns a list of projected matrices of the same length.

See Also

• Wikipedia - Estimation of covariance matrices

• Project Matrix - Equation (6) section of vignette("Theory", package = "gips") or its pkgdown page - A place to learn more about the math behind the gips package and see more examples of project_matrix().

• find_MAP() - Searches for a cyclic group with maximum posterior weight. Heuristic optimizers do not guarantee finding the maximum. The returned permutation can be used to project S onto the selected symmetry model (see examples).

• gips_perm() - Constructor for the perm parameter.

• plot.gips() - For plot(g, type = "MLE"), the project_matrix() is called (see examples).

• summary.gips() - Reports n0 and MLE_available for a gips object; see the MLE availability section.

Examples

p <- 6
my_perm <- "(14)(23)" # permutation (1,4)(2,3)(5)(6)
number_of_observations <- 10
X <- matrix(rnorm(p * number_of_observations), number_of_observations, p)
S <- cov(X)
projected_S <- project_matrix(S, perm = my_perm)
projected_S
# The value in [1,1] is the same as in [4,4]; also, [2,2] and [3,3];
  # also [1,2] and [3,4]; also, [1,5] and [4,5]; and so on

# Plot the projected matrix:
g <- gips(S, number_of_observations, perm = my_perm)
plot(g, type = "MLE")

# Search for the MAP symmetry group:
g_MAP <- find_MAP(g, max_iter = 10, show_progress_bar = FALSE, optimizer = "Metropolis_Hastings")
# Then, project S onto the selected model:
S_MAP <- project_matrix(attr(g, "S"), perm = g_MAP)
S_MAP
plot(g_MAP, type = "heatmap")

Summarizing the gips object

Description

summary method for gips class.

Usage

## S3 method for class 'gips'
summary(object, ...)

## S3 method for class 'summary.gips'
print(x, ...)

Arguments

Value

The function summary.gips() computes and returns a list of summary statistics of the given gips object. Those are:

• For unoptimized gips object:

  1. optimized - FALSE.

  2. start_permutation - the permutation this gips represents.

  3. start_permutation_log_posteriori - the log of the A Posteriori value the start permutation has.

  4. times_more_likely_than_id - how many times more likely the start_permutation is over the identity permutation, ⁠()⁠. It can be less than 1, meaning the identity permutation is more likely. Remember that this number can become very large and overflow to Inf or small and underflow to 0.

  5. log_times_more_likely_than_id - log of times_more_likely_than_id.

  6. likelihood_ratio_test_statistics, likelihood_ratio_test_p_value - statistics and p-value of Likelihood Ratio test, where the H_0 is that the data was drawn from the normal distribution with Covariance matrix invariant under the given permutation. The p-value is calculated from the asymptotic distribution. For interpretation, see the Likelihood-ratio test section in vignette("Theory", package = "gips"). Comparison with the unrestricted model requires n \ge p, or n \ge p + 1 when the mean was estimated.

  7. n0 - the sample-size threshold: the number of permutation cycles, increased by one when was_mean_estimated = TRUE. The projection must also be numerically positive definite; see MLE_available below and C_\sigma and n0 section in vignette("Theory", package = "gips") or in its pkgdown page.

  8. S_matrix - the underlying matrix. This matrix will be used in calculations of the posterior value in log_posteriori_of_gips().

  9. number_of_observations - the number of observations that were observed for the S_matrix to be calculated. This value will be used in calculations of the posterior value in log_posteriori_of_gips().

  10. was_mean_estimated - given by the user while creating the gips object:

      • TRUE means the S parameter was the output of stats::cov() function;

      • FALSE means the S parameter was calculated with S = crossprod(X) / number_of_observations.

  11. delta, D_matrix - the hyperparameters of the Bayesian method. See the Hyperparameters section of gips() documentation.

  12. n_parameters - number of free parameters in the covariance matrix.

  13. AIC, BIC - output of AIC.gips() and BIC.gips() functions.

• For optimized gips object:

  1. optimized - TRUE.

  2. found_permutation - the permutation this gips represents. The visited permutation with the biggest A Posteriori value.

  3. found_permutation_log_posteriori - the log of the A Posteriori value the found permutation has.

  4. start_permutation - the original permutation this gips represented before optimization. It is the first visited permutation.

  5. start_permutation_log_posteriori - the log of the A Posteriori value the start permutation has.

  6. times_more_likely_than_start - how many times more likely the found_permutation is over the start_permutation. It cannot be a number less than 1. Remember that this number can become large and overflow to Inf.

  7. log_times_more_likely_than_start - log of times_more_likely_than_start.

  8. likelihood_ratio_test_statistics, likelihood_ratio_test_p_value - statistics and p-value of Likelihood Ratio test, where the H_0 is that the data was drawn from the normal distribution with Covariance matrix invariant under found_permutation. The p-value is calculated from the asymptotic distribution. For interpretation, see the Likelihood-ratio test section in vignette("Theory", package = "gips"). Comparison with the unrestricted model requires n \ge p, or n \ge p + 1 when the mean was estimated.

  9. n0 - the sample-size threshold: the number of permutation cycles, increased by one when was_mean_estimated = TRUE. The projection must also be numerically positive definite; see MLE_available below and C_\sigma and n0 section in vignette("Theory", package = "gips") or in its pkgdown page.

  10. S_matrix - the underlying matrix. This matrix will be used in calculations of the posterior value in log_posteriori_of_gips().

  11. number_of_observations - the number of observations that were observed for the S_matrix to be calculated. This value will be used in calculations of the posterior value in log_posteriori_of_gips().

  12. was_mean_estimated - given by the user while creating the gips object:

      • TRUE means the S parameter was output of the stats::cov() function;

      • FALSE means the S parameter was calculated with S = crossprod(X) / number_of_observations.

  13. delta, D_matrix - the hyperparameters of the Bayesian method. See the Hyperparameters section of gips() documentation.

  14. n_parameters - number of free parameters in the covariance matrix.

  15. AIC, BIC - output of AIC.gips() and BIC.gips() functions.

  16. optimization_algorithm_used - all used optimization algorithms in order (one could start optimization with "MH", and then do an "HC").

  17. did_converge - a boolean, did the last used algorithm converge.

  18. number_of_log_posteriori_calls - the number of actual posterior evaluations during the recorded optimization. For multiple samples, one evaluation computes the joint objective.

  19. whole_optimization_time - how long was the optimization process; the sum of all optimization times (when there were multiple).

  20. log_posteriori_calls_after_best - the number of posterior evaluations after the best log-posterior value was first obtained. If this value is small, consider running find_MAP() again with optimizer = "continue". For optimizer = "BF", it is NULL.

  21. acceptance_rate - the fraction of eligible MH proposals accepted. Counts are pooled across MH stages, including in mixed histories; HC contributes to neither count. The value is NULL for histories without MH stages.

Both summaries also contain MLE_available: whether the sample-size requirement is met and all projected covariance matrices pass the numerical positive-definiteness check used by logLik().

The function print.summary.gips() returns an invisible NULL.

Methods (by generic)

• print(summary.gips): Printing method for class summary.gips. Prints every interesting information in a form pleasant for humans.

Multi-sample

For multi-sample gips objects the summary includes the number of groups G, the full vector of per-group sample sizes, and the total n_parameters (which equals G * dim_omega(sigma)). The sample-size check is min(n_g) >= n0, and every projected covariance matrix must be numerically positive definite for MLE_available to be TRUE. The Likelihood-Ratio test is not available for multi-sample objects and is reported as NULL.

See Also

• find_MAP() - Usually, the summary.gips() is called on the output of find_MAP().

• log_posteriori_of_gips() - Calculate the likelihood of a permutation.

• AIC.gips(), BIC.gips() - Calculate Akaike's or Bayesian Information Criterion

• project_matrix() - Project the known matrix on the found permutations space.

Examples

require("MASS") # for mvrnorm()

perm_size <- 6
mu <- runif(6, -10, 10) # Assume we don't know the mean
sigma_matrix <- matrix(
  data = c(
    1.1, 0.8, 0.6, 0.4, 0.6, 0.8,
    0.8, 1.1, 0.8, 0.6, 0.4, 0.6,
    0.6, 0.8, 1.1, 0.8, 0.6, 0.4,
    0.4, 0.6, 0.8, 1.1, 0.8, 0.6,
    0.6, 0.4, 0.6, 0.8, 1.1, 0.8,
    0.8, 0.6, 0.4, 0.6, 0.8, 1.1
  ),
  nrow = perm_size, byrow = TRUE
) # sigma_matrix is a matrix invariant under permutation (1,2,3,4,5,6)
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)
unclass(summary(g))

g_map <- find_MAP(g, max_iter = 10, show_progress_bar = FALSE, optimizer = "Metropolis_Hastings")
unclass(summary(g_map))

g_map2 <- find_MAP(g, max_iter = 10, show_progress_bar = FALSE, optimizer = "hill_climbing")
summary(g_map2)
# ================================================================================
S <- matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE)
g <- gips(S, 10)
print(summary(g))