SHOGUN  4.1.0
GaussianDistribution.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Heiko Strathmann
8  */
9
11
12 #ifdef HAVE_EIGEN3
13
14 #include <shogun/base/Parameter.h>
16
17 using namespace shogun;
18 using namespace Eigen;
19
21 {
22  init();
23 }
24
26  SGMatrix<float64_t> cov, bool cov_is_factor) :
27  CProbabilityDistribution(mean.vlen)
28 {
29  REQUIRE(cov.num_rows==cov.num_cols, "Covariance must be square but is "
30  "%dx%d\n", cov.num_rows, cov.num_cols);
31  REQUIRE(mean.vlen==cov.num_cols, "Mean must have same dimension as "
32  "covariance, which is %dx%d, but is %d\n",
33  cov.num_rows, cov.num_cols, mean.vlen);
34
35  init();
36
37  m_mean=mean;
38
39  if (!cov_is_factor)
40  {
41  Map<MatrixXd> eigen_cov(cov.matrix, cov.num_rows, cov.num_cols);
44
45  /* compute cholesky */
46  LLT<MatrixXd> llt(eigen_cov);
47  if (llt.info()==NumericalIssue)
48  {
49  /* try to compute smalles eigenvalue for information */
51  if (solver.info() == Success)
52  {
53  VectorXd ev=solver.eigenvalues();
54  SG_ERROR("Error computing Cholesky of Gaussian's covariance. "
55  "Smallest Eigenvalue is %f.\n", ev[0]);
56  }
57  }
58
59  eigen_L=llt.matrixL();
60  }
61  else
62  m_L=cov;
63 }
64
66 {
67
68 }
69
71  SGMatrix<float64_t> pre_samples) const
72 {
73  REQUIRE(num_samples>0, "Number of samples (%d) must be positive\n",
74  num_samples);
75
76  /* use pre-allocated samples? */
77  SGMatrix<float64_t> samples;
78  if (pre_samples.matrix)
79  {
80  REQUIRE(pre_samples.num_rows==m_dimension, "Dimension of pre-samples"
81  " (%d) does not match dimension of Gaussian (%d)\n",
82  pre_samples.num_rows, m_dimension);
83
84  REQUIRE(pre_samples.num_cols==num_samples, "Number of pre-samples"
85  " (%d) does not desired number of samples (%d)\n",
86  pre_samples.num_cols, num_samples);
87
88  samples=pre_samples;
89  }
90  else
91  {
92  /* allocate memory and sample from std normal */
93  samples=SGMatrix<float64_t>(m_dimension, num_samples);
94  for (index_t i=0; i<m_dimension*num_samples; ++i)
95  samples.matrix[i]=sg_rand->std_normal_distrib();
96  }
97
98  /* map into desired Gaussian covariance */
99  Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows,
100  samples.num_cols);
102  eigen_samples=eigen_L*eigen_samples;
103
105  Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen);
106  eigen_samples.colwise()+=eigen_mean;
107
108  return samples;
109 }
110
112 {
113  REQUIRE(samples.num_cols>0, "Number of samples must be positive, but is %d\n",
114  samples.num_cols);
115  REQUIRE(samples.num_rows==m_dimension, "Sample dimension (%d) does not match"
116  "Gaussian dimension (%d)\n", samples.num_rows, m_dimension);
117
118  /* for easier to read code */
119  index_t num_samples=samples.num_cols;
120
121  float64_t const_part=-0.5 * m_dimension * CMath::log(2 * CMath::PI);
122
123  /* determinant is product of diagonal elements of triangular matrix */
124  float64_t log_det_part=0;
126  VectorXd diag=eigen_L.diagonal();
127  log_det_part=-diag.array().log().sum();
128
129  /* sample based part */
130  Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows,
131  samples.num_cols);
132  Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen);
133
134  /* substract mean from samples (create copy) */
135  SGMatrix<float64_t> centred(m_dimension, num_samples);
136  Map<MatrixXd> eigen_centred(centred.matrix, centred.num_rows,
137  centred.num_cols);
138  for (index_t dim=0; dim<m_dimension; ++dim)
139  {
140  for (index_t sample_idx=0; sample_idx<num_samples; ++sample_idx)
141  centred(dim,sample_idx)=samples(dim,sample_idx)-m_mean[dim];
142  }
143
144  /* solve the linear system based on factorization */
145  MatrixXd solved=eigen_L.triangularView<Lower>().solve(eigen_centred);
146  solved=eigen_L.transpose().triangularView<Upper>().solve(solved);
147
148  /* one quadratic part x^T C^-1 x for each sample x */
149  SGVector<float64_t> result(num_samples);
150  Map<VectorXd> eigen_result(result.vector, result.vlen);
151  for (index_t i=0; i<num_samples; ++i)
152  {
153  /* i-th centred sample */
154  VectorXd left=eigen_centred.block(0, i, m_dimension, 1);
155
156  /* inverted covariance times i-th centred sample */
157  VectorXd right=solved.block(0,i,m_dimension,1);
158  result[i]=-0.5*left.dot(right);
159  }
160
161  /* combine and return */
162  eigen_result=eigen_result.array()+(log_det_part+const_part);
163
164  /* contains everything */
165  return result;
166
167 }
168
169 void CGaussianDistribution::init()
170 {
171  SG_ADD(&m_mean, "mean", "Mean of the Gaussian.", MS_NOT_AVAILABLE);
172  SG_ADD(&m_L, "L", "Lower factor of covariance matrix, "
173  "depending on the factorization type.", MS_NOT_AVAILABLE);
174 }
175
176 #endif // HAVE_EIGEN3
float64_t std_normal_distrib() const
Definition: Random.cpp:238
int32_t index_t
Definition: common.h:62
Definition: SGMatrix.h:20
#define SG_ERROR(...)
Definition: SGIO.h:129
#define REQUIRE(x,...)
Definition: SGIO.h:206
index_t num_cols
Definition: SGMatrix.h:378
CRandom * sg_rand
Definition: init.cpp:39
index_t num_rows
Definition: SGMatrix.h:376
index_t vlen
Definition: SGVector.h:494
A base class for representing n-dimensional probability distribution over the real numbers (64bit) fo...
double float64_t
Definition: common.h:50
virtual SGVector< float64_t > log_pdf_multiple(SGMatrix< float64_t > samples) const
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static float64_t log(float64_t v)
Definition: Math.h:922