Commit dcf595fe authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'math-linear-regression' into 'master'

Math: add a function for linear regression

See merge request !2569
parents 90d98043 894b8adf
......@@ -25,6 +25,7 @@
#define DUMUX_MATH_HH
#include <algorithm>
#include <numeric>
#include <cmath>
#include <utility>
......@@ -881,6 +882,44 @@ vtmv(const Dune::DenseVector<V1>& v1,
return m*(v1*v2);
}
/*!
* \ingroup Common
* \brief Returns the [intercept, slope] of the regression line
* fitted to a set of (x, y) data points.
*
* Note: We use least-square regression method to find
* the regression line.
*
* \param x x-values of the data set
* \param y y-values of the data set
*/
template <class Scalar>
std::array<Scalar, 2> linearRegression(const std::vector<Scalar>& x,
const std::vector<Scalar>& y)
{
if (x.size() != y.size())
DUNE_THROW(Dune::InvalidStateException, "x and y array must have the same length.");
const Scalar averageX = std::accumulate(x.begin(), x.end(), 0.0)/x.size();
const Scalar averageY = std::accumulate(y.begin(), y.end(), 0.0)/y.size();
// calculate temporary variables necessary for slope computation
const Scalar numerator = std::inner_product(
x.begin(), x.end(), y.begin(), 0.0, std::plus<Scalar>(),
[&](auto xx, auto yy) { return (xx - averageX) * (yy - averageY); }
);
const Scalar denominator = std::inner_product(
x.begin(), x.end(), x.begin(), 0.0, std::plus<Scalar>(),
[&](auto xx, auto yy) { return (xx - averageX) * (yy - averageX); }
);
// compute slope and intercept of the regression line
const Scalar slope = numerator / denominator;
const Scalar intercept = averageY - slope * averageX;
return {intercept, slope};
}
} // end namespace Dumux
#endif
......@@ -214,4 +214,25 @@ int main()
DUNE_THROW(Dune::Exception, "[linspace] Not sorted in ascending order!");
}
//////////////////////////////////////////////////////////////////
///// Dumux::linearRegression ////////////////////////////////////
//////////////////////////////////////////////////////////////////
{
const auto x = Dumux::linspace(0.0, 1.0, 100);
const auto [intercept, slope] = Dumux::linearRegression(x, x);
if (std::abs(intercept) > 1e-14)
DUNE_THROW(Dune::Exception, "[linearRegreesion] Wrong intercept " << intercept << ", should be 0.0.");
if (std::abs(slope-1.0) > 1e-14)
DUNE_THROW(Dune::Exception, "[linearRegreesion] Wrong slope " << slope << ", should be 1.0.");
}
{
const auto x = Dumux::linspace(0.0, 1.0, 100);
const auto y = Dumux::linspace(1.0, 3.0, 100);
const auto [intercept, slope] = Dumux::linearRegression(x, y);
if (std::abs(intercept-1.0) > 1e-14)
DUNE_THROW(Dune::Exception, "[linearRegreesion] Wrong intercept " << intercept << ", should be 1.0.");
if (std::abs(slope-2.0) > 1e-14)
DUNE_THROW(Dune::Exception, "[linearRegreesion] Wrong slope " << slope << ", should be 2.0.");
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment