Compute the the maximum likelihood estimates of a Gaussian statistical model. More...

#include <STK_GaussianStatModel.h>

Inheritance diagram for STK::GaussianStatModel:

Public Member Functions
	GaussianStatModel (Matrix const *p_data)
	constructor.
virtual	~GaussianStatModel ()
	destructor.
virtual bool	run ()
	implementation of the Gaussian statistical model
virtual bool	run (Vector const &weights)
	implementation of the weighted Gaussian statistical model
Point const &	mean () const
	get the empirical mean.
MatrixSquare const &	covariance () const
	get the empirical covariance
Protected Member Functions
virtual void	compMean ()
	compute the empirical mean
virtual void	compCovariance ()
	compute the empirical covariance matrix.
virtual void	compWeightedMean (Vector const &weights)
	compute the empirical weighted mean
virtual void	compWeightedCovariance (Vector const &weights)
	compute the empirical weighted covariance matrix.
Protected Attributes
Point	mean_
	Vector of the empirical means.
MatrixSquare	cov_
	Matrix of the empirical covaiance.
Private Types
typedef ITStatModel< Real, Point, Vector, Matrix >	RealStatModel

Detailed Description

Compute the the maximum likelihood estimates of a Gaussian statistical model.

A random vector $X \in; \mathbb{R}^p$ is Gaussian if

$X\ \sim\ \mathcal{N}(\mu,\ \Sigma).$

The likelihood function of a gaussian sample is:

$L(\mu,\Sigma)=[\mbox{constant}] \prod_{i=1}^n \det(\Sigma)^{-1/2} \exp\left(-\frac{1}{2} (X_i-\mu)^T \Sigma^{-1} (X_i-\mu)\right)$

The maximum likelihood estimator can be performed via matrix calculus formulae Re-write the likelihood in the log form using the trace trick:

$\ln L(\mu,\Sigma) = \operatorname{const} -{n \over 2} \ln \det(\Sigma) -\frac{1}{2} \operatorname{tr} \left[ \Sigma^{-1} \sum_{i=1}^n (X_i-\mu) (X_i-\mu)^T \right].$

The differential of this log-likelihood is

$d \ln L(\mu,\Sigma) = -{n \over 2} \operatorname{tr} \left[ \Sigma^{-1} \left\{ d \Sigma \right\} \right] -\frac{1}{2} \operatorname{tr} \left[ - \Sigma^{-1} \{ d \Sigma \} \Sigma^{-1} \sum_{i=1}^n (X_i-\mu)(X_i-\mu)^T - 2 \Sigma^{-1} \sum_{i=1}^n (X_i - \mu) \{ d \mu \}^T \right].$

It naturally breaks down into the part related to the estimation of the mean, and to the part related to the estimation of the variance. The first order condition for maximum, $d \ln L(\mu,\Sigma)=0$ , is satisfied when the terms multiplying $d \mu$ and $d \Sigma$ are identically zero. Assuming (the maximum likelihood estimate of) $\Sigma$ is non-singular, the first order condition for the estimate of the mean vector is

$\sum_{i=1}^n (X_i - \mu) = 0,$

which leads to the maximum likelihood estimator

$\hat{\mu} = \bar{X} = \frac{1}{n} \sum_{i=1}^n X_i.$

This lets us simplify $\sum_{i=1}^n (X_i-\mu)(X_i-\mu)^T = \sum_{i=1}^n (X_i-\bar{X})(X_i-\bar{X})^T = S$ . Then the terms involving $d \Sigma$ in $d \ln L$ can be combined as

$-\frac{1}{2} \operatorname{tr} \left( \Sigma^{-1} \left\{ d \Sigma \right\} \left[ nI_p - \Sigma^{-1} S \right] \right).$

The first order condition $d \ln L(\mu,\Sigma)=0$ will hold when the term in the square bracket is (matrix-valued) zero. Pre-multiplying the latter by $\Sigma$ and dividing by $n$ gives $\hat{\Sigma} = \frac{1}{n} S,$

