glm.hpp

Go to the documentation of this file.

/*
 * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
 * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
 * 
 * This file is part of SHORE.
 * 
 * SHORE is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * SHORE is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with SHORE.  If not, see <http://www.gnu.org/licenses/>.
 */


#ifndef SHORE_STATISTICS_GLM_HPP__
#define SHORE_STATISTICS_GLM_HPP__

#include <cmath>
#include <iostream>
#include <vector>

#include <gsl/gsl_multifit.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

namespace shore {

class glht;

class glm
{
 private:

        friend class glht;

        const size_t m_rows;
        const size_t m_cols;

        void *m_modelmatrix;

        void *m_workspace;

        void *m_observations;

        void *m_coefficients;
        void *m_covariance;
        void *m_residuals;

        double m_chisq;

        size_t m_rank;


        void alloc();

        void set_modelmatrix(size_t i,size_t j,double value);

        size_t modelmatrix_rank(const double eps=1e-7);

        void set_observation(size_t i,double x);

        void fit();

        double get_coefficient(size_t i) const;

        void calc_residuals();
        double get_residual(size_t i) const;

        double get_vcov(size_t i,size_t j) const;

 public:
        
        template<class T,size_t rows,size_t columns>
        glm(const T (&design)[rows][columns])
        :m_rows(rows),
         m_cols(columns)
        {
                alloc();

                for(size_t i=0;i<m_rows;++i)
                {
                        for(size_t j=0;j<m_cols;++j)
                        {
                                set_modelmatrix(i,j,design[i][j]);
                        }
                }

                m_rank=modelmatrix_rank(1e-7);

                std::cerr<<"GLM: columns="<<m_cols
                         <<", rank="<<m_rank<<", res. df="<<get_residual_df()<<std::endl;
        }

        glm(const std::vector<std::vector<double> >& design);

        ~glm();

        size_t get_design_rank() const;

        size_t get_residual_df() const;

        size_t nrows() const;

        size_t ncols() const;

        template<class Iter>
        void fit(Iter obs)
        {
                // init vector
                for(size_t i=0;i<m_rows;++i)
                        set_observation(i,*obs++);

                fit();
        }
        
        template<class Iter>
        void get_coefficients(Iter dest) const
        {
                for(size_t i=0;i<m_cols;++i)
                        (*dest++)=get_coefficient(i);
        }

        template<class Iter>
        void get_residuals(Iter dest)
        {
                calc_residuals();

                for(size_t i=0;i<m_rows;++i)
                        (*dest++)=get_residual(i);
        }

        template<class Twodim>
        void get_vcov(Twodim& res) const
        {
                for(size_t i=0;i<m_cols;++i)
                {
                        for(size_t j=0;j<m_cols;++j)
                        {
                                res[i][j]=get_vcov(i,j);
                        }
                }
        }
};

class glht
{
 private:

        const glm& m_glm;

        const size_t m_rows;
        const size_t m_cols;

        void *m_lh;

        void *m_tcross;
        void *m_covm;

        int m_alternative;

        // results:

        void *m_estimates;

        double *m_stderr;
        double *m_tstats;
        double *m_pvalues;


        void alloc();

        void mmprod(void *const A,void *const B,void *const C);

        void tcrossprod(void *const A,void *const B,void *const C);

        void mvprod(void *const A,void *const x,void *const y);

        void set_lhmatrix(size_t i,size_t j,double x);

        double get_estimate(size_t i) const;

 public:

        enum Alternative
        {
                ALT_TWOSIDED=0,
                ALT_LESS=1,
                ALT_GREATER=2
        };
        
        template<class T,size_t rows,size_t columns>
        glht(const glm*const m,const T (&lh)[rows][columns])
        :m_glm(*m),
         m_rows(rows),m_cols(columns),
         m_alternative(ALT_TWOSIDED)
        {
                if(m_glm.ncols()!=m_cols)
                        throw "column counts of glm and glht must match";

                alloc();

                for(size_t i=0;i<m_rows;++i)
                {
                        for(size_t j=0;j<m_cols;++j)
                        {
                                set_lhmatrix(i,j,lh[i][j]);
                        }
                }
        }

        glht(const glm*const m,const std::vector<std::vector<double> >& lh);

        ~glht();

        size_t nrows() const;

        size_t ncols() const;

        void set_alternative(const Alternative alt);

        void run();

        template<class Iter>
        void get_estimates(Iter dest)
        {
                for(size_t i=0;i<m_rows;++i)
                        (*dest++)=get_estimate(i);
        }

        template<class Iter>
        void get_stderrs(Iter dest)
        {
                for(size_t i=0;i<m_rows;++i)
                        (*dest++)=m_stderr[i];
        }

        template<class Iter>
        void get_tstats(Iter dest)
        {
                for(size_t i=0;i<m_rows;++i)
                        (*dest++)=m_tstats[i];
        }

        template<class Iter>
        void get_pvalues(Iter dest)
        {
                for(size_t i=0;i<m_rows;++i)
                        (*dest++)=m_pvalues[i];
        }
};

} // namespace

#endif // SHORE_STATISTICS_GLM_HPP__