statistics.h

Go to the documentation of this file.

/******************************************************************************
 * Copyright 2014 Micah C Chambers (micahc.vt@gmail.com)
 *
 * NPL is free software: you can redistribute it and/or modify it under the
 * terms of the BSD 2-Clause License available in LICENSE or at
 * http://opensource.org/licenses/BSD-2-Clause
 *
 * @file statistics.h Tools for analyzing data, including PCA, ICA and
 * general linear modeling.
 *
 *****************************************************************************/

#ifndef STATISTICS_H
#define STATISTICS_H

#include <Eigen/Dense>
#include "npltypes.h"
#include "mrimage.h"

namespace npl {

double sample_var(const Ref<const VectorXd> vec);

inline
void regressOutLS(VectorXd& signal, const MatrixXd& X, const MatrixXd& covInv)
{
    VectorXd beta = covInv*X.transpose()*signal;
    signal -= X*beta;
}

template <typename T>
void fillGaussian(Ref<T> m)
{
    static std::random_device rd;
    static std::default_random_engine rng(rd());
    std::normal_distribution<double> rdist(0,1);

    for(size_t cc=0; cc<m.cols(); cc++) {
        for(size_t rr=0; rr<m.rows(); rr++) {
            m(rr, cc) = rdist(rng);
        }
    }
}


/******************************************
 * \defgroup Basic Stats Functions
 * @{
 *****************************************/

double gaussianPDF(double mean, double sd, double x);

double gaussianCDF(double mean, double sd, double x);

double gammaPDF_MS(double mean, double sd, double x);

double mutualInformation(size_t len, double* a, double* b, size_t mbin);

double correlation(size_t len, double* a, double* b);

inline
double sample_var(int count, double sum, double sumsqr)
{
    return (sumsqr-sum*sum/count)/(count-1);
}

inline
double sample_corr(int count, double sum1, double sum2,
        double sumsq1, double sumsq2, double s1s2)
{
    return (count*s1s2-sum1*sum2)/
            sqrt((count*sumsq1-sum1*sum1)*(count*sumsq2-sum2*sum2));
}

struct RegrResult
{
    VectorXd yhat;

    VectorXd bhat;

    double ssres;

    double sigmahat;

    double rsqr;

    double adj_rsqr;

    VectorXd std_err;

    double dof;

    VectorXd t;

    VectorXd p;
};

class StudentsT
{
public:
    StudentsT(int dof = 2, double dt = 0.1, double tmax = 20);

    void setDOF(double dof);

    void setStepT(double dt);

    void setMaxT(double tmax);

    double cumulative(double t) const;

    double cdf(double t) const { return cumulative(t); };

    double density(double t) const;

    double pdf(double t) const { return density(t); };

    double icdf(double t) const;

    double tthresh(double p) const { return icdf(p); };

private:
    void init();

