read_corr_L.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_READ_CORR_L_HPP
#define STAN_MATH_PRIM_FUN_READ_CORR_L_HPP
 
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/prim/fun/log1m.hpp>
#include <stan/math/prim/fun/sqrt.hpp>
#include <stan/math/prim/fun/square.hpp>
#include <stan/math/prim/fun/sum.hpp>
#include <cstddef>
#include <cmath>
 
namespace stan {
namespace math {
 
template <typename T, require_eigen_vector_t<T>* = nullptr>
Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
    const T& CPCs,  // on (-1, 1)
    size_t K) {
  using T_scalar = value_type_t<T>;
  if (K == 0) {
    return {};
  }
  if (K == 1) {
    return Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic>::Identity(1,
                                                                             1);
  }
 
  using std::sqrt;
  Eigen::Array<T_scalar, Eigen::Dynamic, 1> temp;
  Eigen::Array<T_scalar, Eigen::Dynamic, 1> acc(K - 1);
  acc.setOnes();
  // Cholesky factor of correlation matrix
  Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic> L(K, K);
  L.setZero();
 
  size_t position = 0;
  size_t pull = K - 1;
 
  L(0, 0) = 1.0;
  L.col(0).tail(pull) = temp = CPCs.head(pull);
  acc.tail(pull) = T_scalar(1.0) - temp.square();
  for (size_t i = 1; i < (K - 1); i++) {
    position += pull;
    pull--;
    temp = CPCs.segment(position, pull);
    L(i, i) = sqrt(acc(i - 1));
    L.col(i).tail(pull) = temp * acc.tail(pull).sqrt();
    acc.tail(pull) *= T_scalar(1.0) - temp.square();
  }
  L(K - 1, K - 1) = sqrt(acc(K - 2));
  return L;
}
 
template <typename T, require_eigen_vector_t<T>* = nullptr>
Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
    const T& CPCs, size_t K, value_type_t<T>& log_prob) {
  using T_scalar = value_type_t<T>;
  if (K == 0) {
    return {};
  }
  if (K == 1) {
    return Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic>::Identity(1,
                                                                             1);
  }
 
  const Eigen::Ref<const plain_type_t<T>>& CPCs_ref = CPCs;
  size_t pos = 0;
  T_scalar acc = 0;
  // no need to abs() because this Jacobian determinant
  // is strictly positive (and triangular)
  // see inverse of Jacobian in equation 11 of LKJ paper
  for (size_t k = 1; k <= (K - 2); k++) {
    for (size_t i = k + 1; i <= K; i++) {
      acc += (K - k - 1) * log1m(square(CPCs_ref(pos)));
      pos++;
    }
  }
 
  log_prob += 0.5 * acc;
  return read_corr_L(CPCs_ref, K);
}
 
}  // namespace math
}  // namespace stan
 
#endif