Definition at line 103 of file STK_GaussianStatModel.h.

Member Typedef Documentation

typedef ITStatModel<Real, Point, Vector, Matrix> STK::GaussianStatModel::RealStatModel [private]

Definition at line 105 of file STK_GaussianStatModel.h.

Constructor & Destructor Documentation

STK::GaussianStatModel::GaussianStatModel ( Matrix const * p_data )

constructor.

Parameters:

p_data pointer on the data set

Definition at line 47 of file STK_GaussianStatModel.cpp.

References STK::IStatModelBase::nbFreeParameter_, and STK::IStatModelBase::nbVar_.

                                    : RealStatModel(p_data)
                                    , mean_(p_data_->rangeHo())
                                    , cov_(p_data_->rangeHo())
{
  nbFreeParameter_ = nbVar_ + (nbVar_* (nbVar_-1))/2;
}

STK::GaussianStatModel::~GaussianStatModel ( ) [virtual]

destructor.

Definition at line 56 of file STK_GaussianStatModel.cpp.

{}

Member Function Documentation

bool STK::GaussianStatModel::run ( ) [virtual]

implementation of the Gaussian statistical model

Returns:: true if no error occur and false otherwise.

Implements STK::IRunnerBase.

Definition at line 61 of file STK_GaussianStatModel.cpp.

References compCovariance(), compMean(), cov_, STK::IStatModelBase::logLikelihood_, mean_, STK::IRunnerPtr2D< Real, Point, Vector, Matrix >::p_data_, STK::ITStatModel< Real, Point, Vector, Matrix >::p_law_, and STK::Law::ILawBase::update().

Referenced by STK::Gaussian2BlocksStatModel::run().

{
  // compute the mean
  compMean();
  // compute the covariance matrix
  compCovariance();
  // create p_law_ (will be deleted in base class)
  // update gaussian law (will be deleted in base class)
  if (!p_law_) p_law_ = new Law::MultivariateNormal<Point>(mean_, cov_);
  else         p_law_->update();
  // compute log likelihood of the gaussian law
  logLikelihood_ = static_cast<Law::MultivariateNormal<Point>* >(p_law_)->logLikelihood<Vector, Matrix>(*p_data_ );
  // everything ok
  return true;
}

bool STK::GaussianStatModel::run ( Vector const & weights ) [virtual]

implementation of the weighted Gaussian statistical model

Parameters:

weights the weights of the samples

Returns:: true if no error occur and false otherwise.

Implements STK::IRunnerPtr2D< Real, Point, Vector, Matrix >.

Definition at line 81 of file STK_GaussianStatModel.cpp.

References compWeightedCovariance(), compWeightedMean(), cov_, STK::IStatModelBase::logLikelihood_, mean_, STK::IRunnerPtr2D< Real, Point, Vector, Matrix >::p_data_, STK::ITStatModel< Real, Point, Vector, Matrix >::p_law_, and STK::Law::ILawBase::update().

{
  // compute the mean
  compWeightedMean(weights);
  // compute the covariance matrix
  compWeightedCovariance(weights);
  // create p_law_ (will be deleted in base class)
  // update gaussian law (will be deleted in base class)
  if (!p_law_) p_law_ = new Law::MultivariateNormal<Point>(mean_, cov_);
  else         p_law_->update();
  // compute log likelihood of the gaussian law
  logLikelihood_ = static_cast<Law::MultivariateNormal<Point>* >(p_law_)->logLikelihood<Vector, Matrix>(*p_data_ );
  // everything ok
  return true;
}

Point const& STK::GaussianStatModel::mean ( ) const [inline]

get the empirical mean.

Returns:: the empirical mean

Definition at line 130 of file STK_GaussianStatModel.h.

References mean_.

{ return mean_;}

MatrixSquare const& STK::GaussianStatModel::covariance ( ) const [inline]

get the empirical covariance

