IMP logo
IMP Reference Guide  2.17.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-2022 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/Model.h>
15 #include <IMP/constants.h>
16 #include <IMP/Object.h>
17 #include <math.h>
18 #include <Eigen/Dense>
19 #include <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 Object {
87 
88  private:
89  Eigen::VectorXd FM_, Fbar_, epsilon_, Peps_;
90  double JF_, lJF_, norm_, lnorm_;
91  Eigen::MatrixXd P_, W_, Sigma_, FX_, PW_;
92  int N_; // number of repetitions
93  int M_; // number of variables
94  // Eigen::LLT<Eigen::MatrixXd, Eigen::Upper> ldlt_;
95  Eigen::LDLT<Eigen::MatrixXd, 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  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 Eigen::MatrixXd& FX, double JF,
118  const Eigen::VectorXd& FM,
119  const 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
125  observation 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 Eigen::VectorXd& Fbar, double JF,
133  const Eigen::VectorXd& FM, int Nobs,
134  const Eigen::MatrixXd& W,
135  const Eigen::MatrixXd& Sigma,
136  double factor = 1);
137 
138  //! Probability density function
139  double density() const;
140 
141  //! Energy (score) function, aka -log(p)
142  double evaluate() const;
143 
144  //! gradient of the energy wrt the mean F(M)
145  Eigen::VectorXd evaluate_derivative_FM() const;
146 
147  //! gradient of the energy wrt the variance-covariance matrix Sigma
148  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  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  Eigen::MatrixXd evaluate_second_derivative_FM_Sigma(unsigned l) const;
160 
161  //! second derivative wrt Sigma and Sigma(k,l)
162  Eigen::MatrixXd evaluate_second_derivative_Sigma_Sigma(unsigned k,
163  unsigned l) const;
164 
165  /* change of parameters */
166  void set_FX(const Eigen::MatrixXd& f);
167  Eigen::MatrixXd get_FX() const;
168 
169  void set_Fbar(const Eigen::VectorXd& f);
170  Eigen::VectorXd get_Fbar() const;
171 
172  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 Eigen::VectorXd& f);
180  Eigen::VectorXd get_FM() const;
181 
182  void set_W(const Eigen::MatrixXd& f);
183  Eigen::MatrixXd get_W() const;
184 
185  void set_Sigma(const Eigen::MatrixXd& f);
186  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  //! return Sigma's eigenvalues from smallest to biggest
198  Eigen::VectorXd get_Sigma_eigenvalues() const;
199 
200  //! return Sigma's condition number
201  double get_Sigma_condition_number() const;
202 
203  //! Solve for Sigma*X = B, yielding X
204  Eigen::MatrixXd solve(Eigen::MatrixXd B) const;
205 
206  //! return transpose(epsilon)*P*epsilon
207  double get_mean_square_residuals() const;
208 
209  //! return minus exponent
210  /** \f[-\frac{1}{2\sigma^2}
211  \sum_{i=1}^N {}^t(F(\mu) - F(x_i))\Sigma^{-1}(F(\mu)-F(x_i)) \f]
212  */
213  double get_minus_exponent() const;
214 
215  //! return minus log normalization
216  /** \f[\frac{N}{2}\left(\log(2\pi\sigma^2) + \log |\Sigma|\right)
217  -\log J(F) \f]
218  */
219  double get_minus_log_normalization() const;
220 
221  //! return log |\sigma^2 \Sigma|
222  double get_log_generalized_variance() const;
223 
224  /* remaining stuff */
226 
227  private:
228  // conjugate gradient init
229  void setup_cg();
230 
231  // precision matrix
232  Eigen::MatrixXd get_P() const;
233  void set_P(const Eigen::MatrixXd& P);
234 
235  // precision * W
236  Eigen::MatrixXd get_PW() const;
237  Eigen::MatrixXd compute_PW_direct() const;
238  Eigen::MatrixXd compute_PW_cg() const;
239  void set_PW(const Eigen::MatrixXd& PW);
240 
241  // precision * epsilon
242  Eigen::VectorXd get_Peps() const;
243  void set_Peps(const Eigen::VectorXd& Peps);
244 
245  // epsilon = Fbar - FM
246  void set_epsilon(const Eigen::VectorXd& eps);
247 
248  // gets factorization object
249  // Eigen::LLT<Eigen::MatrixXd, Eigen::Upper> get_ldlt() const;
250  Eigen::LDLT<Eigen::MatrixXd, Eigen::Upper> get_ldlt() const;
251  void set_ldlt(
252  const Eigen::LDLT<Eigen::MatrixXd, Eigen::Upper>& ldlt);
253 
254  // compute determinant and norm
255  void set_norms(double norm, double lnorm);
256  std::vector<double> get_norms() const;
257 
258  /* return trace(W.P) */
259  double trace_WP() const;
260 
261  /* return P*epsilon*transpose(P*epsilon) */
262  Eigen::MatrixXd compute_PTP() const;
263 
264  /* return P * W * P, O(M^2) */
265  Eigen::MatrixXd compute_PWP() const;
266 
267  /*computes the discrepancy vector*/
268  void compute_epsilon();
269 };
270 
271 IMPISD_END_NAMESPACE
272 
273 #endif /* IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_H */
Helper functions to check for NaN or infinity.
Various useful constants.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Storage of a model, its restraints, constraints and particles.
Various general useful macros for IMP.
Common base class for heavy weight IMP objects.
Definition: Object.h:106
A shared base class to help in debugging and things.