18 using namespace Eigen;
32 "covariance, which is %dx%d, but is %d\n",
46 LLT<MatrixXd> llt(eigen_cov);
47 if (llt.info()==NumericalIssue)
50 SelfAdjointEigenSolver<MatrixXd> solver(eigen_cov);
51 if (solver.info() == Success)
53 VectorXd ev=solver.eigenvalues();
54 SG_ERROR(
"Error computing Cholesky of Gaussian's covariance. "
55 "Smallest Eigenvalue is %f.\n", ev[0]);
59 eigen_L=llt.matrixL();
73 REQUIRE(num_samples>0,
"Number of samples (%d) must be positive\n",
81 " (%d) does not match dimension of Gaussian (%d)\n",
85 " (%d) does not desired number of samples (%d)\n",
94 for (
index_t i=0; i<m_dimension*num_samples; ++i)
102 eigen_samples=eigen_L*eigen_samples;
106 eigen_samples.colwise()+=eigen_mean;
113 REQUIRE(samples.
num_cols>0,
"Number of samples must be positive, but is %d\n",
126 VectorXd diag=eigen_L.diagonal();
127 log_det_part=-diag.array().log().sum();
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];
145 MatrixXd solved=eigen_L.triangularView<Lower>().solve(eigen_centred);
146 solved=eigen_L.transpose().triangularView<Upper>().solve(solved);
151 for (
index_t i=0; i<num_samples; ++i)
154 VectorXd left=eigen_centred.block(0, i, m_dimension, 1);
157 VectorXd right=solved.block(0,i,m_dimension,1);
158 result[i]=-0.5*left.dot(right);
162 eigen_result=eigen_result.array()+(log_det_part+const_part);
169 void CGaussianDistribution::init()
172 SG_ADD(&
m_L,
"L",
"Lower factor of covariance matrix, "
176 #endif // HAVE_EIGEN3
virtual ~CGaussianDistribution()
float64_t std_normal_distrib() const
A base class for representing n-dimensional probability distribution over the real numbers (64bit) fo...
virtual SGVector< float64_t > log_pdf_multiple(SGMatrix< float64_t > samples) const
all of classes and functions are contained in the shogun namespace
SGVector< float64_t > m_mean
static float64_t log(float64_t v)
SGMatrix< float64_t > m_L
virtual SGVector< float64_t > sample() const
static const float64_t PI