Stan Math Library  2.20.0
reverse mode automatic differentiation
gaussian_dlm_obs_lpdf.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_PROB_GAUSSIAN_DLM_OBS_LPDF_HPP
2 #define STAN_MATH_PRIM_MAT_PROB_GAUSSIAN_DLM_OBS_LPDF_HPP
3 
24 
25 /*
26  TODO: time-varying system matrices
27  TODO: use sequential processing even for non-diagonal obs
28  covariance.
29  TODO: add constant terms in observation.
30 */
31 namespace stan {
32 namespace math {
66 template <bool propto, typename T_y, typename T_F, typename T_G, typename T_V,
67  typename T_W, typename T_m0, typename T_C0>
68 typename return_type<
71  const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
72  const Eigen::Matrix<T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
73  const Eigen::Matrix<T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
74  const Eigen::Matrix<T_V, Eigen::Dynamic, Eigen::Dynamic>& V,
75  const Eigen::Matrix<T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
76  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
77  const Eigen::Matrix<T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
78  static const char* function = "gaussian_dlm_obs_lpdf";
79  typedef
80  typename return_type<T_y, typename return_type<T_F, T_G, T_V, T_W, T_m0,
81  T_C0>::type>::type T_lp;
82  int r = y.rows(); // number of variables
83  int T = y.cols(); // number of observations
84  int n = G.rows(); // number of states
85 
86  check_finite(function, "y", y);
87  check_not_nan(function, "y", y);
88  check_size_match(function, "columns of F", F.cols(), "rows of y", y.rows());
89  check_size_match(function, "rows of F", F.rows(), "rows of G", G.rows());
90  check_finite(function, "F", F);
91  check_square(function, "G", G);
92  check_finite(function, "G", G);
93  check_size_match(function, "rows of V", V.rows(), "rows of y", y.rows());
94  // TODO(anyone): incorporate support for infinite V
95  check_finite(function, "V", V);
96  check_spsd_matrix(function, "V", V);
97  check_size_match(function, "rows of W", W.rows(), "rows of G", G.rows());
98  // TODO(anyone): incorporate support for infinite W
99  check_finite(function, "W", W);
100  check_spsd_matrix(function, "W", W);
101  check_size_match(function, "size of m0", m0.size(), "rows of G", G.rows());
102  check_finite(function, "m0", m0);
103  check_size_match(function, "rows of C0", C0.rows(), "rows of G", G.rows());
104  check_pos_definite(function, "C0", C0);
105  check_finite(function, "C0", C0);
106 
107  if (size_zero(y))
108  return 0;
109 
110  T_lp lp(0);
112  lp -= 0.5 * LOG_TWO_PI * r * T;
113  }
114 
116  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> m(n);
117  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
118 
119  // TODO(anyone): how to recast matrices
120  for (int i = 0; i < m0.size(); i++) {
121  m(i) = m0(i);
122  }
123  for (int i = 0; i < C0.rows(); i++) {
124  for (int j = 0; j < C0.cols(); j++) {
125  C(i, j) = C0(i, j);
126  }
127  }
128 
129  Eigen::Matrix<typename return_type<T_y>::type, Eigen::Dynamic, 1> yi(r);
130  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> a(n);
131  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> R(n, n);
132  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> f(r);
133  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> Q(r, r);
134  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> Q_inv(r, r);
135  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> e(r);
136  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> A(n, r);
137 
138  for (int i = 0; i < y.cols(); i++) {
139  yi = y.col(i);
140  // // Predict state
141  // a_t = G_t m_{t-1}
142  a = multiply(G, m);
143  // R_t = G_t C_{t-1} G_t' + W_t
144  R = add(quad_form_sym(C, transpose(G)), W);
145  // // predict observation
146  // f_t = F_t' a_t
147  f = multiply(transpose(F), a);
148  // Q_t = F'_t R_t F_t + V_t
149  Q = add(quad_form_sym(R, F), V);
150  Q_inv = inverse_spd(Q);
151  // // filtered state
152  // e_t = y_t - f_t
153  e = subtract(yi, f);
154  // A_t = R_t F_t Q^{-1}_t
155  A = multiply(multiply(R, F), Q_inv);
156  // m_t = a_t + A_t e_t
157  m = add(a, multiply(A, e));
158  // C = R_t - A_t Q_t A_t'
159  C = subtract(R, quad_form_sym(Q, transpose(A)));
160  lp -= 0.5 * (log_determinant_spd(Q) + trace_quad_form(Q_inv, e));
161  }
162  }
163  return lp;
164 }
165 
166 template <typename T_y, typename T_F, typename T_G, typename T_V, typename T_W,
167  typename T_m0, typename T_C0>
168 inline typename return_type<
171  const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
172  const Eigen::Matrix<T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
173  const Eigen::Matrix<T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
174  const Eigen::Matrix<T_V, Eigen::Dynamic, Eigen::Dynamic>& V,
175  const Eigen::Matrix<T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
176  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
177  const Eigen::Matrix<T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
178  return gaussian_dlm_obs_lpdf<false>(y, F, G, V, W, m0, C0);
179 }
180 
216 template <bool propto, typename T_y, typename T_F, typename T_G, typename T_V,
217  typename T_W, typename T_m0, typename T_C0>
218 typename return_type<
221  const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
222  const Eigen::Matrix<T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
223  const Eigen::Matrix<T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
224  const Eigen::Matrix<T_V, Eigen::Dynamic, 1>& V,
225  const Eigen::Matrix<T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
226  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
227  const Eigen::Matrix<T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
228  static const char* function = "gaussian_dlm_obs_lpdf";
229  typedef
230  typename return_type<T_y, typename return_type<T_F, T_G, T_V, T_W, T_m0,
231  T_C0>::type>::type T_lp;
232  using std::log;
233 
234  int r = y.rows(); // number of variables
235  int T = y.cols(); // number of observations
236  int n = G.rows(); // number of states
237 
238  check_finite(function, "y", y);
239  check_not_nan(function, "y", y);
240  check_size_match(function, "columns of F", F.cols(), "rows of y", y.rows());
241  check_size_match(function, "rows of F", F.rows(), "rows of G", G.rows());
242  check_finite(function, "F", F);
243  check_not_nan(function, "F", F);
244  check_size_match(function, "rows of G", G.rows(), "columns of G", G.cols());
245  check_finite(function, "G", G);
246  check_not_nan(function, "G", G);
247  check_nonnegative(function, "V", V);
248  check_size_match(function, "size of V", V.size(), "rows of y", y.rows());
249  // TODO(anyone): support infinite V
250  check_finite(function, "V", V);
251  check_not_nan(function, "V", V);
252  check_spsd_matrix(function, "W", W);
253  check_size_match(function, "rows of W", W.rows(), "rows of G", G.rows());
254  // TODO(anyone): support infinite W
255  check_finite(function, "W", W);
256  check_not_nan(function, "W", W);
257  check_size_match(function, "size of m0", m0.size(), "rows of G", G.rows());
258  check_finite(function, "m0", m0);
259  check_not_nan(function, "m0", m0);
260  check_pos_definite(function, "C0", C0);
261  check_size_match(function, "rows of C0", C0.rows(), "rows of G", G.rows());
262  check_finite(function, "C0", C0);
263  check_not_nan(function, "C0", C0);
264 
265  if (y.cols() == 0 || y.rows() == 0)
266  return 0;
267 
268  T_lp lp(0);
270  lp += 0.5 * NEG_LOG_TWO_PI * r * T;
271  }
272 
274  T_lp f;
275  T_lp Q;
276  T_lp Q_inv;
277  T_lp e;
278  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> A(n);
279  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> Fj(n);
280  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> m(n);
281  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
282 
283  // TODO(anyone): how to recast matrices
284  for (int i = 0; i < m0.size(); i++) {
285  m(i) = m0(i);
286  }
287  for (int i = 0; i < C0.rows(); i++) {
288  for (int j = 0; j < C0.cols(); j++) {
289  C(i, j) = C0(i, j);
290  }
291  }
292 
293  for (int i = 0; i < y.cols(); i++) {
294  // Predict state
295  // reuse m and C instead of using a and R
296  m = multiply(G, m);
297  C = add(quad_form_sym(C, transpose(G)), W);
298  for (int j = 0; j < y.rows(); ++j) {
299  // predict observation
300  T_lp yij(y(j, i));
301  // dim Fj = (n, 1)
302  for (int k = 0; k < F.rows(); ++k) {
303  Fj(k) = F(k, j);
304  }
305  // f_{t, i} = F_{t, i}' m_{t, i-1}
306  f = dot_product(Fj, m);
307  Q = trace_quad_form(C, Fj) + V(j);
308  Q_inv = 1.0 / Q;
309  // filtered observation
310  // e_{t, i} = y_{t, i} - f_{t, i}
311  e = yij - f;
312  // A_{t, i} = C_{t, i-1} F_{t, i} Q_{t, i}^{-1}
313  A = multiply(multiply(C, Fj), Q_inv);
314  // m_{t, i} = m_{t, i-1} + A_{t, i} e_{t, i}
315  m += multiply(A, e);
316  // c_{t, i} = C_{t, i-1} - Q_{t, i} A_{t, i} A_{t, i}'
317  // tcrossprod throws an error (ambiguous)
318  // C = subtract(C, multiply(Q, tcrossprod(A)));
319  C -= multiply(Q, multiply(A, transpose(A)));
320  C = 0.5 * add(C, transpose(C));
321  lp -= 0.5 * (log(Q) + pow(e, 2) * Q_inv);
322  }
323  }
324  }
325  return lp;
326 }
327 
328 template <typename T_y, typename T_F, typename T_G, typename T_V, typename T_W,
329  typename T_m0, typename T_C0>
330 inline typename return_type<
333  const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
334  const Eigen::Matrix<T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
335  const Eigen::Matrix<T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
336  const Eigen::Matrix<T_V, Eigen::Dynamic, 1>& V,
337  const Eigen::Matrix<T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
338  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
339  const Eigen::Matrix<T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
340  return gaussian_dlm_obs_lpdf<false>(y, F, G, V, W, m0, C0);
341 }
342 
343 } // namespace math
344 } // namespace stan
345 #endif
auto subtract(const matrix_cl &A, const matrix_cl &B)
Matrix subtraction on the OpenCL device Subtracts the second matrix from the first matrix and stores ...
Definition: subtract.hpp:28
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
fvar< T > trace_quad_form(const Eigen::Matrix< fvar< T >, RA, CA > &A, const Eigen::Matrix< fvar< T >, RB, CB > &B)
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
void check_square(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is square.
Eigen::Matrix< fvar< T >, R1, C1 > multiply(const Eigen::Matrix< fvar< T >, R1, C1 > &m, const fvar< T > &c)
Definition: multiply.hpp:14
Template metaprogram to calculate the base scalar return type resulting from promoting all the scalar...
Definition: return_type.hpp:33
bool size_zero(T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition: size_zero.hpp:18
return_type< T_y, typename return_type< T_F, T_G, T_V, T_W, T_m0, T_C0 >::type >::type gaussian_dlm_obs_lpdf(const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y, const Eigen::Matrix< T_F, Eigen::Dynamic, Eigen::Dynamic > &F, const Eigen::Matrix< T_G, Eigen::Dynamic, Eigen::Dynamic > &G, const Eigen::Matrix< T_V, Eigen::Dynamic, Eigen::Dynamic > &V, const Eigen::Matrix< T_W, Eigen::Dynamic, Eigen::Dynamic > &W, const Eigen::Matrix< T_m0, Eigen::Dynamic, 1 > &m0, const Eigen::Matrix< T_C0, Eigen::Dynamic, Eigen::Dynamic > &C0)
The log of a Gaussian dynamic linear model (GDLM).
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.
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
const double LOG_TWO_PI
Definition: constants.hpp:164
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
Definition: return_type.hpp:36
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
void check_spsd_matrix(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Check if the specified matrix is a square, symmetric, and positive semi-definite. ...
fvar< T > dot_product(const Eigen::Matrix< fvar< T >, R1, C1 > &v1, const Eigen::Matrix< fvar< T >, R2, C2 > &v2)
Definition: dot_product.hpp:14
matrix_cl transpose(const matrix_cl &src)
Takes the transpose of the matrix on the OpenCL device.
Definition: transpose.hpp:20
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
matrix_cl add(const matrix_cl &A, const matrix_cl &B)
Matrix addition on the OpenCL device.
Definition: add.hpp:24
Eigen::Matrix< fvar< T >, CB, CB > quad_form_sym(const Eigen::Matrix< fvar< T >, RA, CA > &A, const Eigen::Matrix< double, RB, CB > &B)
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > inverse_spd(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Returns the inverse of the specified symmetric, pos/neg-definite matrix.
Definition: inverse_spd.hpp:18
const double NEG_LOG_TWO_PI
Definition: constants.hpp:166
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16
void check_pos_definite(const char *function, const char *name, const Eigen::Matrix< T_y, -1, -1 > &y)
Check if the specified square, symmetric matrix is positive definite.
T log_determinant_spd(const Eigen::Matrix< T, R, C > &m)
Returns the log absolute determinant of the specified square matrix.

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