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

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

g

An object of a gips class.

Value

The logarithm of an unnormalized A Posteriori.

Details

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

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

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
#> [1] -134.1615
log_posteriori_of_gips(g2) # -138.1695
#> [1] -138.1695

exp(log_posteriori_of_gips(g1) - log_posteriori_of_gips(g2)) # 55.0
#> [1] 55.03601
# 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 permutation () is 55.036 times more likely than the (1,2) permutation.
# 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
#> [1] -98.40984
log_posteriori_of_gips(g2) # -95.92039, this is the MAP Estimator
#> [1] -95.92039

exp(log_posteriori_of_gips(g2) - log_posteriori_of_gips(g1)) # 12.05
#> [1] 12.0546
# 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 permutation (1,2) is 12.055 times more likely than the () permutation.
# The same result, but presented in a more pleasant way