1 /* _______________________________________________________________________
2
3 PECOS: Parallel Environment for Creation Of Stochastics
4 Copyright (c) 2011, Sandia National Laboratories.
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Pecos directory.
7 _______________________________________________________________________ */
8
9 #ifndef MULTIVARIATE_NORMAL_DISTRIBUTION_HPP
10 #define MULTIVARIATE_NORMAL_DISTRIBUTION_HPP
11
12 #include "MultivariateDistribution.hpp"
13 #include "NormalRandomVariable.hpp"
14
15 namespace Pecos {
16
17
18 /// Derived multivariate distribution class for a multivariate normal (MVN)
19
20 /** An MVN can be compactly defined by a mean vector and a symmetric
21 covariance matric. */
22
23 class MultivariateNormalDistribution: public MultivariateDistribution
24 {
25 public:
26
27 //
28 //- Heading: Constructors and destructor
29 //
30
31 MultivariateNormalDistribution(); ///< constructor
32 MultivariateNormalDistribution(const RealVector& means,
33 const RealSymMatrix& cov); ///< alt constructor
34 ~MultivariateNormalDistribution(); ///< destructor
35
36 //
37 //- Heading: Member functions
38 //
39
40 void update(const RealVector& means, const RealSymMatrix& cov);
41
42 void initialize_correlations();
43
44 /// update vector values
45 void push_parameters(short dist_param, const RealVector& values);
46 /// update symmetric matrix values
47 void push_parameters(short dist_param, const RealSymMatrix& values);
48
49 // migrated from NormalRandomVariable static fns:
50 static Real std_pdf(Real beta, size_t n);
51 static Real std_pdf(const RealVector& u);
52
53 protected:
54
55 //
56 //- Heading: Virtual function redefinitions
57 //
58
59 /// return the multivariate PDF value for full set of random variables
60 Real pdf(const RealVector& pt) const;
61 /// return the multivariate log PDF value for full set of random variables
62 Real log_pdf(const RealVector& pt) const;
63
64 //
65 //- Heading: Data
66 //
67
68 /// vector of mean values for multivariate Normal distribution
69 RealVector mvnMeans;
70 /// symmetric covariance matrix for multivariate Normal distribution
71 RealSymMatrix mvnCovariance;
72 };
73
74
MultivariateNormalDistribution()75 inline MultivariateNormalDistribution::MultivariateNormalDistribution():
76 MultivariateDistribution(BaseConstructor())
77 { }
78
79
80 inline MultivariateNormalDistribution::
MultivariateNormalDistribution(const RealVector & means,const RealSymMatrix & cov)81 MultivariateNormalDistribution(const RealVector& means,
82 const RealSymMatrix& cov):
83 MultivariateDistribution(BaseConstructor()),
84 mvnMeans(means), mvnCovariance(cov)
85 { initialize_correlations(); }
86
87
~MultivariateNormalDistribution()88 inline MultivariateNormalDistribution::~MultivariateNormalDistribution()
89 { }
90
91
92 inline void MultivariateNormalDistribution::
update(const RealVector & means,const RealSymMatrix & cov)93 update(const RealVector& means, const RealSymMatrix& cov)
94 { mvnMeans = means; mvnCovariance = cov; initialize_correlations(); }
95
96
initialize_correlations()97 inline void MultivariateNormalDistribution::initialize_correlations()
98 {
99 size_t i, j, num_rv = mvnCovariance.numRows();
100 correlationFlag = false;
101 for (i=1; i<num_rv; i++)
102 for (j=0; j<i; j++)
103 if (std::abs(mvnCovariance(i,j)) > SMALL_NUMBER)
104 { correlationFlag = true; break; }
105 }
106
107
108 inline void MultivariateNormalDistribution::
push_parameters(short dist_param,const RealVector & values)109 push_parameters(short dist_param, const RealVector& values)
110 {
111 switch (dist_param) {
112 case N_MEAN: case N_LOCATION: mvnMeans = values; break;
113 default:
114 PCerr << "Error: lookup failure for distribution parameter " << dist_param
115 << " in MultivariateNormalDistribution::parameters(RealVector)."
116 << std::endl;
117 abort_handler(-1); break;
118 }
119 }
120
121
122 inline void MultivariateNormalDistribution::
push_parameters(short dist_param,const RealSymMatrix & values)123 push_parameters(short dist_param, const RealSymMatrix& values)
124 {
125 bool err_flag = false;
126 switch (dist_param) {
127 case N_VARIANCE: //case N_STD_DEV: case N_SCALE:
128 mvnCovariance = values; initialize_correlations(); break;
129 default:
130 PCerr << "Error: update failure for distribution parameter " << dist_param
131 << " in MultivariateNormalDistribution::parameters(RealSymMatrix)."
132 << std::endl;
133 abort_handler(-1); break;
134 }
135 }
136
137
pdf(const RealVector & x) const138 inline Real MultivariateNormalDistribution::pdf(const RealVector& x) const
139 {
140 // *** TO DO:
141 // std::exp(- (x-mvnMeans)^T inverse_mvnCovar * (x-mvnMeans) / 2.)
142 // / std::sqrt( determinant(mvnCovar) * (2.*Pi)^d )
143 return 0.; //
144 }
145
146
147 // Multivariate standard normal density function from vector.
log_pdf(const RealVector & x) const148 inline Real MultivariateNormalDistribution::log_pdf(const RealVector& x) const
149 {
150 // *** TO DO:
151 // std::log(term) - (x-mvnMeans)^T inverse_mvnCovar * (x-mvnMeans) / 2.
152 return 0.;
153 }
154
155
156 // *** The following are fns migrated from NormalRandomVariable::mvn_*()
157
158 // Multivariate standard normal density function with aggregate distance.
std_pdf(Real beta,size_t n)159 inline Real MultivariateNormalDistribution::std_pdf(Real beta, size_t n)
160 {
161 // need n instances of 1/sqrt(2Pi), but 1D pdf only includes 1:
162 return (n > 1) ? // correct the 1D pdf for n dimensions
163 NormalRandomVariable::std_pdf(beta) * std::pow(2.*PI, -((Real)(n-1))/2.) :
164 NormalRandomVariable::std_pdf(beta);
165 }
166
167
168 // Multivariate standard normal density function from vector.
std_pdf(const RealVector & u)169 inline Real MultivariateNormalDistribution::std_pdf(const RealVector& u)
170 {
171 return std_pdf(u.normFrobenius(), u.length());
172
173 // Alternate implementation invokes exp() repeatedly:
174 //normal_dist norm(0., 1.);
175 //size_t i, n = u.length(); Real pdf = 1.;
176 //for (i=0; i<n; ++i)
177 // pdf *= bmth::pdf(norm, u[i]);
178 }
179
180 } // namespace Pecos
181
182 #endif
183