Stan Math Library  2.20.0
reverse mode automatic differentiation
gp_exponential_cov.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_GP_EXPONENTIAL_COV_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_GP_EXPONENTIAL_COV_HPP
3 
11 #include <cmath>
12 #include <vector>
13 
14 namespace stan {
15 namespace math {
16 
34 template <typename T_x, typename T_s, typename T_l>
35 inline typename Eigen::Matrix<typename stan::return_type<T_x, T_s, T_l>::type,
36  Eigen::Dynamic, Eigen::Dynamic>
37 gp_exponential_cov(const std::vector<T_x> &x, const T_s &sigma,
38  const T_l &length_scale) {
39  using std::exp;
40  using std::pow;
41 
42  size_t x_size = size_of(x);
43  Eigen::Matrix<typename stan::return_type<T_x, T_s, T_l>::type, Eigen::Dynamic,
44  Eigen::Dynamic>
45  cov(x_size, x_size);
46  if (x_size == 0)
47  return cov;
48 
49  const char *function = "gp_exponential_cov";
50  size_t x_obs_size = size_of(x[0]);
51  for (size_t i = 0; i < x_size; ++i)
52  check_size_match(function, "x row", x_obs_size, "x's other row",
53  size_of(x[i]));
54 
55  for (size_t i = 0; i < x_size; ++i)
56  check_not_nan(function, "x", x[i]);
57 
58  check_positive_finite(function, "magnitude", sigma);
59  check_positive_finite(function, "length scale", length_scale);
60 
61  T_s sigma_sq = square(sigma);
62  T_l neg_inv_l = -1.0 / length_scale;
63 
64  for (size_t i = 0; i < x_size; ++i) {
65  cov(i, i) = sigma_sq;
66  for (size_t j = i + 1; j < x_size; ++j) {
67  cov(i, j) = sigma_sq * exp(neg_inv_l * distance(x[i], x[j]));
68  cov(j, i) = cov(i, j);
69  }
70  }
71  return cov;
72 }
73 
90 template <typename T_x, typename T_s, typename T_l>
91 inline typename Eigen::Matrix<typename stan::return_type<T_x, T_s, T_l>::type,
92  Eigen::Dynamic, Eigen::Dynamic>
93 gp_exponential_cov(const std::vector<Eigen::Matrix<T_x, -1, 1>> &x,
94  const T_s &sigma, const std::vector<T_l> &length_scale) {
95  using std::exp;
96  using std::pow;
97 
98  size_t x_size = size_of(x);
99  Eigen::Matrix<typename stan::return_type<T_x, T_s, T_l>::type, Eigen::Dynamic,
100  Eigen::Dynamic>
101  cov(x_size, x_size);
102  if (x_size == 0)
103  return cov;
104 
105  const char *function = "gp_exponential_cov";
106  for (size_t n = 0; n < x_size; ++n)
107  check_not_nan(function, "x", x[n]);
108 
109  check_positive_finite(function, "magnitude", sigma);
110  check_positive_finite(function, "length scale", length_scale);
111 
112  size_t l_size = length_scale.size();
113  check_size_match(function, "x dimension", x[0].size(),
114  "number of length scales", l_size);
115 
116  std::vector<Eigen::Matrix<typename return_type<T_x, T_l>::type, -1, 1>> x_new
117  = divide_columns(x, length_scale);
118 
119  T_s sigma_sq = square(sigma);
120  for (size_t i = 0; i < x_size; ++i) {
121  cov(i, i) = sigma_sq;
122  for (size_t j = i + 1; j < x_size; ++j) {
123  typename return_type<T_x, T_l>::type dist = distance(x_new[i], x_new[j]);
124  cov(i, j) = sigma_sq * exp(-dist);
125  cov(j, i) = cov(i, j);
126  }
127  }
128  return cov;
129 }
130 
151 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
152 inline typename Eigen::Matrix<
153  typename stan::return_type<T_x1, T_x2, T_s, T_l>::type, Eigen::Dynamic,
154  Eigen::Dynamic>
155 gp_exponential_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
156  const T_s &sigma, const T_l &length_scale) {
157  using std::exp;
158  using std::pow;
159 
160  size_t x1_size = size_of(x1);
161  size_t x2_size = size_of(x2);
162  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_s, T_l>::type,
163  Eigen::Dynamic, Eigen::Dynamic>
164  cov(x1_size, x2_size);
165  if (x1_size == 0 || x2_size == 0)
166  return cov;
167 
168  const char *function = "gp_exponential_cov";
169  size_t x1_obs_size = size_of(x1[0]);
170  for (size_t i = 0; i < x1_size; ++i)
171  check_size_match(function, "x1's row", x1_obs_size, "x1's other row",
172  size_of(x1[i]));
173  for (size_t i = 0; i < x2_size; ++i)
174  check_size_match(function, "x1's row", x1_obs_size, "x2's other row",
175  size_of(x2[i]));
176 
177  for (size_t n = 0; n < x1_size; ++n)
178  check_not_nan(function, "x1", x1[n]);
179  for (size_t n = 0; n < x2_size; ++n)
180  check_not_nan(function, "x2", x2[n]);
181 
182  check_positive_finite(function, "magnitude", sigma);
183  check_positive_finite(function, "length scale", length_scale);
184 
185  T_s sigma_sq = square(sigma);
186  T_l neg_inv_l = -1.0 / length_scale;
187  for (size_t i = 0; i < x1_size; ++i) {
188  for (size_t j = 0; j < x2_size; ++j) {
189  cov(i, j) = sigma_sq * exp(neg_inv_l * distance(x1[i], x2[j]));
190  }
191  }
192  return cov;
193 }
194 
216 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
217 inline typename Eigen::Matrix<typename return_type<T_x1, T_x2, T_s, T_l>::type,
218  Eigen::Dynamic, Eigen::Dynamic>
219 gp_exponential_cov(const std::vector<Eigen::Matrix<T_x1, -1, 1>> &x1,
220  const std::vector<Eigen::Matrix<T_x2, -1, 1>> &x2,
221  const T_s &sigma, const std::vector<T_l> &length_scale) {
222  using std::exp;
223  using std::pow;
224 
225  size_t x1_size = size_of(x1);
226  size_t x2_size = size_of(x2);
227  size_t l_size = size_of(length_scale);
228  Eigen::Matrix<typename stan::return_type<T_x1, T_x2, T_s, T_l>::type,
229  Eigen::Dynamic, Eigen::Dynamic>
230  cov(x1_size, x2_size);
231  if (x1_size == 0 || x2_size == 0)
232  return cov;
233 
234  const char *function = "gp_exponential_cov";
235  for (size_t n = 0; n < x1_size; ++n)
236  check_not_nan(function, "x1", x1[n]);
237  for (size_t n = 0; n < x2_size; ++n)
238  check_not_nan(function, "x2", x2[n]);
239 
240  check_positive_finite(function, "magnitude", sigma);
241  check_positive_finite(function, "length scale", length_scale);
242 
243  for (size_t i = 0; i < x1_size; ++i)
244  check_size_match(function, "x1's row", size_of(x1[i]),
245  "number of length scales", l_size);
246  for (size_t i = 0; i < x2_size; ++i)
247  check_size_match(function, "x2's row", size_of(x2[i]),
248  "number of length scales", l_size);
249 
250  T_s sigma_sq = square(sigma);
251  T_l temp;
252 
253  std::vector<Eigen::Matrix<typename return_type<T_x1, T_l>::type, -1, 1>>
254  x1_new = divide_columns(x1, length_scale);
255  std::vector<Eigen::Matrix<typename return_type<T_x2, T_l>::type, -1, 1>>
256  x2_new = divide_columns(x2, length_scale);
257 
258  for (size_t i = 0; i < x1_size; ++i) {
259  for (size_t j = 0; j < x2_size; ++j) {
260  cov(i, j) = sigma_sq * exp(-distance(x1_new[i], x2_new[j]));
261  }
262  }
263  return cov;
264 }
265 } // namespace math
266 } // namespace stan
267 #endif
Eigen::Matrix< typename stan::return_type< T_x, T_s, T_l >::type, Eigen::Dynamic, Eigen::Dynamic > gp_exponential_cov(const std::vector< T_x > &x, const T_s &sigma, const T_l &length_scale)
Returns a Matern exponential covariance Matrix.
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
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
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.