IMP  2.3.0
The Integrative Modeling Platform
MultivariateFNormalSufficient.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/MultivariateFNormalSufficient.h
3  * \brief Normal distribution of Function
4  *
5  * Copyright 2007-2014 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_H
9 #define IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_H
10 
11 #include <IMP/isd/isd_config.h>
12 #include "internal/timer.h"
13 #include <IMP/macros.h>
14 #include <IMP/kernel/Model.h>
15 #include <IMP/constants.h>
16 #include <IMP/base/Object.h>
17 #include <math.h>
18 #include <IMP/algebra/eigen3/Eigen/Dense>
19 #include <IMP/algebra/eigen3/Eigen/Cholesky>
20 #include <IMP/isd/internal/cg_eigen.h>
21 
22 IMPISD_BEGIN_NAMESPACE
23 
24 // IMP_MVN_TIMER_NFUNCS is the number of functions used by the timer
25 #define IMP_MVN_TIMER_NFUNCS 11
26 
27 //! MultivariateFNormalSufficient
28 /** Probability density function and -log(p) of multivariate normal
29  * distribution of N M-variate observations.
30  *
31  * \f[ p(x_1,\cdots,x_N|\mu,F,\Sigma) =
32  * \left((2\pi\sigma^2)^M|\Sigma|\right)^{-N/2} J(F)
33  * \exp\left(-\frac{1}{2\sigma^2}
34  * \sum_{i=1}^N {}^t(F(\mu) - F(x_i))\Sigma^{-1}(F(\mu)-F(x_i))
35  * \right)
36  * \f]
37  * which is implemented as
38  * \f[ p(x_1,\cdots,x_N|\mu,F,\Sigma) = ((2\pi\sigma^2)^M|\Sigma|)^{-N/2} J(F)
39  * \exp\left(-\frac{N}{2\sigma^2} {}^t\epsilon \Sigma^{-1} \epsilon\right)
40  * \exp\left(-\frac{1}{2\sigma^2} \text{tr}(W\Sigma^{-1})\right)
41  * \f]
42  * where
43  * \f[\epsilon = (F(\mu)- \overline{F(x)}) \quad
44  * \overline{F(x)} = \frac{1}{N} \sum_{i=1}^N F(x_i)\f]
45  * and
46  * \f[W=\sum_{i=1}^N(F(x_i)-\overline{F(x)}){}^t(F(x_i)-\overline{F(x)}) \f]
47  *
48  * \f$\sigma\f$ is a multiplicative scaling factor that factors out of the
49  * \f$\Sigma\f$ covariance matrix. It is set to 1 by default and its intent is
50  * to avoid inverting the \f$\Sigma\f$ matrix unless necessary.
51  *
52  * Set J(F) to 1 if you want the multivariate normal distribution.
53  * The distribution is normalized with respect to the matrix variable X.
54  * The Sufficient statistics are calculated at initialization.
55  *
56  * Example: if F is the log function, the multivariate F-normal distribution
57  * is the multivariate lognormal distribution with mean \f$\mu\f$ and
58  * standard deviation \f$\Sigma\f$.
59  *
60  * \note This is an implementation of the matrix normal distribution for F(X),
61  * where rows of F(X) are independent and homoscedastic (they represent
62  * repetitions of the same experiment), but columns might be correlated,
63  * though the provided matrix.
64  *
65  * \note For now, F must be monotonically increasing, so that J(F) > 0. The
66  * program will not check for that. Uses a Cholesky (\f$LDL^T\f$)
67  * decomposition of \f$\Sigma\f$, which is recomputed when needed.
68  *
69  * \note All observations must be given, so if you have missing data you might
70  * want to do some imputation on it first.
71  *
72  * \note References:
73  * - Multivariate Likelihood:
74  * Box and Tiao, "Bayesian Inference in Statistical Analysis",
75  * Addison-Wesley publishing company, 1973, pp 423.
76  * - Factorization in terms of sufficient statistics:
77  * Daniel Fink, "A Compendium of Conjugate Priors", online, May 1997, p.40-41.
78  * - Matrix calculations for derivatives:
79  * Petersen and Pedersen, "The Matrix Cookbook", 2008, matrixcookbook.com
80  * - Useful reading on missing data (for the case of varying Nobs):
81  * Little and Rubin, "Statistical Analysis with Missing Data", 2nd ed, Wiley,
82  * 2002, Chapters 6,7 and 11.
83  *
84  */
85 
86 class IMPISDEXPORT MultivariateFNormalSufficient : public base::Object {
87 
88  private:
89  IMP_Eigen::VectorXd FM_, Fbar_, epsilon_, Peps_;
90  double JF_, lJF_, norm_, lnorm_;
91  IMP_Eigen::MatrixXd P_, W_, Sigma_, FX_, PW_;
92  int N_; // number of repetitions
93  int M_; // number of variables
94  // IMP_Eigen::LLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> ldlt_;
95  IMP_Eigen::LDLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> ldlt_;
96  // flags are true if the corresponding object is up to date.
97  bool flag_FM_, flag_FX_, flag_Fbar_, flag_W_, flag_Sigma_, flag_epsilon_,
98  flag_PW_, flag_P_, flag_ldlt_, flag_norms_, flag_Peps_;
99  // cg-related variables
100  IMP_Eigen::MatrixXd precond_;
101  bool use_cg_, first_PW_, first_PWP_;
102  double cg_tol_;
103  double factor_;
105 
106  internal::CallTimer<IMP_MVN_TIMER_NFUNCS> timer_;
107 
108  public:
109  /** Initialize with all observed data
110 * \param [in] FX F(X) matrix of observations with M columns and N rows.
111 * \param [in] JF J(F) determinant of Jacobian of F with respect to
112 * observation matrix X.
113 * \param [in] FM F(M) mean vector \f$F(\mu)\f$ of size M.
114 * \param [in] Sigma : MxM variance-covariance matrix \f$\Sigma\f$.
115 * \param [in] factor : multiplicative factor (default 1)
116 * */
117  MultivariateFNormalSufficient(const IMP_Eigen::MatrixXd& FX, double JF,
118  const IMP_Eigen::VectorXd& FM,
119  const IMP_Eigen::MatrixXd& Sigma,
120  double factor = 1);
121 
122  /** Initialize with sufficient statistics
123 * \param [in] Fbar : M-dimensional vector of mean observations.
124 * \param [in] JF J(F) determinant of Jacobian of F with respect to observation
125 * matrix X.
126 * \param [in] FM F(M) : M-dimensional true mean vector \f$\mu\f$.
127 * \param [in] Nobs : number of observations for each variable.
128 * \param [in] W : MxM matrix of sample variance-covariances.
129 * \param [in] Sigma : MxM variance-covariance matrix Sigma.
130 * \param [in] factor : multiplicative factor (default 1)
131 * */
132  MultivariateFNormalSufficient(const IMP_Eigen::VectorXd& Fbar, double JF,
133  const IMP_Eigen::VectorXd& FM, int Nobs,
134  const IMP_Eigen::MatrixXd& W,
135  const IMP_Eigen::MatrixXd& Sigma,
136  double factor = 1);
137 
138  /* probability density function */
139  double density() const;
140 
141  /* energy (score) functions, aka -log(p) */
142  double evaluate() const;
143 
144  /* gradient of the energy wrt the mean F(M) */
145  IMP_Eigen::VectorXd evaluate_derivative_FM() const;
146 
147  /* gradient of the energy wrt the variance-covariance matrix Sigma */
148  IMP_Eigen::MatrixXd evaluate_derivative_Sigma() const;
149 
150  // derivative wrt scalar factor
151  double evaluate_derivative_factor() const;
152 
153  /* second derivative wrt FM and FM */
154  IMP_Eigen::MatrixXd evaluate_second_derivative_FM_FM() const;
155 
156  /* second derivative wrt FM(l) and Sigma
157  * row and column indices in the matrix returned are for Sigma
158  */
159  IMP_Eigen::MatrixXd evaluate_second_derivative_FM_Sigma(unsigned l) const;
160 
161  /* second derivative wrt Sigma and Sigma(k,l) */
162  IMP_Eigen::MatrixXd evaluate_second_derivative_Sigma_Sigma(unsigned k,
163  unsigned l) const;
164 
165  /* change of parameters */
166  void set_FX(const IMP_Eigen::MatrixXd& f);
167  IMP_Eigen::MatrixXd get_FX() const;
168 
169  void set_Fbar(const IMP_Eigen::VectorXd& f);
170  IMP_Eigen::VectorXd get_Fbar() const;
171 
172  IMP_Eigen::VectorXd get_epsilon() const;
173 
174  void set_jacobian(double f);
175  double get_jacobian() const;
176  void set_minus_log_jacobian(double f); //-log(J)
177  double get_minus_log_jacobian() const;
178 
179  void set_FM(const IMP_Eigen::VectorXd& f);
180  IMP_Eigen::VectorXd get_FM() const;
181 
182  void set_W(const IMP_Eigen::MatrixXd& f);
183  IMP_Eigen::MatrixXd get_W() const;
184 
185  void set_Sigma(const IMP_Eigen::MatrixXd& f);
186  IMP_Eigen::MatrixXd get_Sigma() const;
187 
188  void set_factor(double f);
189  double get_factor() const;
190 
191  // if you want to force a recomputation of all stored variables
192  void reset_flags();
193 
194  // use conjugate gradients (default false)
195  void set_use_cg(bool use, double tol);
196 
197  // print runtime statistics
198  void stats() const;
199 
200  // return Sigma's eigenvalues from smallest to biggest
201  IMP_Eigen::VectorXd get_Sigma_eigenvalues() const;
202 
203  // return Sigma's condition number
204  double get_Sigma_condition_number() const;
205 
206  // Solve for Sigma*X = B, yielding X
207  IMP_Eigen::MatrixXd solve(IMP_Eigen::MatrixXd B) const;
208 
209  /* return transpose(epsilon)*P*epsilon */
210  double get_mean_square_residuals() const;
211 
212  /* return minus exponent
213  * \f[-\frac{1}{2\sigma^2}
214  * \sum_{i=1}^N {}^t(F(\mu) - F(x_i))\Sigma^{-1}(F(\mu)-F(x_i)) \f]
215  */
216  double get_minus_exponent() const;
217 
218  /* return minus log normalization
219  * \f[\frac{N}{2}\left(\log(2\pi\sigma^2) + \log |\Sigma|\right)
220  * -\log J(F) \f]
221  */
222  double get_minus_log_normalization() const;
223 
224  /* return log |\sigma^2 \Sigma| */
225  double get_log_generalized_variance() const;
226 
227  /* remaining stuff */
229  /*IMP_OBJECT_INLINE(MultivariateFNormalSufficient,
230  out << "MultivariateFNormalSufficient: "
231  << N_ << " observations of "
232  << M_ << " variables " <<std::endl,
233  {});*/
234 
235  private:
236  // conjugate gradient init
237  void setup_cg();
238 
239  // precision matrix
240  IMP_Eigen::MatrixXd get_P() const;
241  void set_P(const IMP_Eigen::MatrixXd& P);
242 
243  // precision * W
244  IMP_Eigen::MatrixXd get_PW() const;
245  IMP_Eigen::MatrixXd compute_PW_direct() const;
246  IMP_Eigen::MatrixXd compute_PW_cg() const;
247  void set_PW(const IMP_Eigen::MatrixXd& PW);
248 
249  // precision * epsilon
250  IMP_Eigen::VectorXd get_Peps() const;
251  void set_Peps(const IMP_Eigen::VectorXd& Peps);
252 
253  // epsilon = Fbar - FM
254  void set_epsilon(const IMP_Eigen::VectorXd& eps);
255 
256  // gets factorization object
257  // IMP_Eigen::LLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> get_ldlt() const;
258  IMP_Eigen::LDLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper> get_ldlt() const;
259  void set_ldlt(
260  const IMP_Eigen::LDLT<IMP_Eigen::MatrixXd, IMP_Eigen::Upper>& ldlt);
261 
262  // compute determinant and norm
263  void set_norms(double norm, double lnorm);
264  std::vector<double> get_norms() const;
265 
266  /* return trace(W.P) */
267  double trace_WP() const;
268 
269  /* return P*epsilon*transpose(P*epsilon) */
270  IMP_Eigen::MatrixXd compute_PTP() const;
271 
272  /* return P * W * P, O(M^2) */
273  IMP_Eigen::MatrixXd compute_PWP() const;
274 
275  /*computes the discrepancy vector*/
276  void compute_epsilon();
277 };
278 
279 IMPISD_END_NAMESPACE
280 
281 #endif /* IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_H */
Declare an efficient stl-compatible map.
Import IMP/kernel/constants.h in the namespace.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Import IMP/kernel/macros.h in the namespace.
Storage of a model, its restraints, constraints and particles.
Common base class for heavy weight IMP objects.
Definition: Object.h:106
A shared base class to help in debugging and things.