| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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>
6 #include <stan/math/prim/mat/err/check_cov_matrix.hpp>
7 #include <stan/math/prim/scal/err/check_size_match.hpp>
8 #include <stan/math/prim/mat/err/check_spsd_matrix.hpp>
9 #include <stan/math/prim/scal/err/check_finite.hpp>
10 #include <stan/math/prim/scal/err/check_nonnegative.hpp>
11 #include <stan/math/prim/scal/err/check_not_nan.hpp>
12 #include <stan/math/prim/mat/fun/add.hpp>
13 #include <stan/math/prim/mat/fun/dot_product.hpp>
14 #include <stan/math/prim/mat/fun/inverse_spd.hpp>
15 #include <stan/math/prim/mat/fun/log.hpp>
16 #include <stan/math/prim/mat/fun/log_determinant_spd.hpp>
17 #include <stan/math/prim/mat/fun/multiply.hpp>
18 #include <stan/math/prim/mat/fun/quad_form.hpp>
19 #include <stan/math/prim/mat/fun/quad_form_sym.hpp>
20 #include <stan/math/prim/mat/fun/subtract.hpp>
21 #include <stan/math/prim/mat/fun/tcrossprod.hpp>
22 #include <stan/math/prim/mat/fun/trace_quad_form.hpp>
23 #include <stan/math/prim/mat/fun/transpose.hpp>
24 #include <stan/math/prim/scal/fun/constants.hpp>
25 #include <stan/math/prim/scal/meta/include_summand.hpp>
26
27 /*
28 TODO: time-varying system matrices
29 TODO: use sequential processing even for non-diagonal obs
30 covariance.
31 TODO: add constant terms in observation.
32 */
33
34 namespace stan {
35 namespace math {
69 template <bool propto,
70 typename T_y,
71 typename T_F, typename T_G,
72 typename T_V, typename T_W,
73 typename T_m0, typename T_C0
74 >
75 typename return_type<
76 T_y,
77 typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type >::type
78 gaussian_dlm_obs_log(const Eigen::Matrix
79 <T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
80 const Eigen::Matrix
81 <T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
82 const Eigen::Matrix
83 <T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
84 const Eigen::Matrix
85 <T_V, Eigen::Dynamic, Eigen::Dynamic>& V,
86 const Eigen::Matrix
87 <T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
88 const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
89 const Eigen::Matrix
90 <T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
91 static const char* function("gaussian_dlm_obs_log");
92 typedef typename return_type<
93 T_y,
94 typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type>::type T_lp;
95 T_lp lp(0.0);
96
97 int r = y.rows(); // number of variables
98 int T = y.cols(); // number of observations
99 int n = G.rows(); // number of states
100
101 check_finite(function, "y", y);
102 check_not_nan(function, "y", y);
103 check_size_match(function,
104 "columns of F", F.cols(),
105 "rows of y", y.rows());
106 check_size_match(function,
107 "rows of F", F.rows(),
108 "rows of G", G.rows());
109 check_finite(function, "F", F);
110 check_square(function, "G", G);
111 check_finite(function, "G", G);
112 check_size_match(function,
113 "rows of V", V.rows(),
114 "rows of y", y.rows());
115 // TODO(anyone): incorporate support for infinite V
116 check_finite(function, "V", V);
117 check_spsd_matrix(function, "V", V);
118 check_size_match(function,
119 "rows of W", W.rows(),
120 "rows of G", G.rows());
121 // TODO(anyone): incorporate support for infinite W
122 check_finite(function, "W", W);
123 check_spsd_matrix(function, "W", W);
124 check_size_match(function,
125 "size of m0", m0.size(),
126 "rows of G", G.rows());
127 check_finite(function, "m0", m0);
128 check_size_match(function,
129 "rows of C0", C0.rows(),
130 "rows of G", G.rows());
131 check_cov_matrix(function, "C0", C0);
132 check_finite(function, "C0", C0);
133
134 if (y.cols() == 0 || y.rows() == 0)
135 return lp;
136
137 if (include_summand<propto>::value) {
138 lp -= 0.5 * LOG_TWO_PI * r * T;
139 }
140
141 if (include_summand<propto, T_y, T_F, T_G, T_V, T_W, T_m0, T_C0>::value) {
142 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> m(n);
143 Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
144
145 // TODO(anyone): how to recast matrices
146 for (int i = 0; i < m0.size(); i++) {
147 m(i) = m0(i);
148 }
149 for (int i = 0; i < C0.rows(); i++) {
150 for (int j = 0; j < C0.cols(); j++) {
151 C(i, j) = C0(i, j);
152 }
153 }
154
155 Eigen::Matrix<typename return_type<T_y>::type,
156 Eigen::Dynamic, 1> yi(r);
157 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> a(n);
158 Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> R(n, n);
159 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> f(r);
160 Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> Q(r, r);
161 Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> Q_inv(r, r);
162 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> e(r);
163 Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> A(n, r);
164
165 for (int i = 0; i < y.cols(); i++) {
166 yi = y.col(i);
167 // // Predict state
168 // a_t = G_t m_{t-1}
169 a = multiply(G, m);
170 // R_t = G_t C_{t-1} G_t' + W_t
171 R = add(quad_form_sym(C, transpose(G)), W);
172 // // predict observation
173 // f_t = F_t' a_t
174 f = multiply(transpose(F), a);
175 // Q_t = F'_t R_t F_t + V_t
176 Q = add(quad_form_sym(R, F), V);
177 Q_inv = inverse_spd(Q);
178 // // filtered state
179 // e_t = y_t - f_t
180 e = subtract(yi, f);
181 // A_t = R_t F_t Q^{-1}_t
182 A = multiply(multiply(R, F), Q_inv);
183 // m_t = a_t + A_t e_t
184 m = add(a, multiply(A, e));
185 // C = R_t - A_t Q_t A_t'
186 C = subtract(R, quad_form_sym(Q, transpose(A)));
187 lp -= 0.5 * (log_determinant_spd(Q) + trace_quad_form(Q_inv, e));
188 }
189 }
190 return lp;
191 }
192
193 template <typename T_y,
194 typename T_F, typename T_G,
195 typename T_V, typename T_W,
196 typename T_m0, typename T_C0
197 >
198 inline
199 typename return_type<
200 T_y,
201 typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type >::type
202 gaussian_dlm_obs_log(const Eigen::Matrix
203 <T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
204 const Eigen::Matrix
205 <T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
206 const Eigen::Matrix
207 <T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
208 const Eigen::Matrix
209 <T_V, Eigen::Dynamic, Eigen::Dynamic>& V,
210 const Eigen::Matrix
211 <T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
212 const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
213 const Eigen::Matrix
214 <T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
215 return gaussian_dlm_obs_log<false>(y, F, G, V, W, m0, C0);
216 }
217
253 template <bool propto,
254 typename T_y,
255 typename T_F, typename T_G,
256 typename T_V, typename T_W,
257 typename T_m0, typename T_C0
258 >
259 typename return_type<
260 T_y,
261 typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type >::type
262 gaussian_dlm_obs_log(const Eigen::Matrix
263 <T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
264 const Eigen::Matrix
265 <T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
266 const Eigen::Matrix
267 <T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
268 const Eigen::Matrix<T_V, Eigen::Dynamic, 1>& V,
269 const Eigen::Matrix
270 <T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
271 const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
272 const Eigen::Matrix
273 <T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
274 static const char* function("gaussian_dlm_obs_log");
275 typedef
276 typename return_type
277 <T_y, typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type>::type
278 T_lp;
279 T_lp lp(0.0);
280
281 using std::log;
282
283 int r = y.rows(); // number of variables
284 int T = y.cols(); // number of observations
285 int n = G.rows(); // number of states
286
287 check_finite(function, "y", y);
288 check_not_nan(function, "y", y);
289 check_size_match(function,
290 "columns of F", F.cols(),
291 "rows of y", y.rows());
292 check_size_match(function,
293 "rows of F", F.rows(),
294 "rows of G", G.rows());
295 check_finite(function, "F", F);
296 check_not_nan(function, "F", F);
297 check_size_match(function,
298 "rows of G", G.rows(),
299 "columns of G", G.cols());
300 check_finite(function, "G", G);
301 check_not_nan(function, "G", G);
302 check_nonnegative(function, "V", V);
303 check_size_match(function,
304 "size of V", V.size(),
305 "rows of y", y.rows());
306 // TODO(anyone): support infinite V
307 check_finite(function, "V", V);
308 check_not_nan(function, "V", V);
309 check_spsd_matrix(function, "W", W);
310 check_size_match(function,
311 "rows of W", W.rows(),
312 "rows of G", G.rows());
313 // TODO(anyone): support infinite W
314 check_finite(function, "W", W);
315 check_not_nan(function, "W", W);
316 check_size_match(function,
317 "size of m0", m0.size(),
318 "rows of G", G.rows());
319 check_finite(function, "m0", m0);
320 check_not_nan(function, "m0", m0);
321 check_cov_matrix(function, "C0", C0);
322 check_size_match(function,
323 "rows of C0", C0.rows(),
324 "rows of G", G.rows());
325 check_finite(function, "C0", C0);
326 check_not_nan(function, "C0", C0);
327
328 if (y.cols() == 0 || y.rows() == 0)
329 return lp;
330
331 if (include_summand<propto>::value) {
332 lp += 0.5 * NEG_LOG_TWO_PI * r * T;
333 }
334
335 if (include_summand<propto, T_y, T_F, T_G, T_V, T_W, T_m0, T_C0>::value) {
336 T_lp f;
337 T_lp Q;
338 T_lp Q_inv;
339 T_lp e;
340 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> A(n);
341 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> Fj(n);
342 Eigen::Matrix<T_lp, Eigen::Dynamic, 1> m(n);
343 Eigen::Matrix<T_lp, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
344
345 // TODO(anyone): how to recast matrices
346 for (int i = 0; i < m0.size(); i++) {
347 m(i) = m0(i);
348 }
349 for (int i = 0; i < C0.rows(); i++) {
350 for (int j = 0; j < C0.cols(); j++) {
351 C(i, j) = C0(i, j);
352 }
353 }
354
355 for (int i = 0; i < y.cols(); i++) {
356 // Predict state
357 // reuse m and C instead of using a and R
358 m = multiply(G, m);
359 C = add(quad_form_sym(C, transpose(G)), W);
360 for (int j = 0; j < y.rows(); ++j) {
361 // predict observation
362 T_lp yij(y(j, i));
363 // dim Fj = (n, 1)
364 for (int k = 0; k < F.rows(); ++k) {
365 Fj(k) = F(k, j);
366 }
367 // f_{t, i} = F_{t, i}' m_{t, i-1}
368 f = dot_product(Fj, m);
369 Q = trace_quad_form(C, Fj) + V(j);
370 Q_inv = 1.0 / Q;
371 // filtered observation
372 // e_{t, i} = y_{t, i} - f_{t, i}
373 e = yij - f;
374 // A_{t, i} = C_{t, i-1} F_{t, i} Q_{t, i}^{-1}
375 A = multiply(multiply(C, Fj), Q_inv);
376 // m_{t, i} = m_{t, i-1} + A_{t, i} e_{t, i}
377 m += multiply(A, e);
378 // c_{t, i} = C_{t, i-1} - Q_{t, i} A_{t, i} A_{t, i}'
379 // tcrossprod throws an error (ambiguous)
380 // C = subtract(C, multiply(Q, tcrossprod(A)));
381 C -= multiply(Q, multiply(A, transpose(A)));
382 C = 0.5 * add(C, transpose(C));
383 lp -= 0.5 * (log(Q) + pow(e, 2) * Q_inv);
384 }
385 }
386 }
387 return lp;
388 }
389
390 template <typename T_y,
391 typename T_F, typename T_G,
392 typename T_V, typename T_W,
393 typename T_m0, typename T_C0>
394 inline
395 typename return_type
396 <T_y, typename return_type<T_F, T_G, T_V, T_W, T_m0, T_C0>::type>::type
397 gaussian_dlm_obs_log
398 (const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
399 const Eigen::Matrix<T_F, Eigen::Dynamic, Eigen::Dynamic>& F,
400 const Eigen::Matrix<T_G, Eigen::Dynamic, Eigen::Dynamic>& G,
401 const Eigen::Matrix<T_V, Eigen::Dynamic, 1>& V,
402 const Eigen::Matrix<T_W, Eigen::Dynamic, Eigen::Dynamic>& W,
403 const Eigen::Matrix<T_m0, Eigen::Dynamic, 1>& m0,
404 const Eigen::Matrix<T_C0, Eigen::Dynamic, Eigen::Dynamic>& C0) {
405 return gaussian_dlm_obs_log<false>(y, F, G, V, W, m0, C0);
406 }
407
408 }
409 }
410 #endif
check_not_nan.hpp
multiply.hpp
tcrossprod.hpp
log_determinant_spd.hpp
stan::math::trace_quad_form
fvar< T > trace_quad_form(const Eigen::Matrix< fvar< T >, RA, CA > &A, const Eigen::Matrix< fvar< T >, RB, CB > &B)
Definition: trace_quad_form.hpp:18
stan::math::check_not_nan
bool check_not_nan(const char *function, const char *name, const T_y &y)
Return true if y is not NaN.
Definition: check_not_nan.hpp:60
check_size_match.hpp
stan
Definition: log_sum_exp.hpp:8
stan::math::subtract
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
stan::math::gaussian_dlm_obs_log
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).
Definition: gaussian_dlm_obs_log.hpp:78
stan::math::log
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:14
stan::math::multiply
Eigen::Matrix< fvar< T >, R1, C1 > multiply(const Eigen::Matrix< fvar< T >, R1, C1 > &m, const fvar< T > &c)
Definition: multiply.hpp:20
add.hpp
stan::return_type
Metaprogram to calculate the base scalar return type resulting from promoting all the scalar types of...
Definition: return_type.hpp:19
stan::math::include_summand
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
Definition: include_summand.hpp:35
stan::math::check_spsd_matrix
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.
Definition: check_spsd_matrix.hpp:31
log.hpp
check_cov_matrix.hpp
transpose.hpp
stan::return_type::type
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
include_summand.hpp
stan::math::LOG_TWO_PI
const double LOG_TWO_PI
Definition: constants.hpp:190
check_nonnegative.hpp
quad_form_sym.hpp
trace_quad_form.hpp
stan::math::check_size_match
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.
Definition: check_size_match.hpp:30
stan::math::dot_product
fvar< T > dot_product(const Eigen::Matrix< fvar< T >, R1, C1 > &v1, const Eigen::Matrix< fvar< T >, R2, C2 > &v2)
Definition: dot_product.hpp:18
dot_product.hpp
stan::math::check_cov_matrix
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.
Definition: check_cov_matrix.hpp:31
stan::math::e
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:94
stan::math::check_finite
bool check_finite(const char *function, const char *name, const T_y &y)
Return true if y is finite.
Definition: check_finite.hpp:60
stan::math::quad_form_sym
Eigen::Matrix< fvar< T >, CB, CB > quad_form_sym(const Eigen::Matrix< fvar< T >, RA, CA > &A, const Eigen::Matrix< double, RB, CB > &B)
Definition: quad_form_sym.hpp:14
constants.hpp
check_finite.hpp
stan::math::inverse_spd
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:20
stan::math::NEG_LOG_TWO_PI
const double NEG_LOG_TWO_PI
Definition: constants.hpp:192
stan::math::add
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
stan::math::pow
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:17
stan::math::check_nonnegative
bool check_nonnegative(const char *function, const char *name, const T_y &y)
Return true if y is non-negative.
Definition: check_nonnegative.hpp:64
stan::math::transpose
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
quad_form.hpp
subtract.hpp
stan::math::check_square
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.
Definition: check_square.hpp:28
check_spsd_matrix.hpp
inverse_spd.hpp
stan::math::log_determinant_spd
T log_determinant_spd(const Eigen::Matrix< T, R, C > &m)
Returns the log absolute determinant of the specified square matrix.
Definition: log_determinant_spd.hpp:19

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