Stan Math Library  2.9.0
reverse mode automatic differentiation
gaussian_dlm_obs_log.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_PROB_GAUSSIAN_DLM_OBS_LOG_HPP
2 #define STAN_MATH_PRIM_MAT_PROB_GAUSSIAN_DLM_OBS_LOG_HPP
3 
4 #include <boost/random/normal_distribution.hpp>
5 #include <boost/random/variate_generator.hpp>
24 
27 
28 /*
29  TODO: time-varying system matrices
30  TODO: use sequential processing even for non-diagonal obs
31  covariance.
32  TODO: add constant terms in observation.
33 */
34 
35 namespace stan {
36  namespace math {
70  template <bool propto,
71  typename T_y,
72  typename T_F, typename T_G,
73  typename T_V, typename T_W,
74  typename T_m0, typename T_C0
75  >
76  typename return_type<
77  T_y,
78  typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type >::type
79  gaussian_dlm_obs_log(const Eigen::Matrix
80  <T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
81  const Eigen::Matrix
82  <T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
83  const Eigen::Matrix
84  <T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
85  const Eigen::Matrix
86  <T_V, Eigen::Dynamic, Eigen::Dynamic>& V,
87  const Eigen::Matrix
88  <T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
89  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
90  const Eigen::Matrix
91  <T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
92  static const char* function("stan::math::gaussian_dlm_obs_log");
93  typedef typename return_type<
94  T_y,
96  T_lp lp(0.0);
97 
98  using stan::math::add;
107  using stan::math::multiply;
109  using stan::math::subtract;
111  using stan::math::transpose;
112 
113  int r = y.rows(); // number of variables
114  int T = y.cols(); // number of observations
115  int n = G.rows(); // number of states
116 
117  // check y
118  check_finite(function, "y", y);
119  check_not_nan(function, "y", y);
120  // check F
121  check_size_match(function,
122  "columns of F", F.cols(),
123  "rows of y", y.rows());
124  check_size_match(function,
125  "rows of F", F.rows(),
126  "rows of G", G.rows());
127  check_finite(function, "F", F);
128  // check G
129  check_square(function, "G", G);
130  check_finite(function, "G", G);
131  // check V
132  check_size_match(function,
133  "rows of V", V.rows(),
134  "rows of y", y.rows());
135  // TODO(anyone): incorporate support for infinite V
136  check_finite(function, "V", V);
137  check_spsd_matrix(function, "V", V);
138  // check W
139  check_size_match(function,
140  "rows of W", W.rows(),
141  "rows of G", G.rows());
142  // TODO(anyone): incorporate support for infinite W
143  check_finite(function, "W", W);
144  check_spsd_matrix(function, "W", W);
145  // check m0
146  check_size_match(function,
147  "size of m0", m0.size(),
148  "rows of G", G.rows());
149  check_finite(function, "m0", m0);
150  // check C0
151  check_size_match(function,
152  "rows of C0", C0.rows(),
153  "rows of G", G.rows());
154  check_cov_matrix(function, "C0", C0);
155  check_finite(function, "C0", C0);
156 
157  if (y.cols() == 0 || y.rows() == 0)
158  return lp;
159 
161  lp -= 0.5 * LOG_TWO_PI * r * T;
162  }
163 
165  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> m(n);
166  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
167 
168  // TODO(anyone): how to recast matrices
169  for (int i = 0; i < m0.size(); i++) {
170  m(i) = m0(i);
171  }
172  for (int i = 0; i < C0.rows(); i++) {
173  for (int j = 0; j < C0.cols(); j++) {
174  C(i, j) = C0(i, j);
175  }
176  }
177 
178  Eigen::Matrix<typename return_type<T_y>::type,
179  Eigen::Dynamic, 1> yi(r);
180  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> a(n);
181  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> R(n, n);
182  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> f(r);
183  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> Q(r, r);
184  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> Q_inv(r, r);
185  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> e(r);
186  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> A(n, r);
187 
188  for (int i = 0; i < y.cols(); i++) {
189  yi = y.col(i);
190  // // Predict state
191  // a_t = G_t m_{t-1}
192  a = multiply(G, m);
193  // R_t = G_t C_{t-1} G_t' + W_t
194  R = add(quad_form_sym(C, transpose(G)), W);
195  // // predict observation
196  // f_t = F_t' a_t
197  f = multiply(transpose(F), a);
198  // Q_t = F'_t R_t F_t + V_t
199  Q = add(quad_form_sym(R, F), V);
200  Q_inv = inverse_spd(Q);
201  // // filtered state
202  // e_t = y_t - f_t
203  e = subtract(yi, f);
204  // A_t = R_t F_t Q^{-1}_t
205  A = multiply(multiply(R, F), Q_inv);
206  // m_t = a_t + A_t e_t
207  m = add(a, multiply(A, e));
208  // C = R_t - A_t Q_t A_t'
209  C = subtract(R, quad_form_sym(Q, transpose(A)));
210  lp -= 0.5 * (log_determinant_spd(Q) + trace_quad_form(Q_inv, e));
211  }
212  }
213  return lp;
214  }
215 
216  template <typename T_y,
217  typename T_F, typename T_G,
218  typename T_V, typename T_W,
219  typename T_m0, typename T_C0
220  >
221  inline
222  typename return_type<
223  T_y,
225  gaussian_dlm_obs_log(const Eigen::Matrix
226  <T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
227  const Eigen::Matrix
228  <T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
229  const Eigen::Matrix
230  <T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
231  const Eigen::Matrix
232  <T_V, Eigen::Dynamic, Eigen::Dynamic>& V,
233  const Eigen::Matrix
234  <T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
235  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
236  const Eigen::Matrix
237  <T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
238  return gaussian_dlm_obs_log<false>(y, F, G, V, W, m0, C0);
239  }
240 
276  template <bool propto,
277  typename T_y,
278  typename T_F, typename T_G,
279  typename T_V, typename T_W,
280  typename T_m0, typename T_C0
281  >
282  typename return_type<
283  T_y,
285  gaussian_dlm_obs_log(const Eigen::Matrix
286  <T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
287  const Eigen::Matrix
288  <T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
289  const Eigen::Matrix
290  <T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
291  const Eigen::Matrix<T_V, Eigen::Dynamic, 1>& V,
292  const Eigen::Matrix
293  <T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
294  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
295  const Eigen::Matrix
296  <T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
297  static const char* function("stan::math::gaussian_dlm_obs_log");
298  typedef
299  typename return_type
301  T_lp;
302  T_lp lp(0.0);
303 
304  using stan::math::add;
312  using stan::math::multiply;
314  using stan::math::subtract;
317  using stan::math::transpose;
318  using std::log;
319 
320  int r = y.rows(); // number of variables
321  int T = y.cols(); // number of observations
322  int n = G.rows(); // number of states
323 
324  // check y
325  check_finite(function, "y", y);
326  check_not_nan(function, "y", y);
327  // check F
328  check_size_match(function,
329  "columns of F", F.cols(),
330  "rows of y", y.rows());
331  check_size_match(function,
332  "rows of F", F.rows(),
333  "rows of G", G.rows());
334  check_finite(function, "F", F);
335  check_not_nan(function, "F", F);
336  // check G
337  check_size_match(function,
338  "rows of G", G.rows(),
339  "columns of G", G.cols());
340  check_finite(function, "G", G);
341  check_not_nan(function, "G", G);
342  // check V
343  check_nonnegative(function, "V", V);
344  check_size_match(function,
345  "size of V", V.size(),
346  "rows of y", y.rows());
347  // TODO(anyone): support infinite V
348  check_finite(function, "V", V);
349  check_not_nan(function, "V", V);
350  // check W
351  check_spsd_matrix(function, "W", W);
352  check_size_match(function,
353  "rows of W", W.rows(),
354  "rows of G", G.rows());
355  // TODO(anyone): support infinite W
356  check_finite(function, "W", W);
357  check_not_nan(function, "W", W);
358  // check m0
359  check_size_match(function,
360  "size of m0", m0.size(),
361  "rows of G", G.rows());
362  check_finite(function, "m0", m0);
363  check_not_nan(function, "m0", m0);
364  // check C0
365  check_cov_matrix(function, "C0", C0);
366  check_size_match(function,
367  "rows of C0", C0.rows(),
368  "rows of G", G.rows());
369  check_finite(function, "C0", C0);
370  check_not_nan(function, "C0", C0);
371 
372  if (y.cols() == 0 || y.rows() == 0)
373  return lp;
374 
376  lp += 0.5 * NEG_LOG_TWO_PI * r * T;
377  }
378 
380  T_lp f;
381  T_lp Q;
382  T_lp Q_inv;
383  T_lp e;
384  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> A(n);
385  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> Fj(n);
386  Eigen::Matrix<T_lp, Eigen::Dynamic, 1> m(n);
387  Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
388 
389  // TODO(anyone): how to recast matrices
390  for (int i = 0; i < m0.size(); i++) {
391  m(i) = m0(i);
392  }
393  for (int i = 0; i < C0.rows(); i++) {
394  for (int j = 0; j < C0.cols(); j++) {
395  C(i, j) = C0(i, j);
396  }
397  }
398 
399  for (int i = 0; i < y.cols(); i++) {
400  // Predict state
401  // reuse m and C instead of using a and R
402  m = multiply(G, m);
403  C = add(quad_form_sym(C, transpose(G)), W);
404  for (int j = 0; j < y.rows(); ++j) {
405  // predict observation
406  T_lp yij(y(j, i));
407  // dim Fj = (n, 1)
408  for (int k = 0; k < F.rows(); ++k) {
409  Fj(k) = F(k, j);
410  }
411  // // f_{t, i} = F_{t, i}' m_{t, i-1}
412  f = dot_product(Fj, m);
413  Q = trace_quad_form(C, Fj) + V(j);
414  Q_inv = 1.0 / Q;
415  // // filtered observation
416  // // e_{t, i} = y_{t, i} - f_{t, i}
417  e = yij - f;
418  // // A_{t, i} = C_{t, i-1} F_{t, i} Q_{t, i}^{-1}
419  A = multiply(multiply(C, Fj), Q_inv);
420  // // m_{t, i} = m_{t, i-1} + A_{t, i} e_{t, i}
421  m += multiply(A, e);
422  // // c_{t, i} = C_{t, i-1} - Q_{t, i} A_{t, i} A_{t, i}'
423  // // // tcrossprod throws an error (ambiguous)
424  // C = subtract(C, multiply(Q, tcrossprod(A)));
425  C -= multiply(Q, multiply(A, transpose(A)));
426  C = 0.5 * add(C, transpose(C));
427  lp -= 0.5 * (log(Q) + pow(e, 2) * Q_inv);
428  }
429  }
430  }
431  return lp;
432  }
433 
434  template <typename T_y,
435  typename T_F, typename T_G,
436  typename T_V, typename T_W,
437  typename T_m0, typename T_C0>
438  inline
439  typename return_type
442  (const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
443  const Eigen::Matrix<T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
444  const Eigen::Matrix<T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
445  const Eigen::Matrix<T_V, Eigen::Dynamic, 1>& V,
446  const Eigen::Matrix<T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
447  const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
448  const Eigen::Matrix<T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
449  return gaussian_dlm_obs_log<false>(y, F, G, V, W, m0, C0);
450  }
451  }
452 
453 }
454 
455 #endif
fvar< T > trace_quad_form(const Eigen::Matrix< fvar< T >, RA, CA > &A, const Eigen::Matrix< fvar< T >, RB, CB > &B)
bool check_not_nan(const char *function, const char *name, const T_y &y)
Return true if y is not NaN.
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R, C > subtract(const Eigen::Matrix< T1, R, C > &m1, const Eigen::Matrix< T2, R, C > &m2)
Return the result of subtracting the second specified matrix from the first specified matrix...
Definition: subtract.hpp:27
return_type< T_y, typename return_type< T_F, T_G, T_V, T_W, T_m0, T_C0 >::type >::type gaussian_dlm_obs_log(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).
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:15
Eigen::Matrix< fvar< T >, R1, C1 > multiply(const Eigen::Matrix< fvar< T >, R1, C1 > &m, const fvar< T > &c)
Definition: multiply.hpp:20
Metaprogram to calculate the base scalar return type resulting from promoting all the scalar types of...
Definition: return_type.hpp:19
Eigen::Matrix< fvar< T >, R, R > tcrossprod(const Eigen::Matrix< fvar< T >, R, C > &m)
Definition: tcrossprod.hpp:17
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
bool check_spsd_matrix(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified matrix is a square, symmetric, and positive semi-definite.
boost::math::tools::promote_args< typename scalar_type< T1 >::type, typename scalar_type< T2 >::type, typename scalar_type< T3 >::type, typename scalar_type< T4 >::type, typename scalar_type< T5 >::type, typename scalar_type< T6 >::type >::type type
Definition: return_type.hpp:27
const double LOG_TWO_PI
Definition: constants.hpp:193
bool check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Return true if the provided sizes match.
fvar< T > dot_product(const Eigen::Matrix< fvar< T >, R1, C1 > &v1, const Eigen::Matrix< fvar< T >, R2, C2 > &v2)
Definition: dot_product.hpp:20
bool check_cov_matrix(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified matrix is a valid covariance matrix.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95
bool check_finite(const char *function, const char *name, const T_y &y)
Return true if y is finite.
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:19
const double NEG_LOG_TWO_PI
Definition: constants.hpp:195
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R, C > add(const Eigen::Matrix< T1, R, C > &m1, const Eigen::Matrix< T2, R, C > &m2)
Return the sum of the specified matrices.
Definition: add.hpp:27
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:18
bool check_nonnegative(const char *function, const char *name, const T_y &y)
Return true if y is non-negative.
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
bool check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified matrix is square.
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–2015, Stan Development Team.