|| || || || || || || ||

---

MatCovariance.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatCovariance_h
#define Impala_Core_Matrix_MatCovariance_h

#include "Core/Matrix/MatSet.h"
#include "Core/Matrix/MatDiv.h"
#include "Core/Matrix/VecFunc.h"

namespace Impala
{
namespace Core
{
namespace Matrix
{


template<class ArrayT>
inline ArrayT*
MatCovariance(ArrayT* m, ArrayT* meanOut = 0)
{
    typedef typename ArrayT::ArithType ArithT;
    typedef typename ArrayT::StorType StorT;

    int nr = MatNrRow(m);
    int nc = MatNrCol(m);
    ArrayT* cov = MatCreate<ArrayT>(nc, nc);
    MatSet(cov, 0);
    ArrayT* mean = VecCreate<ArrayT>(nc) ; //means
    MatSet(mean, 0);

    for (int k=0; k<nr; k++)
    {
        StorT* meanPtr = VecE(mean, 0);
        StorT* mkPtr = MatE(m, k, 0);
        for (int i=0; i<nc; i++)
        {
            ArithT mkiVal = mkPtr[i];
            meanPtr[i] += mkiVal;
            StorT* coviPtr = MatE(cov, i, 0);
            for (int j=i; j<nc; j++)
            {
                coviPtr[j] += mkiVal * mkPtr[j];
            }
        }
    }

    MatDiv(mean, nr);

    for (int i=0; i<nc; i++)
    {
        StorT* coviPtr = MatE(cov, i, 0);
        StorT* meanPtr = VecE(mean, 0);
        ArithT meani = meanPtr[i];
        for (int j=i; j<nc; j++)
        {
            ArithT res = coviPtr[j] / nr - meani * meanPtr[j];
            coviPtr[j] = res;
            *MatE(cov, j, i) = res;
        }
    }

    if (meanOut != 0)
        MatSet(meanOut, mean);
    delete mean;

    return cov;
}

template<class ArrayT>
inline ArrayT*
MatCovarianceExact(ArrayT* m, ArrayT* meanOut = 0)
{
    typedef typename ArrayT::ArithType ArithT;
    typedef typename ArrayT::StorType StorT;

    int nr = MatNrRow(m);
    int nc = MatNrCol(m);
    ArrayT* cov = MatCreate<ArrayT>(nc, nc);
    MatSet(cov, 0);
    ArrayT* mean = VecCreate<ArrayT>(nc) ; //means
    MatSet(mean, 0);

    int k;
    for ( k=0; k<nr; k++)
    {
        StorT* meanPtr = VecE(mean, 0);
        StorT* mkPtr = MatE(m, k, 0);
        for (int i=0; i<nc; i++)
        {
            ArithT mkiVal = mkPtr[i];
            meanPtr[i] += mkiVal;
            /*
            StorT* coviPtr = MatE(cov, i, 0);
            for (int j=i; j<nc; j++) {
                coviPtr[j] += mkiVal * mkPtr[j];
            }
            */
        }
    }

    MatDiv(mean, nr);

    for ( k=0; k<nr; k++) {
        StorT* meanPtr = VecE(mean, 0);
        StorT* mkPtr = MatE(m, k, 0);
        for (int i=0; i<nc; i++)
        {
            //ArithT mkiVal = mkPtr[i];
            StorT* coviPtr = MatE(cov, i, 0);
            for (int j=i; j<nc; j++)
            {
                coviPtr[j] += (mkPtr[i] - meanPtr[i]) * (mkPtr[j] - meanPtr[j]);
            }
        }
    }

    for (int i=0; i<nc; i++) {
        StorT* coviPtr = MatE(cov, i, 0);
        //StorT* meanPtr = VecE(mean, 0);
        //ArithT meani = meanPtr[i];
        for (int j=i; j<nc; j++)
        {
            ArithT res = coviPtr[j] / (nr - 1);
            coviPtr[j] = res;
            *MatE(cov, j, i) = res;
        }
    }

    if (meanOut != 0)
        MatSet(meanOut, mean);
    delete mean;

    return cov;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---