Automatic Differentiation
 
Loading...
Searching...
No Matches
gp_periodic_cov.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_GP_PERIODIC_COV_HPP
2#define STAN_MATH_PRIM_FUN_GP_PERIODIC_COV_HPP
3
13#include <cmath>
14#include <vector>
15
16namespace stan {
17namespace math {
18
44template <typename T_x, typename T_sigma, typename T_l, typename T_p>
45inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l, T_p>,
46 Eigen::Dynamic, Eigen::Dynamic>
47gp_periodic_cov(const std::vector<T_x> &x, const T_sigma &sigma, const T_l &l,
48 const T_p &p) {
49 using std::exp;
50 using std::sin;
51 const char *fun = "gp_periodic_cov";
52 check_positive(fun, "signal standard deviation", sigma);
53 check_positive(fun, "length-scale", l);
54 check_positive(fun, "period", p);
55 for (size_t n = 0; n < x.size(); ++n) {
56 check_not_nan(fun, "element of x", x[n]);
57 }
58
59 Eigen::Matrix<return_type_t<T_x, T_sigma, T_l, T_p>, Eigen::Dynamic,
60 Eigen::Dynamic>
61 cov(x.size(), x.size());
62
63 size_t x_size = x.size();
64 if (x_size == 0) {
65 return cov;
66 }
67
68 T_sigma sigma_sq = square(sigma);
69 T_l neg_two_inv_l_sq = -2.0 * inv_square(l);
70 T_p pi_div_p = pi() / p;
71 size_t block_size = 10;
72
73 for (size_t jb = 0; jb < x_size; jb += block_size) {
74 for (size_t ib = jb; ib < x_size; ib += block_size) {
75 size_t j_end = std::min(x_size, jb + block_size);
76 for (size_t j = jb; j < j_end; ++j) {
77 cov.coeffRef(j, j) = sigma_sq;
78 size_t i_end = std::min(x_size, ib + block_size);
79 for (size_t i = std::max(ib, j + 1); i < i_end; ++i) {
80 cov.coeffRef(j, i) = cov.coeffRef(i, j)
81 = sigma_sq
82 * exp(square(sin(pi_div_p * distance(x[i], x[j])))
83 * neg_two_inv_l_sq);
84 }
85 }
86 }
87 }
88 return cov;
89}
90
123template <typename T_x1, typename T_x2, typename T_sigma, typename T_l,
124 typename T_p>
125inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l, T_p>,
126 Eigen::Dynamic, Eigen::Dynamic>
127gp_periodic_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
128 const T_sigma &sigma, const T_l &l, const T_p &p) {
129 using std::exp;
130 using std::sin;
131 const char *fun = "gp_periodic_cov";
132 check_positive(fun, "signal standard deviation", sigma);
133 check_positive(fun, "length-scale", l);
134 check_positive(fun, "period", p);
135 for (size_t n = 0; n < x1.size(); ++n) {
136 check_not_nan(fun, "element of x1", x1[n]);
137 }
138 for (size_t n = 0; n < x2.size(); ++n) {
139 check_not_nan(fun, "element of x2", x2[n]);
140 }
141
142 Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l, T_p>, Eigen::Dynamic,
143 Eigen::Dynamic>
144 cov(x1.size(), x2.size());
145 if (x1.size() == 0 || x2.size() == 0) {
146 return cov;
147 }
148
149 T_sigma sigma_sq = square(sigma);
150 T_l neg_two_inv_l_sq = -2.0 * inv_square(l);
151 T_p pi_div_p = pi() / p;
152 size_t block_size = 10;
153
154 for (size_t ib = 0; ib < x1.size(); ib += block_size) {
155 for (size_t jb = 0; jb < x2.size(); jb += block_size) {
156 size_t j_end = std::min(x2.size(), jb + block_size);
157 for (size_t j = jb; j < j_end; ++j) {
158 size_t i_end = std::min(x1.size(), ib + block_size);
159 for (size_t i = ib; i < i_end; ++i) {
160 cov.coeffRef(i, j)
161 = sigma_sq
162 * exp(square(sin(pi_div_p * distance(x1[i], x2[j])))
163 * neg_two_inv_l_sq);
164 }
165 }
166 }
167 }
168 return cov;
169}
170} // namespace math
171} // namespace stan
172
173#endif
fvar< T > sin(const fvar< T > &x)
Definition sin.hpp:14
fvar< T > inv_square(const fvar< T > &x)
auto distance(const T_a &a, const T_b &b)
Returns the distance between the specified vectors.
Definition distance.hpp:33
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
Eigen::Matrix< return_type_t< T_x, T_sigma, T_l, T_p >, Eigen::Dynamic, Eigen::Dynamic > gp_periodic_cov(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &l, const T_p &p)
Returns a periodic covariance matrix using the input .
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
static constexpr double pi()
Return the value of pi.
Definition constants.hpp:36
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...