Click here to Skip to main content
65,938 articles
CodeProject is changing. Read more.
Articles / programming / algorithm

Gaussian Multivariate Distribution -Part 1

5.00/5 (1 vote)
17 Aug 2014CPOL2 min read 10.3K  
Multi Variate Gaussian Distribution - Part 1

Introduction

In this article, we will look at the multivariate Gaussian distribution.

  • The mutivariate Gaussian distribution is parameterized by mean vector µ and covariance matrix Σ
    f x ( x 1 , … , x k ) = 1 √ ( 2 π ) k | Σ | exp ( − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) ) ,
    w h e r e
    μ is a Nx1 mean vector
    Σ is NxN matrix covariance matrix
    | Σ | is the determinant of Σ
  • The mutivariate Gaussian is also used as probability density function of the vector <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow class="MJX-TeXAtom-ORD"> <mi>X
  • A naive way of implementation is to directly perform the computation as in the equation to get the result.
  • A Gaussian is defined
  • Since ill be using a lot of opencv conversion routines are defined for Opencv cv::Mat data structure to Eigen::MatrixXf data structure.
  • When loading the covarince matrix ,the determinant, inverse, etc. are also computed which would be later required for probability calculation:
    C#
    void setMean(Mat &v)
    {
        Map<Eigen::Matrix<float,Dynamic,Dynamic,Eigen::RowMajor> > 
        mappedMat((float *)v.data(),1,v.size()) ;
        _mu=mappedMat;
    }
    
    void setSigma(cv::Mat &v)
    {
        Map<Eigen::Matrix<float,Dynamic,Dynamic,Eigen::RowMajor> > 
        mappedMat((float *)v.data,v.rows,v.cols) ;
        _Sigma=mappedMat;
        _dim=v.rows;
        _det=_Sigma.determinant();
        _scale=1.f/(pow(2*PI*_dim*_det,0.5));
        _invsigma=_Sigma.inverse();
    }
    
    MatrixXf setData(cv::Mat &v)
    {
        Map<Eigen::Matrix<float,Dynamic,Dynamic,Eigen::RowMajor> > 
        mappedMat((float *)v.data,1,v.cols) ;
        MatrixXf  ref=mappedMat;
        return ref;
    }
  • The probability computation is simply writing out the formula:
    C#
       void validate(float &res)
       {
    if(res<0)
               res=0;
           if(res>1)
               res=1;
           if(isnan(res))
               res=0;
        }
        
       float Prob(Mat &x)
       {
           MatrixXf tmp=setData(x);
           float res=0;
           MatrixXf tmp1=(tmp-_mu);
           MatrixXf tmp2=tmp1*_invsigma*tmp1.transpose();
           res=scale*tmp2(0,0);        
    validate(res);   
           return res;
       }
  • The covariance matrix is positive semidefinate and symmetric
  • The Cholesky decomposition of a symmetric, positive definite matrix <math xmlns="http://www.w3.org/1998/Math/MathML"> <mi mathvariant="bold">Σ decomposes A into a product of a lower triangular matrix L and its transpose.
  • Thus the product <math display="block" xmlns="http://www.w3.org/1998/Math/MathML"> <mtable columnalign="right center left" columnspacing="0 thickmathspace" displaystyle="true" rowspacing="3pt"> <mtr> <mtd> <mtd> <mi mathvariant="bold">Σ=LLT
    Σ−1=(LLT)−1=L−TL−1
    yTΣ−1y
    yTΣ−1y
    yTL−TL−1y
    (L−1y)T(L−1y)<mtr><mtd><mrow><mo>
  • Thus, we can only compute <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow> <mo>(L−1y)<mo> and take its transpose
  • Thus, we can perform the Cholskey factorization of the covariance matrix to find <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow class="MJX-TeXAtom-ORD"> <mi mathvariant="bold">L
  • Then, we take the inverse of <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow class="MJX-TeXAtom-ORD"> <mi mathvariant="bold">L:
    C#
    float Prob(Mat &x)
    {
        MatrixXf tmp=setData(x);
    
        float res=0;
        MatrixXf tmp1=(tmp-_mu).transpose();
        tmp1=_LI*tmp1;
        MatrixXf tmp2=tmp1.transpose()*tmp1;
        res=tmp2(0,0);
        validate(res);
        return res;
    }
    

Source Code

The code can be found at git repository https://github.com/pi19404/OpenVisionin files ImgML/gaussian.hpp and ImgML.gaussian.cpp files.

The PDF Version of the document can be found below:

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)