Stan Math Library  2.20.0
reverse mode automatic differentiation
gp_matern32_cov.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_GP_MATERN32_COV_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_GP_MATERN32_COV_HPP
3 
15 #include <cmath>
16 #include <vector>
17 
18 namespace stan {
19 namespace math {
20 
39 template <typename T_x, typename T_s, typename T_l>
40 inline typename Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type,
41  Eigen::Dynamic, Eigen::Dynamic>
42 gp_matern32_cov(const std::vector<T_x> &x, const T_s &sigma,
43  const T_l &length_scale) {
44  using std::exp;
45  using std::pow;
46 
47  size_t x_size = size_of(x);
48  Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type, Eigen::Dynamic,
49  Eigen::Dynamic>
50  cov(x_size, x_size);
51 
52  if (x_size == 0)
53  return cov;
54 
55  const char *function = "gp_matern32_cov";
56  size_t x_obs_size = size_of(x[0]);
57  for (size_t i = 0; i < x_size; ++i)
58  check_size_match(function, "x row", x_obs_size, "x's other row",
59  size_of(x[i]));
60 
61  for (size_t n = 0; n < x_size; ++n)
62  check_not_nan(function, "x", x[n]);
63 
64  check_positive_finite(function, "magnitude", sigma);
65  check_positive_finite(function, "length scale", length_scale);
66 
67  T_s sigma_sq = square(sigma);
68  T_l root_3_inv_l = sqrt(3.0) / length_scale;
69  T_l neg_root_3_inv_l = -1.0 * sqrt(3.0) / length_scale;
70 
71  for (size_t i = 0; i < x_size; ++i) {
72  cov(i, i) = sigma_sq;
73  for (size_t j = i + 1; j < x_size; ++j) {
74  typename return_type<T_x>::type dist = distance(x[i], x[j]);
75  cov(i, j) = sigma_sq * (1.0 + root_3_inv_l * dist)
76  * exp(neg_root_3_inv_l * dist);
77  cov(j, i) = cov(i, j);
78  }
79  }
80  return cov;
81 }
82 
101 template <typename T_x, typename T_s, typename T_l>
102 inline typename Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type,
103  Eigen::Dynamic, Eigen::Dynamic>
104 gp_matern32_cov(const std::vector<Eigen::Matrix<T_x, -1, 1>> &x,
105  const T_s &sigma, const std::vector<T_l> &length_scale) {
106  using std::exp;
107 
108  size_t x_size = size_of(x);
109  Eigen::Matrix<typename return_type<T_x, T_s, T_l>::type, Eigen::Dynamic,
110  Eigen::Dynamic>
111  cov(x_size, x_size);
112 
113  if (x_size == 0)
114  return cov;
115  const char *function = "gp_matern32_cov";
116  size_t l_size = length_scale.size();
117  for (size_t n = 0; n < x_size; ++n)
118  check_not_nan(function, "x", x[n]);
119 
120  check_positive_finite(function, "magnitude", sigma);
121  check_positive_finite(function, "length scale", length_scale);
122 
123  for (size_t n = 0; n < x_size; ++n)
124  check_not_nan(function, "length scale", length_scale[n]);
125 
126  check_size_match(function, "x dimension", size_of(x[0]),
127  "number of length scales", l_size);
128 
129  T_s sigma_sq = square(sigma);
130  double root_3 = sqrt(3.0);
131  double neg_root_3 = -1.0 * sqrt(3.0);
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  for (size_t j = i; j < x_size; ++j) {
138  typename return_type<T_x, T_l>::type dist = distance(x_new[i], x_new[j]);
139  cov(i, j) = sigma_sq * (1.0 + root_3 * dist) * exp(neg_root_3 * dist);
140  cov(j, i) = cov(i, j);
141  }
142  }
143  return cov;
144 }
145 
170 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
171 inline typename Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
172  Eigen::Dynamic, Eigen::Dynamic>
173 gp_matern32_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
174  const T_s &sigma, const T_l &length_scale) {
175  using std::exp;
176 
177  size_t x1_size = size_of(x1);
178  size_t x2_size = size_of(x2);
179  Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
180  Eigen::Dynamic, Eigen::Dynamic>
181  cov(x1_size, x2_size);
182 
183  if (x1_size == 0 || x2_size == 0)
184  return cov;
185 
186  const char *function = "gp_matern32_cov";
187  size_t x1_obs_size = size_of(x1[0]);
188  for (size_t i = 0; i < x1_size; ++i)
189  check_size_match(function, "x1's row", x1_obs_size, "x1's other row",
190  size_of(x1[i]));
191  for (size_t i = 0; i < x2_size; ++i)
192  check_size_match(function, "x1's row", x1_obs_size, "x2's other row",
193  size_of(x2[i]));
194 
195  for (size_t n = 0; n < x1_size; ++n)
196  check_not_nan(function, "x1", x1[n]);
197  for (size_t n = 0; n < x2_size; ++n)
198  check_not_nan(function, "x2", x2[n]);
199 
200  check_positive_finite(function, "magnitude", sigma);
201  check_positive_finite(function, "length scale", length_scale);
202 
203  T_s sigma_sq = square(sigma);
204  T_l root_3_inv_l_sq = sqrt(3.0) / length_scale;
205  T_l neg_root_3_inv_l_sq = -1.0 * sqrt(3.0) / length_scale;
206 
207  for (size_t i = 0; i < x1_size; ++i) {
208  for (size_t j = 0; j < x2_size; ++j) {
209  typename return_type<T_x1, T_x2>::type dist = distance(x1[i], x2[j]);
210  cov(i, j) = sigma_sq * (1.0 + root_3_inv_l_sq * dist)
211  * exp(neg_root_3_inv_l_sq * dist);
212  }
213  }
214  return cov;
215 }
216 
241 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
242 inline typename Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
243  Eigen::Dynamic, Eigen::Dynamic>
244 gp_matern32_cov(const std::vector<Eigen::Matrix<T_x1, -1, 1>> &x1,
245  const std::vector<Eigen::Matrix<T_x2, -1, 1>> &x2,
246  const T_s &sigma, const std::vector<T_l> &length_scale) {
247  using std::exp;
248 
249  size_t x1_size = size_of(x1);
250  size_t x2_size = size_of(x2);
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  const char *function = "gp_matern_32_cov";
259  for (size_t n = 0; n < x1_size; ++n)
260  check_not_nan(function, "x1", x1[n]);
261  for (size_t n = 0; n < x2_size; ++n)
262  check_not_nan(function, "x2", x2[n]);
263 
264  check_positive_finite(function, "magnitude", sigma);
265  check_positive_finite(function, "length scale", length_scale);
266 
267  size_t l_size = length_scale.size();
268  for (size_t n = 0; n < l_size; ++n)
269  check_not_nan(function, "length scale", length_scale[n]);
270 
271  for (size_t i = 0; i < x1_size; ++i)
272  check_size_match(function, "x1's row", size_of(x1[i]),
273  "number of length scales", l_size);
274  for (size_t i = 0; i < x2_size; ++i)
275  check_size_match(function, "x2's row", size_of(x2[i]),
276  "number of length scales", l_size);
277 
278  T_s sigma_sq = square(sigma);
279  double root_3 = sqrt(3.0);
280  double neg_root_3 = -1.0 * sqrt(3.0);
281 
282  std::vector<Eigen::Matrix<typename return_type<T_x1, T_l>::type, -1, 1>>
283  x1_new = divide_columns(x1, length_scale);
284  std::vector<Eigen::Matrix<typename return_type<T_x2, T_l>::type, -1, 1>>
285  x2_new = divide_columns(x2, length_scale);
286 
287  for (size_t i = 0; i < x1_size; ++i) {
288  for (size_t j = 0; j < x2_size; ++j) {
290  = distance(x1_new[i], x2_new[j]);
291  cov(i, j) = sigma_sq * (1.0 + root_3 * dist) * exp(neg_root_3 * dist);
292  }
293  }
294  return cov;
295 }
296 } // namespace math
297 } // namespace stan
298 #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...
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.
size_t size_of(const T &x)
Returns the size of the provided vector or the constant 1 if the input argument is not a vector...
Definition: size_of.hpp:27
Eigen::Matrix< typename return_type< T_x, T_s, T_l >::type, Eigen::Dynamic, Eigen::Dynamic > gp_matern32_cov(const std::vector< T_x > &x, const T_s &sigma, const T_l &length_scale)
Returns a Matern 3/2 covariance matrix.
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
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16

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