Automatic Differentiation
 
Loading...
Searching...
No Matches
factor_U.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_FACTOR_U_HPP
2#define STAN_MATH_PRIM_FUN_FACTOR_U_HPP
3
7#include <cstddef>
8#include <stdexcept>
9#include <vector>
10
11namespace stan {
12namespace math {
13
22template <typename T_U, typename T_CPCs, require_eigen_t<T_U>* = nullptr,
23 require_eigen_vector_t<T_CPCs>* = nullptr,
24 require_vt_same<T_U, T_CPCs>* = nullptr>
25void factor_U(const T_U& U, T_CPCs&& CPCs) {
26 size_t K = U.rows();
27 size_t position = 0;
28 size_t pull = K - 1;
29
30 const Eigen::Ref<const plain_type_t<T_U>>& U_ref = U;
31
32 if (K == 2) {
33 CPCs(0) = atanh(U_ref(0, 1));
34 return;
35 }
36
37 Eigen::Array<value_type_t<T_U>, 1, Eigen::Dynamic> temp
38 = U_ref.row(0).tail(pull);
39
40 CPCs.head(pull) = temp;
41
42 Eigen::Array<value_type_t<T_U>, Eigen::Dynamic, 1> acc(K);
43 acc(0) = -0.0;
44 acc.tail(pull) = 1.0 - temp.square();
45 for (size_t i = 1; i < (K - 1); i++) {
46 position += pull;
47 pull--;
48 temp = U_ref.row(i).tail(pull);
49 temp /= sqrt(acc.tail(pull) / acc(i));
50 CPCs.segment(position, pull) = temp;
51 acc.tail(pull) *= 1.0 - temp.square();
52 }
53 CPCs = 0.5 * ((1.0 + CPCs) / (1.0 - CPCs)).log(); // now unbounded
54}
55
56} // namespace math
57} // namespace stan
58
59#endif
fvar< T > atanh(const fvar< T > &x)
Return inverse hyperbolic tangent of specified value.
Definition atanh.hpp:24
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
fvar< T > sqrt(const fvar< T > &x)
Definition sqrt.hpp:17
void factor_U(const T_U &U, T_CPCs &&CPCs)
This function is intended to make starting values, given a unit upper-triangular matrix U such that U...
Definition factor_U.hpp:25
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...