Stan Math Library  2.20.0
reverse mode automatic differentiation
gp_matern52_cov.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_GP_MATERN52_COV_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_GP_MATERN52_COV_HPP
3 
15 #include <cmath>
16 #include <vector>
17 
18 namespace stan {
19 namespace math {
20 
41 template <typename T_x, typename T_s, typename T_l>
42 inline typename Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type,
43  Eigen::Dynamic, Eigen::Dynamic>
44 gp_matern52_cov(const std::vector<T_x> &x, const T_s &sigma,
45  const T_l &length_scale) {
46  using std::exp;
47 
48  size_t x_size = x.size();
49  Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type, Eigen::Dynamic,
50  Eigen::Dynamic>
51  cov(x_size, x_size);
52 
53  if (x_size == 0)
54  return cov;
55 
56  for (size_t n = 0; n < x_size; ++n)
57  check_not_nan("gp_matern52_cov", "x", x[n]);
58 
59  check_positive_finite("gp_matern52_cov", "magnitude", sigma);
60  check_positive_finite("gp_matern52_cov", "length scale", length_scale);
61 
62  T_s sigma_sq = square(sigma);
63  T_l root_5_inv_l = sqrt(5.0) / length_scale;
64  T_l inv_l_sq_5_3 = 5.0 / (3.0 * square(length_scale));
65  T_l neg_root_5_inv_l = -sqrt(5.0) / length_scale;
66 
67  for (size_t i = 0; i < x_size; ++i) {
68  cov(i, i) = sigma_sq;
69  for (size_t j = i + 1; j < x_size; ++j) {
70  typename return_type<T_x>::type sq_distance
71  = squared_distance(x[i], x[j]);
72  typename return_type<T_x>::type dist = sqrt(sq_distance);
73  cov(i, j) = sigma_sq
74  * (1.0 + root_5_inv_l * dist + inv_l_sq_5_3 * sq_distance)
75  * exp(neg_root_5_inv_l * dist);
76  cov(j, i) = cov(i, j);
77  }
78  }
79  return cov;
80 }
81 
105 template <typename T_x, typename T_s, typename T_l>
106 inline typename Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type,
107  Eigen::Dynamic, Eigen::Dynamic>
108 gp_matern52_cov(const std::vector<Eigen::Matrix<T_x, Eigen::Dynamic, 1>> &x,
109  const T_s &sigma, const std::vector<T_l> &length_scale) {
110  using std::exp;
111 
112  size_t x_size = x.size();
113  Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type, Eigen::Dynamic,
114  Eigen::Dynamic>
115  cov(x_size, x_size);
116  if (x_size == 0)
117  return cov;
118 
119  size_t l_size = length_scale.size();
120  for (size_t n = 0; n < x_size; ++n)
121  check_not_nan("gp_matern52_cov", "x", x[n]);
122 
123  check_positive_finite("gp_matern52_cov", "magnitude", sigma);
124  check_positive_finite("gp_matern52_cov", "length scale", length_scale);
125  check_size_match("gp_matern52_cov", "x dimension", x[0].size(),
126  "number of length scales", l_size);
127 
128  T_s sigma_sq = square(sigma);
129  double root_5 = sqrt(5.0);
130  double five_thirds = 5.0 / 3.0;
131  double neg_root_5 = -root_5;
132 
133  std::vector<Eigen::Matrix<typename return_type<T_x, T_l>::type, -1, 1>> x_new
134  = divide_columns(x, length_scale);
135 
136  for (size_t i = 0; i < x_size; ++i) {
137  cov(i, i) = sigma_sq;
138  for (size_t j = i + 1; j < x_size; ++j) {
139  typename return_type<T_x, T_l>::type sq_distance
140  = squared_distance(x_new[i], x_new[j]);
141  typename return_type<T_x, T_l>::type dist = sqrt(sq_distance);
142  cov(i, j) = sigma_sq * (1.0 + root_5 * dist + five_thirds * sq_distance)
143  * exp(neg_root_5 * dist);
144  cov(j, i) = cov(i, j);
145  }
146  }
147  return cov;
148 }
149 
172 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
173 inline typename Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
174  Eigen::Dynamic, Eigen::Dynamic>
175 gp_matern52_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
176  const T_s &sigma, const T_l &length_scale) {
177  using std::exp;
178 
179  size_t x1_size = x1.size();
180  size_t x2_size = x2.size();
181  Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
182  Eigen::Dynamic, Eigen::Dynamic>
183  cov(x1_size, x2_size);
184 
185  if (x1_size == 0 || x2_size == 0)
186  return cov;
187 
188  for (size_t n = 0; n < x1_size; ++n)
189  check_not_nan("gp_matern52_cov", "x1", x1[n]);
190  for (size_t n = 0; n < x2_size; ++n)
191  check_not_nan("gp_matern52_cov", "x1", x2[n]);
192 
193  check_positive_finite("gp_matern52_cov", "magnitude", sigma);
194  check_positive_finite("gp_matern52_cov", "length scale", length_scale);
195 
196  T_s sigma_sq = square(sigma);
197  T_l root_5_inv_l = sqrt(5.0) / length_scale;
198  T_l inv_l_sq_5_3 = 5.0 / (3.0 * square(length_scale));
199  T_l neg_root_5_inv_l = -sqrt(5.0) / length_scale;
200 
201  for (size_t i = 0; i < x1_size; ++i) {
202  for (size_t j = 0; j < x2_size; ++j) {
203  typename return_type<T_x1, T_x2>::type sq_distance
204  = squared_distance(x1[i], x2[j]);
205  typename return_type<T_x1, T_x2>::type dist = sqrt(sq_distance);
206  cov(i, j) = sigma_sq
207  * (1.0 + root_5_inv_l * dist + inv_l_sq_5_3 * sq_distance)
208  * exp(neg_root_5_inv_l * dist);
209  }
210  }
211  return cov;
212 }
213 
240 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
241 inline typename Eigen::Matrix<
242  typename stan::return_type<T_x1, T_x2, T_s, T_l>::type, Eigen::Dynamic,
243  Eigen::Dynamic>
244 gp_matern52_cov(const std::vector<Eigen::Matrix<T_x1, Eigen::Dynamic, 1>> &x1,
245  const std::vector<Eigen::Matrix<T_x2, Eigen::Dynamic, 1>> &x2,
246  const T_s &sigma, const std::vector<T_l> &length_scale) {
247  using std::exp;
248 
249  size_t x1_size = x1.size();
250  size_t x2_size = x2.size();
251  Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
252  Eigen::Dynamic, Eigen::Dynamic>
253  cov(x1_size, x2_size);
254 
255  if (x1_size == 0 || x2_size == 0)
256  return cov;
257 
258  size_t l_size = length_scale.size();
259 
260  for (size_t n = 0; n < x1_size; ++n)
261  check_not_nan("gp_matern52_cov", "x1", x1[n]);
262  for (size_t n = 0; n < x2_size; ++n)
263  check_not_nan("gp_matern52_cov", "x1", x2[n]);
264 
265  check_positive_finite("gp_matern52_cov", "magnitude", sigma);
266  check_positive_finite("gp_matern52_cov", "length scale", length_scale);
267  check_size_match("gp_matern52_cov", "x dimension", x1[0].size(),
268  "number of length scales", l_size);
269  check_size_match("gp_matern52_cov", "x dimension", x2[0].size(),
270  "number of length scales", l_size);
271 
272  T_s sigma_sq = square(sigma);
273  double root_5 = sqrt(5.0);
274  double five_thirds = 5.0 / 3.0;
275  double neg_root_5 = -root_5;
276 
277  std::vector<Eigen::Matrix<typename return_type<T_x1, T_l>::type, -1, 1>>
278  x1_new = divide_columns(x1, length_scale);
279  std::vector<Eigen::Matrix<typename return_type<T_x2, T_l>::type, -1, 1>>
280  x2_new = divide_columns(x2, length_scale);
281 
282  for (size_t i = 0; i < x1_size; ++i) {
283  for (size_t j = 0; j < x2_size; ++j) {
284  typename return_type<T_x1, T_x2, T_l>::type sq_distance
285  = squared_distance(x1_new[i], x2_new[j]);
286  typename return_type<T_x1, T_x2, T_l>::type dist = sqrt(sq_distance);
287  cov(i, j) = sigma_sq * (1.0 + root_5 * dist + five_thirds * sq_distance)
288  * exp(neg_root_5 * dist);
289  }
290  }
291  return cov;
292 }
293 } // namespace math
294 } // namespace stan
295 #endif
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:13
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
std::vector< Eigen::Matrix< typename return_type< T_x, T_v, double >::type, Eigen::Dynamic, 1 > > divide_columns(const std::vector< Eigen::Matrix< T_x, Eigen::Dynamic, 1 >> &x, const std::vector< T_v > &vec)
Takes Stan data type vector[n] x[D] and divides column vector in x element-wise by the values in vec...
fvar< T > squared_distance(const Eigen::Matrix< fvar< T >, R, C > &v1, const Eigen::Matrix< double, R, C > &v2)
Returns the squared distance between the specified vectors of the same dimensions.
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.
Eigen::Matrix< typename return_type< T_x, T_s, T_l >::type, Eigen::Dynamic, Eigen::Dynamic > gp_matern52_cov(const std::vector< T_x > &x, const T_s &sigma, const T_l &length_scale)
Returns a Matern 5/2 covariance matrix with one input vector.
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17

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