Tagged: gaussian Toggle Comment Threads | Keyboard Shortcuts

  • wangxinxi 00:47 on March 26, 2011 Permalink | Reply
    Tags: gaussian, linear algebra,   

    I find the previous code for the multivariate gaussian totally crap. https://xinxiwang.wordpress.com/2011/03/23/multivariate-gaussian/

    Since the covariance matrix of gaussian distribution is symmetric positive definite which can be decomposed into LL’ by Cholesky decomposition, where L is a lower triangular matrix. Cholesky decomposition is much faster and stable than the general LU decomposition algorithm.

    Another place can be improved is the determinant. After the LU decomposition is done, we can easily get the determinant by multiply together the diagonal elements of L and U, which can be done together with inversion in one function.

    Advertisements
     
  • wangxinxi 19:34 on March 23, 2011 Permalink | Reply
    Tags: gaussian, matrix   

    Reinvent the wheel.

    #pragma once
    #include <boost/numeric/ublas/vector.hpp>
    #include <boost/numeric/ublas/vector_proxy.hpp>
    #include <boost/numeric/ublas/matrix.hpp>
    #include <boost/numeric/ublas/triangular.hpp>
    #include <boost/numeric/ublas/lu.hpp>
    #include <boost/numeric/ublas/io.hpp>
    #include <math.h>
    
    #define PI 3.14159265358979323846
    
    template<class matrix_T>
    double determinant(const boost::numeric::ublas::matrix_expression<matrix_T> & mat_r)
    {
        double det = 1.0;
    
        matrix_T mLu(mat_r() );
        boost::numeric::ublas::permutation_matrix<std::size_t> pivots(mat_r().size1() );
    
        int is_singular = lu_factorize(mLu, pivots);
    
        if (!is_singular)
        {
    	for (std::size_t i=0; i < pivots.size(); ++i)
    	{
    	    if (pivots(i) != i)
    		det *= -1.0;
    
    	    det *= mLu(i,i);
    	}
        }
        else
    	det = 0.0;
    
        return det;
    }
    
     /* Matrix inversion routine.
        Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
     template<class T>
     bool invertMatrix(const boost::numeric::ublas::matrix<T>& input, 
    		   boost::numeric::ublas::matrix<T>& inverse) {
     	using namespace boost::numeric::ublas;
     	typedef permutation_matrix<std::size_t> pmatrix;
     	// create a working copy of the input
     	matrix<T> A(input);
     	// create a permutation matrix for the LU-factorization
     	pmatrix pm(A.size1());
     	// perform LU-factorization
     	int res = lu_factorize(A,pm);
            if( res != 0 ) return false;
     	// create identity matrix of "inverse"
     	inverse.assign(identity_matrix<T>(A.size1()));
     	// backsubstitute to get the inverse
     	lu_substitute(A, pm, inverse);
     	return true;
     }
    
    template<class T>
    double  gaussian(const boost::numeric::ublas::vector<T> & f,
    		 const boost::numeric::ublas::vector<T> & mu,
    		 const boost::numeric::ublas::matrix<T> & sigma)
    {
        boost::numeric::ublas::matrix<T> sigma_inverse(sigma.size1(), sigma.size2());
        invertMatrix(sigma, sigma_inverse);
    
        boost::numeric::ublas::vector<T> t = f-mu;
    
        double e = inner_prod(prod(trans(t), sigma_inverse), t);
        double d = pow(2*PI, f.size()/2.0)*sqrt(determinant(sigma));
    
        return 1/d*exp(-0.5*e);
    }

     

     
c
Compose new post
j
Next post/Next comment
k
Previous post/Previous comment
r
Reply
e
Edit
o
Show/Hide comments
t
Go to top
l
Go to login
h
Show/Hide help
shift + esc
Cancel