134 lines · plain
1/*2 * Copyright Nick Thompson, 20193 * Copyright Matt Borland, 20214 * Use, modification and distribution are subject to the5 * Boost Software License, Version 1.0. (See accompanying file6 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7 */8 9#ifndef BOOST_MATH_STATISTICS_LINEAR_REGRESSION_HPP10#define BOOST_MATH_STATISTICS_LINEAR_REGRESSION_HPP11 12#include <cmath>13#include <algorithm>14#include <utility>15#include <tuple>16#include <stdexcept>17#include <type_traits>18#include <boost/math/statistics/univariate_statistics.hpp>19#include <boost/math/statistics/bivariate_statistics.hpp>20 21namespace boost { namespace math { namespace statistics { namespace detail {22 23 24template<class ReturnType, class RandomAccessContainer>25ReturnType simple_ordinary_least_squares_impl(RandomAccessContainer const & x,26 RandomAccessContainer const & y)27{28 using Real = typename std::tuple_element<0, ReturnType>::type;29 if (x.size() <= 1)30 {31 throw std::domain_error("At least 2 samples are required to perform a linear regression.");32 }33 34 if (x.size() != y.size())35 {36 throw std::domain_error("The same number of samples must be in the independent and dependent variable.");37 }38 std::tuple<Real, Real, Real> temp = boost::math::statistics::means_and_covariance(x, y);39 Real mu_x = std::get<0>(temp);40 Real mu_y = std::get<1>(temp);41 Real cov_xy = std::get<2>(temp);42 43 Real var_x = boost::math::statistics::variance(x);44 45 if (var_x <= 0) {46 throw std::domain_error("Independent variable has no variance; this breaks linear regression.");47 }48 49 50 Real c1 = cov_xy/var_x;51 Real c0 = mu_y - c1*mu_x;52 53 return std::make_pair(c0, c1);54}55 56template<class ReturnType, class RandomAccessContainer>57ReturnType simple_ordinary_least_squares_with_R_squared_impl(RandomAccessContainer const & x,58 RandomAccessContainer const & y)59{60 using Real = typename std::tuple_element<0, ReturnType>::type;61 if (x.size() <= 1)62 {63 throw std::domain_error("At least 2 samples are required to perform a linear regression.");64 }65 66 if (x.size() != y.size())67 {68 throw std::domain_error("The same number of samples must be in the independent and dependent variable.");69 }70 std::tuple<Real, Real, Real> temp = boost::math::statistics::means_and_covariance(x, y);71 Real mu_x = std::get<0>(temp);72 Real mu_y = std::get<1>(temp);73 Real cov_xy = std::get<2>(temp);74 75 Real var_x = boost::math::statistics::variance(x);76 77 if (var_x <= 0) {78 throw std::domain_error("Independent variable has no variance; this breaks linear regression.");79 }80 81 82 Real c1 = cov_xy/var_x;83 Real c0 = mu_y - c1*mu_x;84 85 Real squared_residuals = 0;86 Real squared_mean_deviation = 0;87 for(decltype(y.size()) i = 0; i < y.size(); ++i) {88 squared_mean_deviation += (y[i] - mu_y)*(y[i]-mu_y);89 Real ei = (c0 + c1*x[i]) - y[i];90 squared_residuals += ei*ei;91 }92 93 Real Rsquared;94 if (squared_mean_deviation == 0) {95 // Then y = constant, so the linear regression is perfect.96 Rsquared = 1;97 } else {98 Rsquared = 1 - squared_residuals/squared_mean_deviation;99 }100 101 return std::make_tuple(c0, c1, Rsquared);102}103} // namespace detail104 105template<typename RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, 106 typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>107inline auto simple_ordinary_least_squares(RandomAccessContainer const & x, RandomAccessContainer const & y) -> std::pair<double, double>108{109 return detail::simple_ordinary_least_squares_impl<std::pair<double, double>>(x, y);110}111 112template<typename RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, 113 typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>114inline auto simple_ordinary_least_squares(RandomAccessContainer const & x, RandomAccessContainer const & y) -> std::pair<Real, Real>115{116 return detail::simple_ordinary_least_squares_impl<std::pair<Real, Real>>(x, y);117}118 119template<typename RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, 120 typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>121inline auto simple_ordinary_least_squares_with_R_squared(RandomAccessContainer const & x, RandomAccessContainer const & y) -> std::tuple<double, double, double>122{123 return detail::simple_ordinary_least_squares_with_R_squared_impl<std::tuple<double, double, double>>(x, y);124}125 126template<typename RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, 127 typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>128inline auto simple_ordinary_least_squares_with_R_squared(RandomAccessContainer const & x, RandomAccessContainer const & y) -> std::tuple<Real, Real, Real>129{130 return detail::simple_ordinary_least_squares_with_R_squared_impl<std::tuple<Real, Real, Real>>(x, y);131}132}}} // namespace boost::math::statistics133#endif134