LinearRegression.hpp

#pragma once
/* (C) 2020 Roman Werpachowski. */
#include <string>
#include <Eigen/Core>
#include "dll.hpp"
#include "Crossvalidation.hpp"
 
namespace ml 
{
    namespace LinearRegression {
 
        struct Result
        {
            unsigned int n; 
            unsigned int dof; 
            double rss; 
            double tss; 
            double var_y() const {
                if (dof) {
                    return rss / static_cast<double>(dof);
                }
                else {
                    return std::numeric_limits<double>::quiet_NaN();
                }
            }
 
            double r2() const {
                return 1 - rss / tss;
            }
            
            double adjusted_r2() const {
                if (dof) {
                    return 1 - (rss / static_cast<double>(dof)) / (tss / static_cast<double>(n - 1));
                }
                else {
                    return std::numeric_limits<double>::quiet_NaN();
                }
            }
        };
 
        struct UnivariateOLSResult: public Result
        {
            double slope; 
            double intercept; 
            double var_slope; 
            double var_intercept; 
            double cov_slope_intercept; 
            DLL_DECLSPEC std::string to_string() const;
 
            DLL_DECLSPEC Eigen::VectorXd predict(Eigen::Ref<const Eigen::VectorXd> x) const;
 
            double predict(double x) const
            {
                return x * slope + intercept;
            }
        };
 
        struct MultivariateOLSResult : public Result
        {
            Eigen::VectorXd beta; 
            Eigen::MatrixXd cov; 
            DLL_DECLSPEC std::string to_string() const;
 
            DLL_DECLSPEC Eigen::VectorXd predict(Eigen::Ref<const Eigen::MatrixXd> X) const;
 
            DLL_DECLSPEC double predict_single(Eigen::Ref<const Eigen::VectorXd> x) const;
        };
 
        struct RegularisedRegressionResult : public Result
        {
            Eigen::VectorXd beta; 
            double effective_dof; 
            DLL_DECLSPEC Eigen::VectorXd predict(Eigen::Ref<const Eigen::MatrixXd> X) const;
 
            DLL_DECLSPEC double predict_single(Eigen::Ref<const Eigen::VectorXd> x) const;
        };
 
        struct RidgeRegressionResult : public RegularisedRegressionResult
        {
            Eigen::MatrixXd cov;  
            DLL_DECLSPEC std::string to_string() const;
 
            using RegularisedRegressionResult::predict;         
        };
 
        struct LassoRegressionResult : public RegularisedRegressionResult
        {
            DLL_DECLSPEC std::string to_string() const;
 
            using RegularisedRegressionResult::predict;         
        };
 
        DLL_DECLSPEC UnivariateOLSResult univariate(Eigen::Ref<const Eigen::VectorXd> x, Eigen::Ref<const Eigen::VectorXd> y);
 
        DLL_DECLSPEC UnivariateOLSResult univariate(double x0, double dx, Eigen::Ref<const Eigen::VectorXd> y);
 
        DLL_DECLSPEC UnivariateOLSResult univariate_without_intercept(Eigen::Ref<const Eigen::VectorXd> x, Eigen::Ref<const Eigen::VectorXd> y);
 
        DLL_DECLSPEC MultivariateOLSResult multivariate(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y);
 
        template <bool DoStandardise> RidgeRegressionResult ridge(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda);
 
        template <> DLL_DECLSPEC RidgeRegressionResult ridge<true>(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda);
 
        template <> DLL_DECLSPEC RidgeRegressionResult ridge<false>(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda);
 
        inline RidgeRegressionResult ridge(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda, bool do_standardise)
        {
            if (do_standardise) {
                return ridge<true>(X, y, lambda);
            } else {
                return ridge<false>(X, y, lambda);
            }
        }
 
        template <bool DoStandardise> LassoRegressionResult lasso(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda);
 
        template <> DLL_DECLSPEC LassoRegressionResult lasso<true>(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda);
 
        template <> DLL_DECLSPEC LassoRegressionResult lasso<false>(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda);
 
        inline LassoRegressionResult lasso(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, double lambda, bool do_standardise)
        {
            if (do_standardise) {
                return lasso<true>(X, y, lambda);
            } else {
                return lasso<false>(X, y, lambda);
            }
        }
 
        template <class Regression> double press(Eigen::Ref<const Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> y, Regression regression)
        {
            const auto trainer = [&regression](const Eigen::Ref<const Eigen::MatrixXd> X, const Eigen::Ref<const Eigen::VectorXd> y) -> auto {
                return regression(X, y);
            };
            const auto tester = [](const decltype(regression(X, y))& result, const Eigen::Ref<const Eigen::MatrixXd> X, const Eigen::Ref<const Eigen::VectorXd> y) -> double {
                return (y - result.predict(X)).squaredNorm() / static_cast<double>(y.size());
            };
            return Crossvalidation::leave_one_out(X, y, trainer, tester) * static_cast<double>(y.size());
        }
 
        template <bool WithIntercept> double press_univariate(Eigen::Ref<const Eigen::VectorXd> x, Eigen::Ref<const Eigen::VectorXd> y)
        {
            const auto trainer = [](const Eigen::Ref<const Eigen::VectorXd> x, const Eigen::Ref<const Eigen::VectorXd> y) -> UnivariateOLSResult {
                if (WithIntercept) {
                    return univariate(x, y);
                } else {
                    return univariate_without_intercept(x, y);
                }
            };
            const auto tester = [](const UnivariateOLSResult& result, const Eigen::Ref<const Eigen::VectorXd> x, const Eigen::Ref<const Eigen::VectorXd> y) -> double {
                return (y - result.predict(x)).squaredNorm() / static_cast<double>(y.size());
            };
            return Crossvalidation::leave_one_out_scalar(x, y, trainer, tester) * static_cast<double>(y.size());
        }
 
 
 
        DLL_DECLSPEC Eigen::MatrixXd add_ones(Eigen::Ref<const Eigen::MatrixXd> X);
 
        DLL_DECLSPEC void standardise(Eigen::Ref<Eigen::MatrixXd> X);
 
        DLL_DECLSPEC void standardise(Eigen::Ref<Eigen::MatrixXd> X, Eigen::VectorXd& means, Eigen::VectorXd& standard_deviations);
 
        DLL_DECLSPEC void unstandardise(Eigen::Ref<Eigen::MatrixXd> X, Eigen::Ref<const Eigen::VectorXd> means, Eigen::Ref<const Eigen::VectorXd> standard_deviations);
 
        Eigen::VectorXd calculate_XXt_beta(const Eigen::Ref<const Eigen::MatrixXd> X, const Eigen::Ref<const Eigen::VectorXd> y, Eigen::Ref<Eigen::MatrixXd> XXt, Eigen::LDLT<Eigen::MatrixXd>& xxt_decomp, const Eigen::Ref<const Eigen::VectorXd> lambda);
    }
}