1 #ifndef STAN_MATH_PRIM_FUN_GP_MATERN32_COV_HPP
2 #define STAN_MATH_PRIM_FUN_GP_MATERN32_COV_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/Eigen.hpp>
7 #include <stan/math/prim/fun/exp.hpp>
8 #include <stan/math/prim/fun/divide_columns.hpp>
9 #include <stan/math/prim/fun/distance.hpp>
10 #include <stan/math/prim/fun/size.hpp>
11 #include <stan/math/prim/fun/square.hpp>
12 #include <cmath>
13 #include <vector>
14
15 namespace stan {
16 namespace math {
17
18 /**
19 * Returns a Matern 3/2 covariance matrix
20 *
21 * \f[ k(x, x') = \sigma^2(1 +
22 * \frac{\sqrt{3}d(x, x')}{l})exp(-\frac{\sqrt{3}d(x, x')}{l})
23 * \f]
24 *
25 * where \f$ d(x, x') \f$ is the Euclidean distance.
26 *
27 * @tparam T_x type for each scalar
28 * @tparam T_s type of parameter of sigma
29 * @tparam T_l type of parameter length scale
30 *
31 * @param x std::vector of scalars that can be used in squared distance
32 * @param length_scale length scale
33 * @param sigma marginal standard deviation or magnitude
34 * @throw std::domain error if sigma <= 0, l <= 0, or x is nan or inf
35 */
36 template <typename T_x, typename T_s, typename T_l>
37 inline typename Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic,
38 Eigen::Dynamic>
gp_matern32_cov(const std::vector<T_x> & x,const T_s & sigma,const T_l & length_scale)39 gp_matern32_cov(const std::vector<T_x> &x, const T_s &sigma,
40 const T_l &length_scale) {
41 using std::exp;
42 using std::pow;
43
44 size_t x_size = stan::math::size(x);
45 Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic, Eigen::Dynamic>
46 cov(x_size, x_size);
47
48 if (x_size == 0) {
49 return cov;
50 }
51
52 const char *function = "gp_matern32_cov";
53 size_t x_obs_size = stan::math::size(x[0]);
54 for (size_t i = 0; i < x_size; ++i) {
55 check_size_match(function, "x row", x_obs_size, "x's other row",
56 stan::math::size(x[i]));
57 }
58
59 for (size_t n = 0; n < x_size; ++n) {
60 check_not_nan(function, "x", x[n]);
61 }
62
63 check_positive_finite(function, "magnitude", sigma);
64 check_positive_finite(function, "length scale", length_scale);
65
66 T_s sigma_sq = square(sigma);
67 T_l root_3_inv_l = std::sqrt(3.0) / length_scale;
68
69 size_t block_size = 10;
70 for (size_t jb = 0; jb < x_size; jb += block_size) {
71 for (size_t ib = jb; ib < x_size; ib += block_size) {
72 size_t j_end = std::min(x_size, jb + block_size);
73 for (size_t j = jb; j < j_end; ++j) {
74 cov.coeffRef(j, j) = sigma_sq;
75 size_t i_end = std::min(x_size, ib + block_size);
76 for (size_t i = std::max(ib, j + 1); i < i_end; ++i) {
77 return_type_t<T_x> dist = distance(x[i], x[j]);
78 cov.coeffRef(j, i) = cov.coeffRef(i, j)
79 = sigma_sq * (1.0 + root_3_inv_l * dist)
80 * exp(-root_3_inv_l * dist);
81 }
82 }
83 }
84 }
85 return cov;
86 }
87
88 /**
89 * Returns a Matern 3/2 covariance matrix
90 *
91 * \f[ k(x, x') = \sigma^2(1 + \sqrt{3}
92 * \sqrt{\sum_{k=1}^{K}\frac{d(x, x')^2}{l_k^2}})
93 * exp(-\sqrt{3}\sqrt{\sum_{k=1}^{K}\frac{d(x, x')^2}{l_k^2}}) \f]
94 *
95 * where \f$d(x, x')\f$ is the Euclidean distance.
96 *
97 * @tparam T_x type for each scalar
98 * @tparam T_s type of element of parameter sigma
99 * @tparam T_l type of each length scale parameter
100 *
101 * @param x std::vector of Eigen vectors of scalars.
102 * @param length_scale std::vector of length scales
103 * @param sigma standard deviation that can be used in stan::math::square
104 * @throw std::domain error if sigma <= 0, l <= 0, or x is nan or inf
105 */
106 template <typename T_x, typename T_s, typename T_l>
107 inline typename Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic,
108 Eigen::Dynamic>
gp_matern32_cov(const std::vector<Eigen::Matrix<T_x,-1,1>> & x,const T_s & sigma,const std::vector<T_l> & length_scale)109 gp_matern32_cov(const std::vector<Eigen::Matrix<T_x, -1, 1>> &x,
110 const T_s &sigma, const std::vector<T_l> &length_scale) {
111 using std::exp;
112
113 size_t x_size = stan::math::size(x);
114 Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic, Eigen::Dynamic>
115 cov(x_size, x_size);
116
117 if (x_size == 0) {
118 return cov;
119 }
120 const char *function = "gp_matern32_cov";
121 size_t l_size = length_scale.size();
122 for (size_t n = 0; n < x_size; ++n) {
123 check_not_nan(function, "x", x[n]);
124 }
125
126 check_positive_finite(function, "magnitude", sigma);
127 check_positive_finite(function, "length scale", length_scale);
128
129 for (size_t n = 0; n < x_size; ++n) {
130 check_not_nan(function, "length scale", length_scale[n]);
131 }
132
133 check_size_match(function, "x dimension", stan::math::size(x[0]),
134 "number of length scales", l_size);
135
136 T_s sigma_sq = square(sigma);
137 double root_3 = std::sqrt(3.0);
138
139 std::vector<Eigen::Matrix<return_type_t<T_x, T_l>, -1, 1>> x_new
140 = divide_columns(x, length_scale);
141
142 size_t block_size = 10;
143 for (size_t jb = 0; jb < x_size; jb += block_size) {
144 for (size_t ib = jb; ib < x_size; ib += block_size) {
145 size_t j_end = std::min(x_size, jb + block_size);
146 for (size_t j = jb; j < j_end; ++j) {
147 size_t i_end = std::min(x_size, ib + block_size);
148 for (size_t i = std::max(ib, j); i < i_end; ++i) {
149 return_type_t<T_x, T_l> dist = distance(x_new[i], x_new[j]);
150 cov.coeffRef(j, i) = cov.coeffRef(i, j)
151 = sigma_sq * (1.0 + root_3 * dist) * exp(-root_3 * dist);
152 }
153 }
154 }
155 }
156 return cov;
157 }
158
159 /**
160 * Returns a Matern 3/2 cross covariance matrix
161 *
162 * \f[ k(x, x') = \sigma^2(1 +
163 * \frac{\sqrt{3}d(x, x')}{l})exp(-\sqrt{3}\frac{d(x, x')}{l})
164 * \f]
165 *
166 * where \f$d(x, x')\f$ is the Euclidean distance.
167 *
168 * This function is for the cross covariance matrix needed
169 * to compute the posterior predictive density.
170 *
171 * @tparam T_x1 first type of scalars contained in vector x1
172 * @tparam T_x2 second type of scalars contained in vector x2
173 * @tparam T_s type of parameter sigma, marginal standard deviation
174 * @tparam T_l type of parameter length scale
175 *
176 * @param x1 std::vector of scalars that can be used in squared_distance
177 * @param x2 std::vector of scalars that can be used in squared_distance
178 * @param length_scale length scale
179 * @param sigma standard deviation that can be used in stan::math::square
180 * @throw std::domain error if sigma <= 0, l <= 0, or x1, x2 are nan or inf
181 *
182 */
183 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
184 inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>,
185 Eigen::Dynamic, Eigen::Dynamic>
gp_matern32_cov(const std::vector<T_x1> & x1,const std::vector<T_x2> & x2,const T_s & sigma,const T_l & length_scale)186 gp_matern32_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
187 const T_s &sigma, const T_l &length_scale) {
188 using std::exp;
189
190 size_t x1_size = stan::math::size(x1);
191 size_t x2_size = stan::math::size(x2);
192 Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>, Eigen::Dynamic,
193 Eigen::Dynamic>
194 cov(x1_size, x2_size);
195
196 if (x1_size == 0 || x2_size == 0) {
197 return cov;
198 }
199
200 const char *function = "gp_matern32_cov";
201 size_t x1_obs_size = stan::math::size(x1[0]);
202 for (size_t i = 0; i < x1_size; ++i) {
203 check_size_match(function, "x1's row", x1_obs_size, "x1's other row",
204 stan::math::size(x1[i]));
205 }
206 for (size_t i = 0; i < x2_size; ++i) {
207 check_size_match(function, "x1's row", x1_obs_size, "x2's other row",
208 stan::math::size(x2[i]));
209 }
210
211 for (size_t n = 0; n < x1_size; ++n) {
212 check_not_nan(function, "x1", x1[n]);
213 }
214 for (size_t n = 0; n < x2_size; ++n) {
215 check_not_nan(function, "x2", x2[n]);
216 }
217
218 check_positive_finite(function, "magnitude", sigma);
219 check_positive_finite(function, "length scale", length_scale);
220
221 T_s sigma_sq = square(sigma);
222 T_l root_3_inv_l_sq = std::sqrt(3.0) / length_scale;
223
224 size_t block_size = 10;
225 for (size_t ib = 0; ib < x1.size(); ib += block_size) {
226 for (size_t jb = 0; jb < x2.size(); jb += block_size) {
227 size_t j_end = std::min(x2_size, jb + block_size);
228 for (size_t j = jb; j < j_end; ++j) {
229 size_t i_end = std::min(x1_size, ib + block_size);
230 for (size_t i = ib; i < i_end; ++i) {
231 return_type_t<T_x1, T_x2> dist = distance(x1[i], x2[j]);
232 cov(i, j) = sigma_sq * (1.0 + root_3_inv_l_sq * dist)
233 * exp(-root_3_inv_l_sq * dist);
234 }
235 }
236 }
237 }
238 return cov;
239 }
240
241 /**
242 * Returns a Matern 3/2 cross covariance matrix
243 *
244 * \f[ k(x, x') = \sigma^2(1 + \sqrt{3}
245 * \sqrt{\sum_{k=1}^{K}\frac{d(x, x')^2}{l_k^2}})
246 * exp(-\sqrt{3}\sqrt{\sum_{k=1}^{K}\frac{d(x, x')^2}{l_k^2}}) \f]
247 *
248 * where \f$d(x, x')\f$ is the Euclidean distance
249 *
250 * This function is for the cross covariance matrix needed
251 * to compute the posterior predictive density.
252 *
253 * @tparam T_x1 first type of std::vector of scalars
254 * @tparam T_x2 second type of std::vector of scalars
255 * @tparam T_s type of parameter sigma, marginal standard deviation
256 * @tparam T_l type of parameter length scale
257 *
258 * @param x1 std::vector of Eigen vectors of scalars
259 * @param x2 std::vector of Eigen vectors of scalars
260 * @param length_scale parameter length scale
261 * @param sigma standard deviation that can be used in stan::math::square
262 * @throw std::domain error if sigma <= 0, l <= 0, or x1, x2 are nan or inf
263 *
264 */
265 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
266 inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>,
267 Eigen::Dynamic, Eigen::Dynamic>
gp_matern32_cov(const std::vector<Eigen::Matrix<T_x1,-1,1>> & x1,const std::vector<Eigen::Matrix<T_x2,-1,1>> & x2,const T_s & sigma,const std::vector<T_l> & length_scale)268 gp_matern32_cov(const std::vector<Eigen::Matrix<T_x1, -1, 1>> &x1,
269 const std::vector<Eigen::Matrix<T_x2, -1, 1>> &x2,
270 const T_s &sigma, const std::vector<T_l> &length_scale) {
271 using std::exp;
272
273 size_t x1_size = stan::math::size(x1);
274 size_t x2_size = stan::math::size(x2);
275 Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>, Eigen::Dynamic,
276 Eigen::Dynamic>
277 cov(x1_size, x2_size);
278
279 if (x1_size == 0 || x2_size == 0) {
280 return cov;
281 }
282
283 const char *function = "gp_matern_32_cov";
284 for (size_t n = 0; n < x1_size; ++n) {
285 check_not_nan(function, "x1", x1[n]);
286 }
287 for (size_t n = 0; n < x2_size; ++n) {
288 check_not_nan(function, "x2", x2[n]);
289 }
290
291 check_positive_finite(function, "magnitude", sigma);
292 check_positive_finite(function, "length scale", length_scale);
293
294 size_t l_size = length_scale.size();
295 for (size_t n = 0; n < l_size; ++n) {
296 check_not_nan(function, "length scale", length_scale[n]);
297 }
298
299 for (size_t i = 0; i < x1_size; ++i) {
300 check_size_match(function, "x1's row", stan::math::size(x1[i]),
301 "number of length scales", l_size);
302 }
303 for (size_t i = 0; i < x2_size; ++i) {
304 check_size_match(function, "x2's row", stan::math::size(x2[i]),
305 "number of length scales", l_size);
306 }
307
308 T_s sigma_sq = square(sigma);
309 double root_3 = std::sqrt(3.0);
310
311 std::vector<Eigen::Matrix<return_type_t<T_x1, T_l>, -1, 1>> x1_new
312 = divide_columns(x1, length_scale);
313 std::vector<Eigen::Matrix<return_type_t<T_x2, T_l>, -1, 1>> x2_new
314 = divide_columns(x2, length_scale);
315
316 size_t block_size = 10;
317
318 for (size_t ib = 0; ib < x1.size(); ib += block_size) {
319 for (size_t jb = 0; jb < x2.size(); jb += block_size) {
320 size_t j_end = std::min(x2_size, jb + block_size);
321 for (size_t j = jb; j < j_end; ++j) {
322 size_t i_end = std::min(x1_size, ib + block_size);
323 for (size_t i = ib; i < i_end; ++i) {
324 return_type_t<T_x1, T_x2, T_l> dist = distance(x1_new[i], x2_new[j]);
325 cov(i, j) = sigma_sq * (1.0 + root_3 * dist) * exp(-root_3 * dist);
326 }
327 }
328 }
329 }
330 return cov;
331 }
332 } // namespace math
333 } // namespace stan
334 #endif
335