    double m_dt;
    double m_tmax;
    int m_dof;
    std::vector<double> m_cdf;
    std::vector<double> m_pdf;
    std::vector<double> m_tvals;
};


void regress(RegrResult* out,
        const Ref<const VectorXd> y,
        const Ref<const MatrixXd> X,
        const Ref<const VectorXd> covInv,
        const Ref<const MatrixXd> Xinv,
        const StudentsT& distrib);

void regress(RegrResult* out,
        const Ref<const VectorXd> y,
        const Ref<const MatrixXd> X);

/******************************************
 * @}
 * \defgroup Matrix Decompositions
 * @{
 *****************************************/

void randomizePowerIterationSVD(const Ref<const MatrixXd> A,
        size_t subsize, size_t poweriters, MatrixXd& U, VectorXd& E,
        MatrixXd& V);

void randomizePowerIterationSVD(const Ref<const MatrixXd> A,
        double tol, size_t startrank, size_t maxrank, size_t poweriters,
        MatrixXd& U, VectorXd& E, MatrixXd& V);

MatrixXd pca(const Ref<const MatrixXd> X, double varth = 1, int odim = -1);

MatrixXd rpiPCA(const Ref<const MatrixXd> X, double varth, int odim);

MatrixXd symICA(const Ref<const MatrixXd> Xin, MatrixXd* unmix = NULL);

MatrixXd asymICA(const Ref<const MatrixXd> Xin, MatrixXd* unmix = NULL);

MatrixXd pseudoInverse(const Ref<const MatrixXd> X);

VectorXd shootingRegr(const Ref<const MatrixXd> X,
        const Ref<const VectorXd> y, double gamma);

VectorXd activeShootingRegr(const Ref<const MatrixXd> X,
        const Ref<const VectorXd> y, double gamma);


MatrixXd pcacov(const Ref<const MatrixXd> cov, double varth);

/* @}
 * \defgroup Clustering algorithms
 * @{
 */

void expMax1D(const Ref<const VectorXd> data,
        vector<std::function<double(double,double,double)>> pdfs,
        Ref<VectorXd> mean, Ref<VectorXd> sd, Ref<VectorXd> prior,
        std::string plotfile = "");

void gaussGammaMixtureModel(const Ref<const VectorXd> data,
        Ref<VectorXd> mu, Ref<VectorXd> sd, Ref<VectorXd> prior,
        std::string plotfile);

void approxKMeans(const Ref<const MatrixXd> samples,
        size_t nclass, MatrixXd& means);

class Classifier
{
public:
    Classifier(size_t rank) : ndim(rank), maxit(-1), m_valid(false) {};

    virtual
    Eigen::VectorXi classify(const Ref<const MatrixXd> samples) = 0;

    virtual
    size_t classify(const Ref<const MatrixXd> samples, Ref<VectorXi> oclass) = 0;

    virtual
        int update(const Ref<const MatrixXd> samples, bool reinit = false) = 0;

    inline
    void compute(const Ref<const MatrixXd> samples)
    {
        update(samples, true);
    };

    const int ndim;

    int maxit;
protected:
    bool m_valid;

};

class KMeans : public Classifier
{
public:
    KMeans(size_t rank, size_t k = 2);

    void setk(size_t ngroups);

    void updateMeans(const Ref<const MatrixXd> newmeans);

    void updateMeans(const Ref<const MatrixXd> samples,
            const Ref<const Eigen::VectorXi> classes);

    VectorXi classify(const Ref<const MatrixXd> samples);

    size_t classify(const Ref<const MatrixXd> samples, Ref<VectorXi> oclass);

    int update(const Ref<const MatrixXd> samples, bool reinit = false);

    const MatrixXd& getMeans() { return m_mu; };

private:
    size_t m_k;

    MatrixXd m_mu;
};

class ExpMax : public Classifier
{
public:
    ExpMax(size_t rank, size_t k = 2);

    void setk(size_t ngroups);

    void updateMeanCovTau(const Ref<const MatrixXd> newmeans, const Ref<const MatrixXd> newcovs,
            const Ref<const VectorXd> tau);

    void updateMeanCovTau(const Ref<const MatrixXd> samples, Ref<MatrixXd> prob);

    double expectation(const Ref<const MatrixXd> samples, Ref<MatrixXd> prob);

    Eigen::VectorXi classify(const Ref<const MatrixXd> samples);

    size_t classify(const Ref<const MatrixXd> samples, Ref<VectorXi> oclass);

    int update(const Ref<const MatrixXd> samples, bool reinit = false);


    const MatrixXd& getMeans() { return m_mu; };

    const MatrixXd& getCovs() { return m_cov; };
private:
    size_t m_k;

    MatrixXd m_mu;

    MatrixXd m_cov;

    MatrixXd m_covinv;

    VectorXd m_tau;

    double m_ll;
};

int fastSearchFindDP(const Eigen::MatrixXf& samples, double thresh,
        double outthresh, Eigen::VectorXi& classes, bool brute = false);

int findDensityPeaks(const Eigen::MatrixXf& samples, double thresh,
        Eigen::VectorXf& rho, Eigen::VectorXf& delta,
        Eigen::VectorXi& parent);

int findDensityPeaks_brute(const Eigen::MatrixXf& samples, double thresh,
        Eigen::VectorXf& rho, Eigen::VectorXf& delta,
        Eigen::VectorXi& parent);

}

#endif // STATISTICS_H