IMP  2.0.0
The Integrative Modeling Platform
MultivariateFNormalSufficientSparse.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/MultivariateFNormalSufficientSparse.h
3  * \brief Normal distribution of Function
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_SPARSE_H
9 #define IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_SPARSE_H
10 
11 #include <IMP/isd/isd_config.h>
12 
13 #ifdef IMP_ISD_USE_CHOLMOD
14 
15 #include <IMP/macros.h>
16 #include <IMP/Model.h>
17 #include <IMP/constants.h>
18 #include <math.h>
19 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
20 #include <Eigen/Dense>
21 #include <Eigen/Sparse>
22 #include <unsupported/Eigen/CholmodSupport>
23 #include <ufsparse/cholmod.h>
24 
25 
26 IMPISD_BEGIN_NAMESPACE
27 using Eigen::SparseMatrix;
28 using Eigen::MatrixXd;
29 using Eigen::VectorXd;
30 
31 //! MultivariateFNormalSufficientSparse
32 /** Probability density function and -log(p) of multivariate normal
33  distribution of N M-variate observations.
34 
35  \f[ p(x_1,\cdots,x_N|\mu,F,\Sigma)
36  = \left((2\pi)^M|\Sigma|\right)^{-N/2} J(F)
37  \exp\left(-\frac{1}{2}
38  \sum_{i=1}^N {}^t(F(\mu) - F(x_i))\Sigma^{-1}(F(\mu)-F(x_i))
39  \right)
40  \f]
41  which is implemented as
42  \f[ p(x_1,\cdots,x_N|\mu,F,\Sigma) = ((2\pi)^M|\Sigma|)^{-N/2} J(F)
43  \exp\left(-\frac{N}{2} {}^t\epsilon \Sigma^{-1} \epsilon\right)
44  \exp\left(-\frac{1}{2} \operatorname{tr}(W\Sigma^{-1})\right)
45  \f]
46  where
47  \f[\epsilon = (F(\mu)- \overline{F(x)}) \quad
48  \overline{F(x)} = \frac{1}{N} \sum_{i=1}^N F(x_i)\f]
49  \f( W = \sum_{i=1}^N (F(x_i) - \overline{F(x)}){}^t(F(x_i)
50  - \overline{F(x)}) \f)
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, though
63  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. The inverse of \f$\Sigma\f$ is computed at
67  creation time, along with the sufficient statistics and \f$\epsilon\f$.
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 It uses cholmod and expects a cholmod_common struct as last argument
73  of constructor. cholmod_start should be called beforehand, and
74  cholmod_finish after destruction of all instances of this class.
75 
76  \note References:
77  - Multivariate Likelihood:
78  - Box and Tiao, "Bayesian Inference in Statistical Analysis",
79  Addison-Wesley publishing company, 1973, pp 423.
80  - Factorization in terms of sufficient statistics:
81  - Daniel Fink, "A Compendium of Conjugate Priors", online, May 1997,
82  p.40-41.
83  - Matrix calculations for derivatives:
84  - Petersen and Pedersen, "The Matrix Cookbook", 2008, matrixcookbook.com
85  - Useful reading on missing data (for the case of varying Nobs):
86  - Little and Rubin, "Statistical Analysis with Missing Data",
87  2nd ed, Wiley, 2002, Chapters 6,7 and 11.
88  */
89 class IMPISDEXPORT MultivariateFNormalSufficientSparse : public Object
90 {
91 public:
92  //! Initialize with all observed data
93  /** \param(in) FX matrix of observations with M columns and N rows.
94  \param(in) JF determinant of Jacobian of F with respect to
95  observation matrix X.
96  \param(in) FM mean vector \f$F(\mu)\f$ of size M.
97  \param(in) Sigma MxM variance-covariance matrix \f$\Sigma\f$ in sparse
98  format. Will use symmetry by only considering the upper corner.
99  \param(in) c the control struct for cholmod.
100  \param(in) cutoff when to consider that coefficients of the W matrix are
101  zero.
102  */
103  MultivariateFNormalSufficientSparse(const MatrixXd& FX, double JF,
104  const VectorXd& FM, const SparseMatrix<double>& Sigma,
105  cholmod_common *c, double cutoff=1e-7);
106 
107  //! Initialize with sufficient statistics
108  /** \param(in) Fbar M-dimensional vector of mean observations.
109  \param(in) JF determinant of Jacobian of F with respect to
110  observation matrix X.
111  \param(in) FM M-dimensional true mean vector \f$\mu\f$.
112  \param(in) Nobs number of observations for each variable.
113  \param(in) W MxM matrix of sample variance-covariances, sparse.
114  \param(in) Sigma MxM variance-covariance matrix Sigma, sparse.
115  \param(in) c the control struct for cholmod.
116  */
117  MultivariateFNormalSufficientSparse(const VectorXd& Fbar, double JF,
118  const VectorXd& FM, int Nobs, const SparseMatrix<double>& W,
119  const SparseMatrix<double>& Sigma, cholmod_common *c);
120 
121  //! probability density function
122  double density() const;
123 
124  //! energy (score) functions, aka -log(p)
125  double evaluate() const;
126 
127  //! gradient of the energy wrt the mean F(M)
128  cholmod_dense *evaluate_derivative_FM() const;
129 
130  //! gradient of the energy wrt the variance-covariance matrix Sigma
131  cholmod_sparse *evaluate_derivative_Sigma() const;
132 
133  //! change of parameters
134  void set_FX(const MatrixXd& f, double cutoff=1e-7);
135 
136  void set_JF(double f);
137 
138  void set_FM(const VectorXd& f);
139 
140  void set_Fbar(const VectorXd& f);
141 
142  void set_W(const SparseMatrix<double>& f);
143 
144  void set_Sigma(const SparseMatrix<double>& f);
145 
146  /* remaining stuff */
147  IMP_OBJECT_INLINE(MultivariateFNormalSufficientSparse,
148  out << "MultivariateFNormalSufficientSparse: "
149  << N_ << " observations of "
150  << M_ << " variables " <<std::endl,
151  {
152  cholmod_free_sparse(&W_,c_);
153  cholmod_free_sparse(&P_,c_);
154  cholmod_free_sparse(&Sigma_,c_);
155  cholmod_free_dense(&epsilon_,c_);
156  cholmod_free_factor(&L_,c_);
157  cholmod_finish(c_);
158  });
159 
160  private:
161 
162  /* return trace(W.P), O(M^2) */
163  double trace_WP() const;
164 
165  /* return transpose(epsilon)*P*epsilon, O(M^2) */
166  double mean_dist() const;
167 
168  /* return P*epsilon*transpose(P*epsilon), O(M^2) */
169  cholmod_sparse *compute_PTP() const;
170 
171  /* return P * W * P, O(M^3) */
172  cholmod_sparse *compute_PWP() const;
173 
174  /* compute epsilon, W and Fbar, O(N*M^2) */
175  void compute_sufficient_statistics(double cutoff);
176 
177  /*computes the discrepancy vector*/
178  void compute_epsilon();
179 
180  VectorXd FM_, Fbar_;
181  double JF_,lJF_,norm_,lnorm_;
182  cholmod_sparse *W_, *Sigma_, *P_, *PW_; // matrices
183  cholmod_dense *epsilon_; // vectors
184  cholmod_common *c_;
185  cholmod_factor *L_;
186  MatrixXd FX_;
187 
188  int N_; //number of repetitions
189  int M_; //number of variables
190 };
191 
192 IMPISD_END_NAMESPACE
193 
194 #endif /* IMP_ISD_USE_CHOLMOD */
195 
196 #endif /* IMPISD_MULTIVARIATE_FNORMAL_SUFFICIENT_SPARSE_H */