Stan Math Library  2.20.0
reverse mode automatic differentiation
gp_periodic_cov.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_GP_PERIODIC_COV_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_GP_PERIODIC_COV_HPP
3 
4 #include <math.h>
16 #include <cmath>
17 
18 #include <vector>
19 
20 namespace stan {
21 namespace math {
22 
48 template <typename T_x, typename T_sigma, typename T_l, typename T_p>
49 inline typename Eigen::Matrix<
50  typename stan::return_type<T_x, T_sigma, T_l, T_p>::type, Eigen::Dynamic,
51  Eigen::Dynamic>
52 gp_periodic_cov(const std::vector<T_x> &x, const T_sigma &sigma, const T_l &l,
53  const T_p &p) {
54  using std::exp;
55  const char *fun = "gp_periodic_cov";
56  check_positive(fun, "signal standard deviation", sigma);
57  check_positive(fun, "length-scale", l);
58  check_positive(fun, "period", p);
59  for (size_t n = 0; n < x.size(); ++n)
60  check_not_nan(fun, "element of x", x[n]);
61 
62  Eigen::Matrix<typename stan::return_type<T_x, T_sigma, T_l, T_p>::type,
63  Eigen::Dynamic, Eigen::Dynamic>
64  cov(x.size(), x.size());
65 
66  size_t x_size = x.size();
67  if (x_size == 0)
68  return cov;
69 
70  T_sigma sigma_sq = square(sigma);
71  T_l neg_two_inv_l_sq = -2.0 * inv_square(l);
72  T_p pi_div_p = pi() / p;
73 
74  for (size_t j = 0; j < x_size; ++j) {
75  cov(j, j) = sigma_sq;
76  for (size_t i = j + 1; i < x_size; ++i) {
77  cov(i, j) = sigma_sq
78  * exp(square(sin(pi_div_p * distance(x[i], x[j])))
79  * neg_two_inv_l_sq);
80  cov(j, i) = cov(i, j);
81  }
82  }
83  return cov;
84 }
85 
117 template <typename T_x1, typename T_x2, typename T_sigma, typename T_l,
118  typename T_p>
119 inline typename Eigen::Matrix<
121  Eigen::Dynamic, Eigen::Dynamic>
122 gp_periodic_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
123  const T_sigma &sigma, const T_l &l, const T_p &p) {
124  using std::exp;
125  const char *fun = "gp_periodic_cov";
126  check_positive(fun, "signal standard deviation", sigma);
127  check_positive(fun, "length-scale", l);
128  check_positive(fun, "period", p);
129  for (size_t n = 0; n < x1.size(); ++n)
130  check_not_nan(fun, "element of x1", x1[n]);
131  for (size_t n = 0; n < x2.size(); ++n)
132  check_not_nan(fun, "element of x2", x2[n]);
133 
134  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_sigma, T_l, T_p>::type,
135  Eigen::Dynamic, Eigen::Dynamic>
136  cov(x1.size(), x2.size());
137  if (x1.size() == 0 || x2.size() == 0)
138  return cov;
139 
140  T_sigma sigma_sq = square(sigma);
141  T_l neg_two_inv_l_sq = -2.0 * inv_square(l);
142  T_p pi_div_p = pi() / p;
143 
144  for (size_t i = 0; i < x1.size(); ++i) {
145  for (size_t j = 0; j < x2.size(); ++j) {
146  cov(i, j) = sigma_sq
147  * exp(square(sin(pi_div_p * distance(x1[i], x2[j])))
148  * neg_two_inv_l_sq);
149  }
150  }
151  return cov;
152 }
153 } // namespace math
154 } // namespace stan
155 
156 #endif
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:11
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
Definition: return_type.hpp:36
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
fvar< T > inv_square(const fvar< T > &x)
Definition: inv_square.hpp:12
boost::math::tools::promote_args< T1, T2 >::type distance(const Eigen::Matrix< T1, R1, C1 > &v1, const Eigen::Matrix< T2, R2, C2 > &v2)
Returns the distance between the specified vectors.
Definition: distance.hpp:23
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
double pi()
Return the value of pi.
Definition: constants.hpp:80
Eigen::Matrix< typename stan::return_type< T_x, T_sigma, T_l, T_p >::type, 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 .

     [ Stan Home Page ] © 2011–2018, Stan Development Team.