1 #ifndef STAN_MATH_PRIM_FUN_CORR_MATRIX_FREE_HPP
2 #define STAN_MATH_PRIM_FUN_CORR_MATRIX_FREE_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/Eigen.hpp>
7 #include <stan/math/prim/fun/factor_cov_matrix.hpp>
8 #include <cmath>
9
10 namespace stan {
11 namespace math {
12
13 /**
14 * Return the vector of unconstrained partial correlations that
15 * define the specified correlation matrix when transformed.
16 *
17 * <p>The constraining transform is defined as for
18 * <code>corr_matrix_constrain(Matrix, size_t)</code>. The
19 * inverse transform in this function is simpler in that it only
20 * needs to compute the \f$k \choose 2\f$ partial correlations
21 * and then free those.
22 *
23 * @tparam T type of the matrix (must be derived from \c Eigen::MatrixBase)
24 * @param y The correlation matrix to free.
25 * @return Vector of unconstrained values that produce the
26 * specified correlation matrix when transformed.
27 * @throw std::domain_error if the correlation matrix has no
28 * elements or is not a square matrix.
29 * @throw std::runtime_error if the correlation matrix cannot be
30 * factorized by factor_cov_matrix() or if the sds returned by
31 * factor_cov_matrix() on log scale are unconstrained.
32 */
33 template <typename T, require_eigen_t<T>* = nullptr>
corr_matrix_free(const T & y)34 Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, 1> corr_matrix_free(const T& y) {
35 using Eigen::Array;
36 using Eigen::Dynamic;
37
38 check_square("corr_matrix_free", "y", y);
39 check_nonzero_size("corr_matrix_free", "y", y);
40
41 Eigen::Index k = y.rows();
42 Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
43 Eigen::Matrix<value_type_t<T>, Dynamic, 1> x(k_choose_2);
44 Array<value_type_t<T>, Dynamic, 1> sds(k);
45 bool successful = factor_cov_matrix(y, x.array(), sds);
46 if (!successful) {
47 throw_domain_error("corr_matrix_free", "factor_cov_matrix failed on y", y,
48 "");
49 }
50 check_bounded("corr_matrix_free", "log(sd)", sds, -CONSTRAINT_TOLERANCE,
51 CONSTRAINT_TOLERANCE);
52 return x;
53 }
54
55 /**
56 * Overload of `corr_matrix_free()` to untransform each matrix
57 * in a standard vector.
58 * @tparam T A standard vector with with a `value_type` which inherits from
59 * `Eigen::MatrixBase`.
60 * @param x The standard vector to untransform.
61 */
62 template <typename T, require_std_vector_t<T>* = nullptr>
corr_matrix_free(const T & x)63 auto corr_matrix_free(const T& x) {
64 return apply_vector_unary<T>::apply(
65 x, [](auto&& v) { return corr_matrix_free(v); });
66 }
67
68 } // namespace math
69 } // namespace stan
70
71 #endif
72