Returns:: the empirical covariance

Definition at line 134 of file STK_GaussianStatModel.h.

References cov_.

{ return cov_;}

void STK::GaussianStatModel::compMean ( ) [protected, virtual]

compute the empirical mean

Definition at line 98 of file STK_GaussianStatModel.cpp.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::col(), STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), mean_, STK::IRunnerPtr2D< Real, Point, Vector, Matrix >::p_data_, STK::IContainer2D::rangeHo(), and STK::IContainer1D::resize().

Referenced by run().

{
  // resize mean
  mean_.resize(p_data_->rangeHo());
  // get dimensions
  const Integer first = p_data_->firstCol(), last = p_data_->lastCol();
  for (Integer j= first; j <= last; ++j)
  {
    mean_[j] = Stat::mean<Vector>(p_data_->col(j));
  }
}

void STK::GaussianStatModel::compCovariance ( ) [protected, virtual]

compute the empirical covariance matrix.

Reimplemented in STK::Gaussian2BlocksStatModel.

Definition at line 110 of file STK_GaussianStatModel.cpp.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::col(), cov_, STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), mean_, STK::IRunnerPtr2D< Real, Point, Vector, Matrix >::p_data_, STK::IContainer2D::rangeHo(), and STK::MatrixSquare::resize().

Referenced by run().

{
  // resize mean
  cov_.resize(p_data_->rangeHo());
  // get dimensions
  const Integer first = p_data_->firstCol(), last = p_data_->lastCol();
  for (Integer i= first; i <= last; ++i)
  {
    cov_(i, i) = Stat::varianceWithFixedMean<Vector>(p_data_->col(i), mean_[i]);
    for (Integer j= first; j < i; ++j)
    {
      cov_(i, j) = Stat::covarianceWithFixedMean<Vector>(p_data_->col(i), p_data_->col(j), mean_[i], mean_[j]);
      cov_(j, i) = cov_(i,j);
    }
  }
}

void STK::GaussianStatModel::compWeightedMean ( Vector const & weights ) [protected, virtual]

compute the empirical weighted mean

Parameters:

weights the weights of the samples

Definition at line 128 of file STK_GaussianStatModel.cpp.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::col(), STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), mean_, STK::IRunnerPtr2D< Real, Point, Vector, Matrix >::p_data_, STK::IContainer2D::rangeHo(), and STK::IContainer1D::resize().

Referenced by run().

{
  // resize mean
  mean_.resize(p_data_->rangeHo());
  // get dimensions
  const Integer firstCol = p_data_->firstCol(), lastCol = p_data_->lastCol();
  for (Integer j= firstCol; j <= lastCol; ++j)
  {
    mean_[j] = Stat::mean<Vector>(p_data_->col(j), weights);
  }
}

void STK::GaussianStatModel::compWeightedCovariance ( Vector const & weights ) [protected, virtual]

compute the empirical weighted covariance matrix.

Parameters:

weights the weights of the samples

Reimplemented in STK::Gaussian2BlocksStatModel.

Definition at line 140 of file STK_GaussianStatModel.cpp.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::col(), cov_, STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), mean_, STK::IRunnerPtr2D< Real, Point, Vector, Matrix >::p_data_, STK::IContainer2D::rangeHo(), and STK::MatrixSquare::resize().

Referenced by run().

{
  // resize mean
  cov_.resize(p_data_->rangeHo());
  // get dimensions
  const Integer first = p_data_->firstCol(), last = p_data_->lastCol();
  for (Integer i= first; i <= last; ++i)
  {
    cov_(i, i) = Stat::varianceWithFixedMean<Vector>(p_data_->col(i), weights, mean_[i]);
    for (Integer j= first; j < i; ++j)
    {
      cov_(i, j) = Stat::covarianceWithFixedMean<Vector>(p_data_->col(i), p_data_->col(j), weights, mean_[i], mean_[j]);
      cov_(j, i) = cov_(i,j);
    }
  }
}