Gaussian Multivariate Distribution -Part 1

pi19404

5.00/5 (1 vote)

Aug 18, 2014

CPOL

2 min read

10336

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:

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:

   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">Σ=LL^T
Σ⁻¹=(LL^T)⁻¹=L^−TL⁻¹
y^TΣ⁻¹y
y^TΣ⁻¹y
y^TL^−TL⁻¹y
(L⁻¹y)^T(L⁻¹y)<mtr><mtd><mrow><mo>
Thus, we can only compute <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow> <mo>(L⁻¹y)<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:

   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:

Gaussian Multivariate Distribution -Part 1 by pi194043