SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
GaussianDistribution.cpp
Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Written (W) 2013 Heiko Strathmann
00008  */
00009 
00010 #ifdef HAVE_EIGEN3
00011 
00012 #include <shogun/distributions/classical/GaussianDistribution.h>
00013 #include <shogun/base/Parameter.h>
00014 #include <shogun/mathematics/eigen3.h>
00015 
00016 using namespace shogun;
00017 using namespace Eigen;
00018 
00019 CGaussianDistribution::CGaussianDistribution() : CProbabilityDistribution()
00020 {
00021     init();
00022 }
00023 
00024 CGaussianDistribution::CGaussianDistribution(SGVector<float64_t> mean,
00025         SGMatrix<float64_t> cov, bool cov_is_factor) :
00026                 CProbabilityDistribution(mean.vlen)
00027 {
00028     REQUIRE(cov.num_rows==cov.num_cols, "Covariance must be square but is "
00029             "%dx%d\n", cov.num_rows, cov.num_cols);
00030     REQUIRE(mean.vlen==cov.num_cols, "Mean must have same dimension as "
00031             "covariance, which is %dx%d, but is %d\n",
00032             cov.num_rows, cov.num_cols, mean.vlen);
00033 
00034     init();
00035 
00036     m_mean=mean;
00037 
00038     if (!cov_is_factor)
00039     {
00040         Map<MatrixXd> eigen_cov(cov.matrix, cov.num_rows, cov.num_cols);
00041         m_L=SGMatrix<float64_t>(cov.num_rows, cov.num_cols);
00042         Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00043 
00044         /* compute cholesky */
00045         LLT<MatrixXd> llt(eigen_cov);
00046         if (llt.info()==NumericalIssue)
00047         {
00048             /* try to compute smalles eigenvalue for information */
00049             SelfAdjointEigenSolver<MatrixXd> solver(eigen_cov);
00050             if (solver.info() == Success)
00051             {
00052                 VectorXd ev=solver.eigenvalues();
00053                 SG_ERROR("Error computing Cholesky of Gaussian's covariance. "
00054                         "Smallest Eigenvalue is %f.\n", ev[0]);
00055             }
00056         }
00057 
00058         eigen_L=llt.matrixL();
00059     }
00060     else
00061         m_L=cov;
00062 }
00063 
00064 CGaussianDistribution::~CGaussianDistribution()
00065 {
00066 
00067 }
00068 
00069 SGMatrix<float64_t> CGaussianDistribution::sample(int32_t num_samples,
00070         SGMatrix<float64_t> pre_samples) const
00071 {
00072     REQUIRE(num_samples>0, "Number of samples (%d) must be positive\n",
00073             num_samples);
00074 
00075     /* use pre-allocated samples? */
00076     SGMatrix<float64_t> samples;
00077     if (pre_samples.matrix)
00078     {
00079         REQUIRE(pre_samples.num_rows==m_dimension, "Dimension of pre-samples"
00080                 " (%d) does not match dimension of Gaussian (%d)\n",
00081                 pre_samples.num_rows, m_dimension);
00082 
00083         REQUIRE(pre_samples.num_cols==num_samples, "Number of pre-samples"
00084                 " (%d) does not desired number of samples (%d)\n",
00085                 pre_samples.num_cols, num_samples);
00086 
00087         samples=pre_samples;
00088     }
00089     else
00090     {
00091         /* allocate memory and sample from std normal */
00092         samples=SGMatrix<float64_t>(m_dimension, num_samples);
00093         for (index_t i=0; i<m_dimension*num_samples; ++i)
00094             samples.matrix[i]=sg_rand->std_normal_distrib();
00095     }
00096 
00097     /* map into desired Gaussian covariance */
00098     Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows,
00099             samples.num_cols);
00100     Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00101     eigen_samples=eigen_L*eigen_samples;
00102 
00103     /* add mean */
00104     Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen);
00105     eigen_samples.colwise()+=eigen_mean;
00106 
00107     return samples;
00108 }
00109 
00110 SGVector<float64_t> CGaussianDistribution::log_pdf_multiple(SGMatrix<float64_t> samples) const
00111 {
00112     REQUIRE(samples.num_cols>0, "Number of samples must be positive, but is %d\n",
00113             samples.num_cols);
00114     REQUIRE(samples.num_rows==m_dimension, "Sample dimension (%d) does not match"
00115             "Gaussian dimension (%d)\n", samples.num_rows, m_dimension);
00116 
00117     /* for easier to read code */
00118     index_t num_samples=samples.num_cols;
00119 
00120     float64_t const_part=-0.5 * m_dimension * CMath::log(2 * CMath::PI);
00121 
00122     /* determinant is product of diagonal elements of triangular matrix */
00123     float64_t log_det_part=0;
00124     Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
00125     VectorXd diag=eigen_L.diagonal();
00126     log_det_part=-diag.array().log().sum();
00127 
00128     /* sample based part */
00129     Map<MatrixXd> eigen_samples(samples.matrix, samples.num_rows,
00130             samples.num_cols);
00131     Map<VectorXd> eigen_mean(m_mean.vector, m_mean.vlen);
00132 
00133     /* substract mean from samples (create copy) */
00134     SGMatrix<float64_t> centred(m_dimension, num_samples);
00135     Map<MatrixXd> eigen_centred(centred.matrix, centred.num_rows,
00136             centred.num_cols);
00137     for (index_t dim=0; dim<m_dimension; ++dim)
00138     {
00139         for (index_t sample_idx=0; sample_idx<num_samples; ++sample_idx)
00140             centred(dim,sample_idx)=samples(dim,sample_idx)-m_mean[dim];
00141     }
00142 
00143     /* solve the linear system based on factorization */
00144     MatrixXd solved=eigen_L.triangularView<Lower>().solve(eigen_centred);
00145     solved=eigen_L.transpose().triangularView<Upper>().solve(solved);
00146 
00147     /* one quadratic part x^T C^-1 x for each sample x */
00148     SGVector<float64_t> result(num_samples);
00149     Map<VectorXd> eigen_result(result.vector, result.vlen);
00150     for (index_t i=0; i<num_samples; ++i)
00151     {
00152         /* i-th centred sample */
00153         VectorXd left=eigen_centred.block(0, i, m_dimension, 1);
00154 
00155         /* inverted covariance times i-th centred sample */
00156         VectorXd right=solved.block(0,i,m_dimension,1);
00157         result[i]=-0.5*left.dot(right);
00158     }
00159 
00160     /* combine and return */
00161     eigen_result=eigen_result.array()+(log_det_part+const_part);
00162 
00163     /* contains everything */
00164     return result;
00165 
00166 }
00167 
00168 void CGaussianDistribution::init()
00169 {
00170     SG_ADD(&m_mean, "mean", "Mean of the Gaussian.", MS_NOT_AVAILABLE);
00171     SG_ADD(&m_L, "L", "Lower factor of covariance matrix, "
00172             "depending on the factorization type.", MS_NOT_AVAILABLE);
00173 }
00174 
00175 #endif // HAVE_EIGEN3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation