| Title: | Training DA Models Utilizing 'gips' |
|---|---|
| Description: | Extends classical linear and quadratic discriminant analysis by incorporating permutation group symmetries into covariance matrix estimation. The package leverages methodology from the 'gips' framework to identify and impose permutation structures that act as a form of regularization, improving stability and interpretability in settings with symmetric or exchangeable features. Several discriminant analysis variants are provided, including pooled and class-specific covariance models, as well as multi-class extensions with shared or independent symmetry structures. For more details about 'gips' methodology see and Graczyk et al. (2022) <doi:10.1214/22-AOS2174> and Chojecki, Morgen, Kołodziejek (2025, <doi:10.18637/jss.v112.i07>). |
| Authors: | Antoni Zbigniew Kingston [aut], Norbert Maksymilian Frydrysiak [aut, cre], Adam Przemysław Chojecki [ctb] (ORCID: <https://orcid.org/0009-0008-2902-4096>) |
| Maintainer: | Norbert Maksymilian Frydrysiak <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.2 |
| Built: | 2026-06-08 06:34:53 UTC |
| Source: | https://github.com/antonikingston/gipsda |
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.
find_MAP( g, max_iter = NA, optimizer = NA, show_progress_bar = TRUE, save_all_perms = FALSE, return_probabilities = FALSE )find_MAP( g, max_iter = NA, optimizer = NA, show_progress_bar = TRUE, save_all_perms = FALSE, return_probabilities = FALSE )
g |
Object of a |
max_iter |
The number of iterations for an algorithm to perform.
At least 2. For |
optimizer |
The optimizer for the search of the maximum posteriori:
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:
|
save_all_perms |
A boolean. |
return_probabilities |
A boolean.
These additional calculations are costly, so a second and third
progress bar is shown (when To examine probabilities after optimization,
call |
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
Returns an optimized object of a gipsmult class.
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.
"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 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.
gipsmult() - The constructor of a gipsmult class.
The gipsmult object is used as the g parameter of find_MAP().
plot.gipsmult() - Practical plotting function for
visualizing the optimization process.
get_probabilities_from_gipsmult() - When
find_MAP(return_probabilities = TRUE) was called,
probabilities can be extracted with this function.
log_posteriori_of_gipsmult() - The function that the optimizers
of find_MAP() tries to find the argmax of.
require("MASS") # for mvrnorm() perm_size <- 6 mu1 <- runif(6, -10, 10) mu2 <- runif(6, -10, 10) # Assume we don't know the means sigma1 <- 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 ) sigma2 <- matrix( data = c( 1.0, 0.5, 0.2, 0.0, 0.2, 0.5, 0.5, 1.0, 0.5, 0.2, 0.0, 0.2, 0.2, 0.5, 1.0, 0.5, 0.2, 0.0, 0.0, 0.2, 0.5, 1.0, 0.5, 0.2, 0.2, 0.0, 0.2, 0.5, 1.0, 0.5, 0.5, 0.2, 0.0, 0.2, 0.5, 1.0 ), nrow = perm_size, byrow = TRUE ) # sigma1 and sigma2 are matrices invariant under permutation (1,2,3,4,5,6) numbers_of_observations <- c(21, 37) Z1 <- MASS::mvrnorm(numbers_of_observations[1], mu = mu1, Sigma = sigma1) Z2 <- MASS::mvrnorm(numbers_of_observations[2], mu = mu2, Sigma = sigma2) S1 <- cov(Z1) S2 <- cov(Z2) # Assume we have to estimate the mean g <- gipsmult(list(S1, S2), numbers_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")require("MASS") # for mvrnorm() perm_size <- 6 mu1 <- runif(6, -10, 10) mu2 <- runif(6, -10, 10) # Assume we don't know the means sigma1 <- 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 ) sigma2 <- matrix( data = c( 1.0, 0.5, 0.2, 0.0, 0.2, 0.5, 0.5, 1.0, 0.5, 0.2, 0.0, 0.2, 0.2, 0.5, 1.0, 0.5, 0.2, 0.0, 0.0, 0.2, 0.5, 1.0, 0.5, 0.2, 0.2, 0.0, 0.2, 0.5, 1.0, 0.5, 0.5, 0.2, 0.0, 0.2, 0.5, 1.0 ), nrow = perm_size, byrow = TRUE ) # sigma1 and sigma2 are matrices invariant under permutation (1,2,3,4,5,6) numbers_of_observations <- c(21, 37) Z1 <- MASS::mvrnorm(numbers_of_observations[1], mu = mu1, Sigma = sigma1) Z2 <- MASS::mvrnorm(numbers_of_observations[2], mu = mu2, Sigma = sigma2) S1 <- cov(Z1) S2 <- cov(Z2) # Assume we have to estimate the mean g <- gipsmult(list(S1, S2), numbers_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")
gipsmult object optimized with return_probabilities = TRUE
After the gipsmult object was optimized with
the find_MAP(return_probabilities = TRUE) function, then
those calculated probabilities can be extracted with this function.
get_probabilities_from_gipsmult(g)get_probabilities_from_gipsmult(g)
g |
An object of class |
Returns a numeric vector, calculated values of probabilities.
Names contain permutations this probabilities represent.
For gipsmult object optimized with find_MAP(return_probabilities = FALSE),
it returns a NULL object.
It is sorted according to the probability.
find_MAP() - The get_probabilities_from_gipsmult()
is called on the output of
find_MAP(return_probabilities = TRUE, save_all_perms = TRUE).
Ss <- list( matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE), matrix(c(2, 1, 3, 7), nrow = 2, byrow = TRUE) ) noo <- c(10, 13) g <- gipsmult(Ss, noo) g_map <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE, return_probabilities = TRUE, save_all_perms = TRUE ) get_probabilities_from_gipsmult(g_map)Ss <- list( matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE), matrix(c(2, 1, 3, 7), nrow = 2, byrow = TRUE) ) noo <- c(10, 13) g <- gipsmult(Ss, noo) g_map <- find_MAP(g, optimizer = "BF", show_progress_bar = FALSE, return_probabilities = TRUE, save_all_perms = TRUE ) get_probabilities_from_gipsmult(g_map)
Linear discriminant analysis (LDA) using covariance matrices projected via the gips framework to enforce permutation symmetry and improve numerical stability.
gipslda(x, ...) ## S3 method for class 'formula' gipslda(formula, data, ..., subset, na.action) ## Default S3 method: gipslda(x, grouping, prior = proportions, tol = 1e-4, weighted_avg = FALSE, MAP = TRUE, optimizer = NULL, max_iter = NULL, ...) ## S3 method for class 'data.frame' gipslda(x, ...) ## S3 method for class 'matrix' gipslda(x, grouping, ..., subset, na.action)gipslda(x, ...) ## S3 method for class 'formula' gipslda(formula, data, ..., subset, na.action) ## Default S3 method: gipslda(x, grouping, prior = proportions, tol = 1e-4, weighted_avg = FALSE, MAP = TRUE, optimizer = NULL, max_iter = NULL, ...) ## S3 method for class 'data.frame' gipslda(x, ...) ## S3 method for class 'matrix' gipslda(x, grouping, ..., subset, na.action)
x |
(required if no formula is given as the principal argument) a matrix or data frame or Matrix containing the explanatory variables. |
... |
Arguments passed to or from other methods. |
formula |
A formula of the form |
data |
An optional data frame, list or environment from which variables
specified in |
grouping |
(required if no formula principal argument is given) a factor specifying the class for each observation. |
prior |
The prior probabilities of class membership. If unspecified, the class proportions for the training set are used. |
tol |
A tolerance to decide if a matrix is singular; variables whose
variance is less than |
subset |
An index vector specifying the cases to be used in the training sample. (NOTE: must be named.) |
na.action |
A function specifying the action for |
MAP |
Logical; whether to compute a Maximum A Posteriori gips projection of the covariance matrix. |
optimizer |
Character; optimization method used by gips
(e.g. |
max_iter |
Maximum number of iterations for the optimizer. |
weighted_avg |
Logical; Whether to compute scatter from all classes at once or to compute them within classes and compute the main one as average weighted by class proportions. |
This function is a minor modification of lda, replacing
the classical sample covariance estimators by projected covariance matrices
obtained using project_covs().
Unlike classical LDA, the within-class covariance matrix is first projected onto a permutation-invariant structure using the gips framework. This can stabilize covariance estimation in high dimensions or when symmetry assumptions are justified.
The choice of optimizer and MAP estimation affects both the covariance estimate and the resulting discriminant directions.
See Chojecki et al. (2025) for theoretical background.
An object of class "gipslda" containing:
prior: prior class probabilities
counts: number of observations per class
means: group means
scaling: linear discriminant coefficients
svd: singular values of the between-class scatter
N: number of observations
optimization_info: information about the gips optimization
call: matched call
This function is inspired by lda but is not
a drop-in replacement. The covariance estimator, optimization
procedure, and returned object differ substantially.
Chojecki, A., et al. (2025). Learning Permutation Symmetry of a Gaussian Vector with gips in R. Journal of Statistical Software, 112(7), 1–38. doi:10.18637/jss.v112.i07
Iris <- data.frame(rbind(iris3[, , 1], iris3[, , 2], iris3[, , 3]), Sp = rep(c("s", "c", "v"), rep(50, 3)) ) train <- sample(1:150, 75) z <- gipslda(Sp ~ ., Iris, prior = c(1, 1, 1) / 3, subset = train) predict(z, Iris[-train, ])$class (z1 <- update(z, . ~ . - Petal.W.))Iris <- data.frame(rbind(iris3[, , 1], iris3[, , 2], iris3[, , 3]), Sp = rep(c("s", "c", "v"), rep(50, 3)) ) train <- sample(1:150, 75) z <- gipslda(Sp ~ ., Iris, prior = c(1, 1, 1) / 3, subset = train) predict(z, Iris[-train, ])$class (z1 <- update(z, . ~ . - Petal.W.))
gipsmult class.Create a gipsmult 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.
gipsmult( Ss, numbers_of_observations, delta = 3, D_matrices = NULL, was_mean_estimated = TRUE, perm = "" ) new_gipsmult( list_of_gips_perm, Ss, numbers_of_observations, delta, D_matrices, was_mean_estimated, optimization_info )gipsmult( Ss, numbers_of_observations, delta = 3, D_matrices = NULL, was_mean_estimated = TRUE, perm = "" ) new_gipsmult( list_of_gips_perm, Ss, numbers_of_observations, delta, D_matrices, was_mean_estimated, optimization_info )
Ss |
A list of matrices; empirical covariance matrices.
When
|
numbers_of_observations |
Numbers of data points
that |
delta |
A number, hyper-parameter of a Bayesian model. It has to be strictly bigger than 1. See the Hyperparameters section below. |
D_matrices |
A list of symmetric, positive-definite matrices of the same size as matrices in |
was_mean_estimated |
A boolean.
|
perm |
An optional permutation to be the base for the |
list_of_gips_perm |
A list with a single element of
a |
optimization_info |
For internal use only. |
gipsmult() returns an object of
a gipsmult class after the safety checks.
new_gipsmult() returns an object of
a gipsmult class without the safety checks.
new_gipsmult(): Constructor. It is only intended for low-level use.
gipsmult classWe encourage the user to try D_matrix = d * I, where I is an identity matrix of a size
p x p and d > 0 for some different d.
When d is small compared to the data (e.g., d = 0.1 * mean(diag(S))),
bigger structures will be found.
When d is big compared to the data (e.g., d = 100 * mean(diag(S))),
the posterior distribution does not depend on the data.
Taking D_matrix = d * I is equivalent to setting S <- S / d.
The default for D_matrix is D_matrix = d * I, where
d = mean(diag(S)), which is equivalent to modifying S
so that the mean value on the diagonal is 1.
In the Bayesian model, the prior distribution for the covariance matrix is a generalized case of Wishart distribution.
stats::cov() – The Ss parameter, as a list of empirical covariance matrices,
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 gipsmult object.
gips::gips_perm() – The constructor of a gips_perm class.
The gips_perm object is used as the base object for
the gipsmult object.
perm_size <- 5 numbers_of_observations <- c(15, 18, 19) Sigma <- diag(rep(1, perm_size)) n_matrices <- 3 df <- 20 Ss <- rWishart(n = n_matrices, df = df, Sigma = Sigma) Ss <- lapply(1:n_matrices, function(x) Ss[, , x]) g <- gipsmult(Ss, numbers_of_observations) g_map <- find_MAP(g, show_progress_bar = FALSE, optimizer = "brute_force") g_map print(g_map) if (require("graphics")) { plot(g_map, type = "MLE", logarithmic_x = TRUE) }perm_size <- 5 numbers_of_observations <- c(15, 18, 19) Sigma <- diag(rep(1, perm_size)) n_matrices <- 3 df <- 20 Ss <- rWishart(n = n_matrices, df = df, Sigma = Sigma) Ss <- lapply(1:n_matrices, function(x) Ss[, , x]) g <- gipsmult(Ss, numbers_of_observations) g_map <- find_MAP(g, show_progress_bar = FALSE, optimizer = "brute_force") g_map print(g_map) if (require("graphics")) { plot(g_map, type = "MLE", logarithmic_x = TRUE) }
Quadratic Discriminant Analysis (QDA) in which each class covariance matrix is projected using the gipsmult framework, allowing for structured permutation symmetry across multiple covariance matrices.
gipsmultqda(x, ...) ## S3 method for class 'formula' gipsmultqda(formula, data, ..., subset, na.action) ## Default S3 method: gipsmultqda(x, grouping, prior = proportions, nu = 5, MAP = TRUE, optimizer = NULL, max_iter = NULL, ...) ## S3 method for class 'data.frame' gipsmultqda(x, ...) ## S3 method for class 'matrix' gipsmultqda(x, grouping, ..., subset, na.action)gipsmultqda(x, ...) ## S3 method for class 'formula' gipsmultqda(formula, data, ..., subset, na.action) ## Default S3 method: gipsmultqda(x, grouping, prior = proportions, nu = 5, MAP = TRUE, optimizer = NULL, max_iter = NULL, ...) ## S3 method for class 'data.frame' gipsmultqda(x, ...) ## S3 method for class 'matrix' gipsmultqda(x, grouping, ..., subset, na.action)
x |
(required if no formula is given as the principal argument) a matrix or data frame containing the explanatory variables. |
... |
Arguments passed to or from other methods. |
formula |
A formula of the form |
data |
An optional data frame, list or environment from which variables
specified in |
grouping |
A factor specifying the class for each observation. |
prior |
Prior probabilities of class membership. Must sum to one. |
nu |
Degrees of freedom parameter used internally during covariance projection. |
MAP |
Logical; if |
optimizer |
Character string specifying the optimization method used
for covariance projection. If |
max_iter |
Maximum number of iterations for stochastic optimizers. |
subset |
An index vector specifying the cases to be used in the training sample. (NOTE: must be named.) |
na.action |
A function specifying the action to be taken if |
This function is a modification of qda in which the
class-specific covariance matrices are jointly projected to improve
numerical stability and exploit shared symmetry assumptions.
In contrast to classical QDA, which estimates each class covariance matrix
independently, gipsmultqda performs a joint projection of all class
covariance matrices using the gipsmult framework. This allows the
incorporation of shared permutation symmetries and can improve classification
performance in high-dimensional or small-sample regimes.
Several classification rules are available via
predict.gipsmultqda, including plug-in, predictive, debiased,
and leave-one-out cross-validation.
An object of class "gipsmultqda" containing:
prior: prior probabilities of the groups
counts: number of observations per group
means: group means
scaling: array of group-specific scaling matrices derived
from the projected covariance matrices
ldet: log-determinants of the projected covariance matrices
lev: class labels
N: total number of observations
optimization_info: information returned by the covariance
projection optimizer
call: the matched call
This function is not a drop-in replacement for qda.
The covariance estimation, returned object, and classification rules
differ substantially.
The theoretical background and details of the covariance projection are
documented in the gipsmult package.
qda, predict.gipsmultqda,
gipsqda, gipslda
tr <- sample(1:50, 25) train <- rbind(iris3[tr, , 1], iris3[tr, , 2], iris3[tr, , 3]) test <- rbind(iris3[-tr, , 1], iris3[-tr, , 2], iris3[-tr, , 3]) cl <- factor(c(rep("s", 25), rep("c", 25), rep("v", 25))) z <- gipsmultqda(train, cl) predict(z, test)$classtr <- sample(1:50, 25) train <- rbind(iris3[tr, , 1], iris3[tr, , 2], iris3[tr, , 3]) test <- rbind(iris3[-tr, , 1], iris3[-tr, , 2], iris3[-tr, , 3]) cl <- factor(c(rep("s", 25), rep("c", 25), rep("v", 25))) z <- gipsmultqda(train, cl) predict(z, test)$class
Quadratic discriminant analysis (QDA) using covariance matrices projected via the gips framework to enforce permutation symmetry and improve numerical stability.
gipsqda(x, ...) ## S3 method for class 'formula' gipsqda(formula, data, ..., subset, na.action) ## Default S3 method: gipsqda(x, grouping, prior = proportions, nu = 5, MAP = TRUE, optimizer = NULL, max_iter = NULL, ...) ## S3 method for class 'data.frame' gipsqda(x, ...) ## S3 method for class 'matrix' gipsqda(x, grouping, ..., subset, na.action)gipsqda(x, ...) ## S3 method for class 'formula' gipsqda(formula, data, ..., subset, na.action) ## Default S3 method: gipsqda(x, grouping, prior = proportions, nu = 5, MAP = TRUE, optimizer = NULL, max_iter = NULL, ...) ## S3 method for class 'data.frame' gipsqda(x, ...) ## S3 method for class 'matrix' gipsqda(x, grouping, ..., subset, na.action)
x |
(required if no formula is given as the principal argument) a matrix or data frame containing the explanatory variables. |
... |
Arguments passed to or from other methods. |
formula |
A formula of the form |
data |
An optional data frame, list or environment from which variables
specified in |
grouping |
(required if no formula is given) a factor specifying the class for each observation. |
prior |
The prior probabilities of class membership. Must sum to one and have length equal to the number of groups. |
nu |
Degrees of freedom parameter used internally by covariance projection. |
MAP |
Logical; if |
optimizer |
Character string specifying the optimization method used
for covariance projection. If |
max_iter |
Maximum number of iterations for stochastic optimizers. |
subset |
An index vector specifying the cases to be used in the training sample. (NOTE: must be named.) |
na.action |
A function specifying the action to be taken if |
This function is a minor modification of qda, replacing
the classical sample covariance estimators by projected covariance matrices
obtained using project_covs().
Quadratic discriminant analysis models each class with its own covariance
matrix. In gipsqda, these covariance matrices are projected using the
gips framework, which enforces permutation symmetry and mitigates
singularity and overfitting in high-dimensional or small-sample settings.
Classification can be performed using plug-in, predictive, debiased,
or leave-one-out cross-validation rules via predict.gipsqda.
An object of class "gipsqda" containing the following components:
prior: prior probabilities of the groups
counts: number of observations in each group
means: group means
scaling: group-specific scaling matrices derived from the
projected covariance matrices
ldet: log-determinants of the projected covariance matrices
lev: class labels
N: total number of observations
optimization_info: information returned by the covariance
projection optimizer
call: the matched call
The function may be called with either a formula interface or with a matrix
and grouping factor. Arguments subset and na.action, if used,
must be named.
Chojecki, A., et al. (2025). Learning Permutation Symmetry of a Gaussian Vector with gips in R. Journal of Statistical Software, 112(7), 1–38. doi:10.18637/jss.v112.i07
Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with S. Fourth edition. Springer.
qda, predict.gipsqda,
gipslda, lda
tr <- sample(1:50, 25) train <- rbind(iris3[tr, , 1], iris3[tr, , 2], iris3[tr, , 3]) test <- rbind(iris3[-tr, , 1], iris3[-tr, , 2], iris3[-tr, , 3]) cl <- factor(c(rep("s", 25), rep("c", 25), rep("v", 25))) z <- gipsqda(train, cl) predict(z, test)$classtr <- sample(1:50, 25) train <- rbind(iris3[tr, , 1], iris3[tr, , 2], iris3[tr, , 3]) test <- rbind(iris3[-tr, , 1], iris3[-tr, , 2], iris3[-tr, , 3]) cl <- factor(c(rep("s", 25), rep("c", 25), rep("v", 25))) z <- gipsqda(train, cl) predict(z, test)$class
More precisely, it 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.
log_posteriori_of_gipsmult(g)log_posteriori_of_gipsmult(g)
g |
An object of a |
It is calculated using formulas (33) and (27) from references.
If Inf or NaN is reached, it produces a warning.
Returns a value of the logarithm of an unnormalized A Posteriori.
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
find_MAP() - The function that optimizes
the log_posteriori_of_gips function.
gips::compare_posteriories_of_perms() - Uses log_posteriori_of_gips()
to compare a posteriori of two permutations.
# In the space with p = 2, there is 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) S2 <- matrix(c(2, 1, 3, 7), nrow = 2, byrow = TRUE) g1 <- gipsmult(list(S1, S2), c(100, 100), perm = perm1) g2 <- gipsmult(list(S1, S2), c(100, 100), perm = perm2) log_posteriori_of_gipsmult(g1) # -354.4394, this is the MAP Estimator log_posteriori_of_gipsmult(g2) # -380.0079 exp(log_posteriori_of_gipsmult(g1) - log_posteriori_of_gipsmult(g2)) # 127131902082 # g1 is 127131902082 times more likely than g2. # This is the expected outcome because S1[1,1] and S2[1,1] # differ significantly from S1[2,2] and S2[2,2] respectively. # ======================================================================== S3 <- matrix(c(1, 0.5, 0.5, 1.1), nrow = 2, byrow = TRUE) S4 <- matrix(c(2, 1, 3, 2.137), nrow = 2, byrow = TRUE) g1 <- gipsmult(list(S3, S4), c(100, 100), perm = perm1) g2 <- gipsmult(list(S3, S4), c(100, 100), perm = perm2) log_posteriori_of_gipsmult(g1) # -148.6485 log_posteriori_of_gipsmult(g2) # -145.3019, this is the MAP Estimator exp(log_posteriori_of_gipsmult(g2) - log_posteriori_of_gipsmult(g1)) # 28.406 # g2 is 28.406 times more likely than g1. # This is the expected outcome because S1[1,1] and S2[1,1] # are very close to S1[2,2] and S2[2,2] respectively.# In the space with p = 2, there is 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) S2 <- matrix(c(2, 1, 3, 7), nrow = 2, byrow = TRUE) g1 <- gipsmult(list(S1, S2), c(100, 100), perm = perm1) g2 <- gipsmult(list(S1, S2), c(100, 100), perm = perm2) log_posteriori_of_gipsmult(g1) # -354.4394, this is the MAP Estimator log_posteriori_of_gipsmult(g2) # -380.0079 exp(log_posteriori_of_gipsmult(g1) - log_posteriori_of_gipsmult(g2)) # 127131902082 # g1 is 127131902082 times more likely than g2. # This is the expected outcome because S1[1,1] and S2[1,1] # differ significantly from S1[2,2] and S2[2,2] respectively. # ======================================================================== S3 <- matrix(c(1, 0.5, 0.5, 1.1), nrow = 2, byrow = TRUE) S4 <- matrix(c(2, 1, 3, 2.137), nrow = 2, byrow = TRUE) g1 <- gipsmult(list(S3, S4), c(100, 100), perm = perm1) g2 <- gipsmult(list(S3, S4), c(100, 100), perm = perm2) log_posteriori_of_gipsmult(g1) # -148.6485 log_posteriori_of_gipsmult(g2) # -145.3019, this is the MAP Estimator exp(log_posteriori_of_gipsmult(g2) - log_posteriori_of_gipsmult(g1)) # 28.406 # g2 is 28.406 times more likely than g1. # This is the expected outcome because S1[1,1] and S2[1,1] # are very close to S1[2,2] and S2[2,2] respectively.
gipsmult objectPlot heatmaps of the MAP covariance matrices estimator
or the convergence of the optimization method.
The plot depends on the type argument.
## S3 method for class 'gipsmult' plot( x, type = NA, logarithmic_y = TRUE, logarithmic_x = FALSE, color = NULL, title_text = "Convergence plot", xlabel = NULL, ylabel = NULL, show_legend = TRUE, ylim = NULL, xlim = NULL, ... )## S3 method for class 'gipsmult' plot( x, type = NA, logarithmic_y = TRUE, logarithmic_x = FALSE, color = NULL, title_text = "Convergence plot", xlabel = NULL, ylabel = NULL, show_legend = TRUE, ylim = NULL, xlim = NULL, ... )
x |
Object of a |
type |
A character vector of length 1. One of
The default value is |
logarithmic_y, logarithmic_x
|
A boolean. Sets the axis of the plots in logarithmic scale. |
color |
Vector of colors to be used to plot lines. |
title_text |
Text to be in the title of the plot. |
xlabel |
Text to be on the bottom of the plot. |
ylabel |
Text to be on the left of the plot. |
show_legend |
A boolean. Whether or not to show a legend. |
ylim |
Limits of the y axis. When |
xlim |
Limits of the x axis. When |
... |
Additional arguments passed to other various elements of the plot. |
When type is one of "best", "all", "both" or "n0",
returns an invisible NULL.
When type is one of "heatmap", "MLE" or "block_heatmap",
returns an object of class ggplot.
find_MAP() - Usually, the plot.gipsmult()
is called on the output of find_MAP().
gipsmult() - The constructor of a gipsmult class.
The gipsmult object is used as the x parameter.
require("MASS") # for mvrnorm() perm_size <- 6 mu1 <- runif(6, -10, 10) mu2 <- runif(6, -10, 10) # Assume we don't know the means sigma1 <- 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 ) sigma2 <- matrix( data = c( 1.0, 0.5, 0.2, 0.0, 0.2, 0.5, 0.5, 1.0, 0.5, 0.2, 0.0, 0.2, 0.2, 0.5, 1.0, 0.5, 0.2, 0.0, 0.0, 0.2, 0.5, 1.0, 0.5, 0.2, 0.2, 0.0, 0.2, 0.5, 1.0, 0.5, 0.5, 0.2, 0.0, 0.2, 0.5, 1.0 ), nrow = perm_size, byrow = TRUE ) # sigma1 and sigma2 are matrices invariant under permutation (1,2,3,4,5,6) numbers_of_observations <- c(21, 37) Z1 <- MASS::mvrnorm(numbers_of_observations[1], mu = mu1, Sigma = sigma1) Z2 <- MASS::mvrnorm(numbers_of_observations[2], mu = mu2, Sigma = sigma2) S1 <- cov(Z1) S2 <- cov(Z2) # Assume we have to estimate the mean g <- gipsmult(list(S1, S2), numbers_of_observations) if (require("graphics")) { plot(g, type = "MLE") } g_map <- find_MAP(g, max_iter = 30, show_progress_bar = FALSE, optimizer = "hill_climbing") if (require("graphics")) { plot(g_map, type = "both", logarithmic_x = TRUE) } if (require("graphics")) { 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") if (require("graphics")) { plot(g_map_MH, type = "n0") }require("MASS") # for mvrnorm() perm_size <- 6 mu1 <- runif(6, -10, 10) mu2 <- runif(6, -10, 10) # Assume we don't know the means sigma1 <- 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 ) sigma2 <- matrix( data = c( 1.0, 0.5, 0.2, 0.0, 0.2, 0.5, 0.5, 1.0, 0.5, 0.2, 0.0, 0.2, 0.2, 0.5, 1.0, 0.5, 0.2, 0.0, 0.0, 0.2, 0.5, 1.0, 0.5, 0.2, 0.2, 0.0, 0.2, 0.5, 1.0, 0.5, 0.5, 0.2, 0.0, 0.2, 0.5, 1.0 ), nrow = perm_size, byrow = TRUE ) # sigma1 and sigma2 are matrices invariant under permutation (1,2,3,4,5,6) numbers_of_observations <- c(21, 37) Z1 <- MASS::mvrnorm(numbers_of_observations[1], mu = mu1, Sigma = sigma1) Z2 <- MASS::mvrnorm(numbers_of_observations[2], mu = mu2, Sigma = sigma2) S1 <- cov(Z1) S2 <- cov(Z2) # Assume we have to estimate the mean g <- gipsmult(list(S1, S2), numbers_of_observations) if (require("graphics")) { plot(g, type = "MLE") } g_map <- find_MAP(g, max_iter = 30, show_progress_bar = FALSE, optimizer = "hill_climbing") if (require("graphics")) { plot(g_map, type = "both", logarithmic_x = TRUE) } if (require("graphics")) { 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") if (require("graphics")) { plot(g_map_MH, type = "n0") }
gipsmult objectPrinting function for a gipsmult class.
## S3 method for class 'gipsmult' print( x, digits = 3, compare_to_original = TRUE, log_value = FALSE, oneline = FALSE, ... )## S3 method for class 'gipsmult' print( x, digits = 3, compare_to_original = TRUE, log_value = FALSE, oneline = FALSE, ... )
x |
An object of a |
digits |
The number of digits after the comma
for a posteriori to be presented. It can be negative.
By default, |
compare_to_original |
A logical. Whether to print how many times more likely is the current permutation compared to:
|
log_value |
A logical. Whether to print the logarithmic value.
Default to |
oneline |
A logical. Whether to print in
one or multiple lines. Default to |
... |
The additional arguments passed to |
Returns an invisible NULL.
find_MAP() - The function that makes
an optimized gipsmult object out of the unoptimized one.
Ss <- list( matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE), matrix(c(2, 1, 3, 7), nrow = 2, byrow = TRUE) ) noo <- c(10, 13) g <- gipsmult(Ss, noo, perm = "(12)") print(g, digits = 4, oneline = TRUE)Ss <- list( matrix(c(1, 0.5, 0.5, 2), nrow = 2, byrow = TRUE), matrix(c(2, 1, 3, 7), nrow = 2, byrow = TRUE) ) noo <- c(10, 13) g <- gipsmult(Ss, noo, perm = "(12)") print(g, digits = 4, oneline = TRUE)