15,032,755 members
Articles / General Programming / Algorithms
Article
Posted 17 Aug 2014

8K views
5 bookmarked

# Gaussian Multivariate Distribution -Part 1

Rate me:
17 Aug 2014CPOL2 min read
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: