[ceph.git] / ceph / src / boost / libs / math / test / linear_regression_test.cpp

/*
 * Copyright Nick Thompson, 2019
 * Copyright Matt Borland, 2021
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0. (See accompanying file
 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 */

#include "math_unit_test.hpp"
#include <vector>
#include <random>
#include <utility>
#include <tuple>
#include <algorithm>
#include <boost/math/statistics/linear_regression.hpp>

using boost::math::statistics::simple_ordinary_least_squares;
using boost::math::statistics::simple_ordinary_least_squares_with_R_squared;

template<typename Real>
void test_line()
{
    std::vector<Real> x(128);
    std::vector<Real> y(128);
    Real expected_c0 = 7;
    Real expected_c1 = 12;
    for (size_t i = 0; i < x.size(); ++i) {
        x[i] = i;
        y[i] = expected_c0 + expected_c1*x[i];
    }

    std::pair<Real, Real> temp = simple_ordinary_least_squares(x, y);
    Real computed_c0 = std::get<0>(temp);
    Real computed_c1 = std::get<1>(temp);

    CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);

    std::tuple<Real, Real, Real> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
    Real computed_c0_R = std::get<0>(temp_with_R);
    Real computed_c1_R = std::get<1>(temp_with_R);
    Real Rsquared = std::get<2>(temp_with_R);

    Real expected_Rsquared = 1;
    CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
    CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);

}

template<typename Z>
void test_integer_line()
{
    std::vector<Z> x(128);
    std::vector<Z> y(128);
    double expected_c0 = 7;
    double expected_c1 = 12;
    for (size_t i = 0; i < x.size(); ++i) {
        x[i] = i;
        y[i] = expected_c0 + expected_c1*x[i];
    }

    std::pair<double, double> temp = simple_ordinary_least_squares(x, y);
    double computed_c0 = std::get<0>(temp);
    double computed_c1 = std::get<1>(temp);

    CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);

    std::tuple<double, double, double> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
    double computed_c0_R = std::get<0>(temp_with_R);
    double computed_c1_R = std::get<1>(temp_with_R);
    double Rsquared = std::get<2>(temp_with_R);

    double expected_Rsquared = 1;
    CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
    CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);

}

template<typename Real>
void test_constant()
{
    std::vector<Real> x(128);
    std::vector<Real> y(128);
    Real expected_c0 = 7;
    Real expected_c1 = 0;
    for (size_t i = 0; i < x.size(); ++i) {
        x[i] = i;
        y[i] = expected_c0 + expected_c1*x[i];
    }

    std::pair<Real, Real> temp = simple_ordinary_least_squares(x, y);
    Real computed_c0 = std::get<0>(temp);
    Real computed_c1 = std::get<1>(temp);

    CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);

    std::tuple<Real, Real, Real> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
    Real computed_c0_R = std::get<0>(temp_with_R);
    Real computed_c1_R = std::get<1>(temp_with_R);
    Real Rsquared = std::get<2>(temp_with_R);

    Real expected_Rsquared = 1;
    CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
    CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);

}

template<typename Z>
void test_integer_constant()
{
    std::vector<Z> x(128);
    std::vector<Z> y(128);
    double expected_c0 = 7;
    double expected_c1 = 0;
    for (size_t i = 0; i < x.size(); ++i) {
        x[i] = i;
        y[i] = expected_c0 + expected_c1*x[i];
    }

    std::pair<double, double> temp = simple_ordinary_least_squares(x, y);
    double computed_c0 = std::get<0>(temp);
    double computed_c1 = std::get<1>(temp);

    CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);

    std::tuple<double, double, double> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
    double computed_c0_R = std::get<0>(temp_with_R);
    double computed_c1_R = std::get<1>(temp_with_R);
    double Rsquared = std::get<2>(temp_with_R);

    double expected_Rsquared = 1;
    CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
    CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
    CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);

}

template<typename Real>
void test_permutation_invariance()
{
    std::vector<Real> x(256);
    std::vector<Real> y(256);
    std::mt19937_64 gen{123456};
    std::normal_distribution<Real> dis(0, 0.1);

    Real expected_c0 = -7.2;
    Real expected_c1 = -13.5;

    x[0] = 0;
    y[0] = expected_c0 + dis(gen);
    for(size_t i = 1; i < x.size(); ++i) {
        Real t = dis(gen);
        x[i] = x[i-1] + t*t;
        y[i] = expected_c0 + expected_c1*x[i] + dis(gen);
    }

    std::tuple<Real, Real, Real> temp = simple_ordinary_least_squares_with_R_squared(x, y);
    Real c0 = std::get<0>(temp);
    Real c1 = std::get<1>(temp);
    Real Rsquared = std::get<2>(temp);
    CHECK_MOLLIFIED_CLOSE(expected_c0, c0, 0.002);
    CHECK_MOLLIFIED_CLOSE(expected_c1, c1, 0.002);

    int j = 0;
    std::mt19937_64 gen1{12345};
    std::mt19937_64 gen2{12345};
    while(j++ < 10) {
        std::shuffle(x.begin(), x.end(), gen1);
        std::shuffle(y.begin(), y.end(), gen2);

        std::tuple<Real, Real, Real> temp_ = simple_ordinary_least_squares_with_R_squared(x, y);
        Real c0_ = std::get<0>(temp_);
        Real c1_ = std::get<1>(temp_);
        Real Rsquared_ = std::get<2>(temp_);

        CHECK_ULP_CLOSE(c0, c0_, 100);
        CHECK_ULP_CLOSE(c1, c1_, 100);
        CHECK_ULP_CLOSE(Rsquared, Rsquared_, 65);
    }
}

template<typename Real>
void test_scaling_relations()
{
    std::vector<Real> x(256);
    std::vector<Real> y(256);
    std::mt19937_64 gen{123456};
    std::normal_distribution<Real> dis(0, 0.1);

    Real expected_c0 = 3.2;
    Real expected_c1 = -13.5;

    x[0] = 0;
    y[0] = expected_c0 + dis(gen);
    for(size_t i = 1; i < x.size(); ++i) {
        Real t = dis(gen);
        x[i] = x[i-1] + t*t;
        y[i] = expected_c0 + expected_c1*x[i] + dis(gen);
    }

    std::tuple<Real, Real, Real> temp = simple_ordinary_least_squares_with_R_squared(x, y);
    Real c0 = std::get<0>(temp);
    Real c1 = std::get<1>(temp);
    Real Rsquared = std::get<2>(temp);
    CHECK_MOLLIFIED_CLOSE(expected_c0, c0, 0.005);
    CHECK_MOLLIFIED_CLOSE(expected_c1, c1, 0.005);

    // If y -> lambda y, then c0 -> lambda c0 and c1 -> lambda c1.
    Real lambda = 6;

    for (auto& s : y) {
        s *= lambda;
    }

    temp = simple_ordinary_least_squares_with_R_squared(x, y);
    Real c0_lambda = std::get<0>(temp);
    Real c1_lambda = std::get<1>(temp);
    Real Rsquared_lambda = std::get<2>(temp);

    CHECK_ULP_CLOSE(lambda*c0, c0_lambda, 30);
    CHECK_ULP_CLOSE(lambda*c1, c1_lambda, 30);
    CHECK_ULP_CLOSE(Rsquared, Rsquared_lambda, 3);

    // If x -> lambda x, then c0 -> c0 and c1 -> c1/lambda
    for (auto& s : x) {
        s *= lambda;
    }
    // Put y back into it's original state:
    for (auto& s : y) {
        s /= lambda;
    }

    temp = simple_ordinary_least_squares_with_R_squared(x, y);
    Real c0_ = std::get<0>(temp);
    Real c1_ = std::get<1>(temp);
    Real Rsquared_ = std::get<2>(temp);

    CHECK_ULP_CLOSE(c0, c0_, 70);
    CHECK_ULP_CLOSE(c1, c1_*lambda, 50);
    CHECK_ULP_CLOSE(Rsquared, Rsquared_, 50);

}


int main()
{
    test_line<float>();
    test_line<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_line<long double>();
#endif
    test_integer_line<int>();
    test_integer_line<int32_t>();
    test_integer_line<int64_t>();
    test_integer_line<uint32_t>();

    test_constant<float>();
    test_constant<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_constant<long double>();
#endif
    test_integer_constant<int>();
    test_integer_constant<int32_t>();
    test_integer_constant<int64_t>();
    test_integer_constant<uint32_t>();

    test_permutation_invariance<float>();
    test_permutation_invariance<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_permutation_invariance<long double>();
#endif

    test_scaling_relations<float>();
    test_scaling_relations<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_scaling_relations<long double>();
#endif
    return boost::math::test::report_errors();
}
Commit	Line	Data
f67539c2 TL	1	/*
f67539c2 TL	2	* Copyright Nick Thompson, 2019
1e59de90	3	* Copyright Matt Borland, 2021
f67539c2 TL	4	* Use, modification and distribution are subject to the
	5	* Boost Software License, Version 1.0. (See accompanying file
	6	* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
	7	*/
	8
	9	#include "math_unit_test.hpp"
	10	#include <vector>
	11	#include <random>
1e59de90 TL	12	#include <utility>
	13	#include <tuple>
	14	#include <algorithm>
f67539c2 TL	15	#include <boost/math/statistics/linear_regression.hpp>
	16
	17	using boost::math::statistics::simple_ordinary_least_squares;
	18	using boost::math::statistics::simple_ordinary_least_squares_with_R_squared;
	19
	20	template<typename Real>
	21	void test_line()
	22	{
	23	std::vector<Real> x(128);
	24	std::vector<Real> y(128);
	25	Real expected_c0 = 7;
	26	Real expected_c1 = 12;
	27	for (size_t i = 0; i < x.size(); ++i) {
	28	x[i] = i;
	29	y[i] = expected_c0 + expected_c1*x[i];
	30	}
	31
1e59de90 TL	32	std::pair<Real, Real> temp = simple_ordinary_least_squares(x, y);
	33	Real computed_c0 = std::get<0>(temp);
	34	Real computed_c1 = std::get<1>(temp);
f67539c2 TL	35
	36	CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
	37	CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);
	38
1e59de90 TL	39	std::tuple<Real, Real, Real> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
	40	Real computed_c0_R = std::get<0>(temp_with_R);
	41	Real computed_c1_R = std::get<1>(temp_with_R);
	42	Real Rsquared = std::get<2>(temp_with_R);
f67539c2 TL	43
f67539c2 TL	44	Real expected_Rsquared = 1;
1e59de90 TL	45	CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
	46	CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
	47	CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);
	48
	49	}
	50
	51	template<typename Z>
	52	void test_integer_line()
	53	{
	54	std::vector<Z> x(128);
	55	std::vector<Z> y(128);
	56	double expected_c0 = 7;
	57	double expected_c1 = 12;
	58	for (size_t i = 0; i < x.size(); ++i) {
	59	x[i] = i;
	60	y[i] = expected_c0 + expected_c1*x[i];
	61	}
	62
	63	std::pair<double, double> temp = simple_ordinary_least_squares(x, y);
	64	double computed_c0 = std::get<0>(temp);
	65	double computed_c1 = std::get<1>(temp);
	66
f67539c2 TL	67	CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
f67539c2 TL	68	CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);
1e59de90 TL	69
	70	std::tuple<double, double, double> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
	71	double computed_c0_R = std::get<0>(temp_with_R);
	72	double computed_c1_R = std::get<1>(temp_with_R);
	73	double Rsquared = std::get<2>(temp_with_R);
	74
	75	double expected_Rsquared = 1;
	76	CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
	77	CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
f67539c2 TL	78	CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);
	79
	80	}
	81
	82	template<typename Real>
	83	void test_constant()
	84	{
	85	std::vector<Real> x(128);
	86	std::vector<Real> y(128);
	87	Real expected_c0 = 7;
	88	Real expected_c1 = 0;
	89	for (size_t i = 0; i < x.size(); ++i) {
	90	x[i] = i;
	91	y[i] = expected_c0 + expected_c1*x[i];
	92	}
	93
1e59de90 TL	94	std::pair<Real, Real> temp = simple_ordinary_least_squares(x, y);
	95	Real computed_c0 = std::get<0>(temp);
	96	Real computed_c1 = std::get<1>(temp);
f67539c2 TL	97
	98	CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
	99	CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);
	100
1e59de90 TL	101	std::tuple<Real, Real, Real> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
	102	Real computed_c0_R = std::get<0>(temp_with_R);
	103	Real computed_c1_R = std::get<1>(temp_with_R);
	104	Real Rsquared = std::get<2>(temp_with_R);
f67539c2 TL	105
f67539c2 TL	106	Real expected_Rsquared = 1;
1e59de90 TL	107	CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
	108	CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
	109	CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);
	110
	111	}
	112
	113	template<typename Z>
	114	void test_integer_constant()
	115	{
	116	std::vector<Z> x(128);
	117	std::vector<Z> y(128);
	118	double expected_c0 = 7;
	119	double expected_c1 = 0;
	120	for (size_t i = 0; i < x.size(); ++i) {
	121	x[i] = i;
	122	y[i] = expected_c0 + expected_c1*x[i];
	123	}
	124
	125	std::pair<double, double> temp = simple_ordinary_least_squares(x, y);
	126	double computed_c0 = std::get<0>(temp);
	127	double computed_c1 = std::get<1>(temp);
	128
f67539c2 TL	129	CHECK_ULP_CLOSE(expected_c0, computed_c0, 0);
f67539c2 TL	130	CHECK_ULP_CLOSE(expected_c1, computed_c1, 0);
1e59de90 TL	131
	132	std::tuple<double, double, double> temp_with_R = simple_ordinary_least_squares_with_R_squared(x, y);
	133	double computed_c0_R = std::get<0>(temp_with_R);
	134	double computed_c1_R = std::get<1>(temp_with_R);
	135	double Rsquared = std::get<2>(temp_with_R);
	136
	137	double expected_Rsquared = 1;
	138	CHECK_ULP_CLOSE(expected_c0, computed_c0_R, 0);
	139	CHECK_ULP_CLOSE(expected_c1, computed_c1_R, 0);
f67539c2 TL	140	CHECK_ULP_CLOSE(expected_Rsquared, Rsquared, 0);
	141
	142	}
	143
	144	template<typename Real>
	145	void test_permutation_invariance()
	146	{
	147	std::vector<Real> x(256);
	148	std::vector<Real> y(256);
	149	std::mt19937_64 gen{123456};
	150	std::normal_distribution<Real> dis(0, 0.1);
	151
	152	Real expected_c0 = -7.2;
	153	Real expected_c1 = -13.5;
	154
	155	x[0] = 0;
	156	y[0] = expected_c0 + dis(gen);
	157	for(size_t i = 1; i < x.size(); ++i) {
	158	Real t = dis(gen);
	159	x[i] = x[i-1] + t*t;
	160	y[i] = expected_c0 + expected_c1*x[i] + dis(gen);
	161	}
	162
1e59de90 TL	163	std::tuple<Real, Real, Real> temp = simple_ordinary_least_squares_with_R_squared(x, y);
	164	Real c0 = std::get<0>(temp);
	165	Real c1 = std::get<1>(temp);
	166	Real Rsquared = std::get<2>(temp);
f67539c2 TL	167	CHECK_MOLLIFIED_CLOSE(expected_c0, c0, 0.002);
	168	CHECK_MOLLIFIED_CLOSE(expected_c1, c1, 0.002);
	169
	170	int j = 0;
	171	std::mt19937_64 gen1{12345};
	172	std::mt19937_64 gen2{12345};
	173	while(j++ < 10) {
	174	std::shuffle(x.begin(), x.end(), gen1);
	175	std::shuffle(y.begin(), y.end(), gen2);
1e59de90 TL	176
	177	std::tuple<Real, Real, Real> temp_ = simple_ordinary_least_squares_with_R_squared(x, y);
	178	Real c0_ = std::get<0>(temp_);
	179	Real c1_ = std::get<1>(temp_);
	180	Real Rsquared_ = std::get<2>(temp_);
f67539c2 TL	181
	182	CHECK_ULP_CLOSE(c0, c0_, 100);
	183	CHECK_ULP_CLOSE(c1, c1_, 100);
	184	CHECK_ULP_CLOSE(Rsquared, Rsquared_, 65);
	185	}
	186	}
	187
	188	template<typename Real>
	189	void test_scaling_relations()
	190	{
	191	std::vector<Real> x(256);
	192	std::vector<Real> y(256);
	193	std::mt19937_64 gen{123456};
	194	std::normal_distribution<Real> dis(0, 0.1);
	195
	196	Real expected_c0 = 3.2;
	197	Real expected_c1 = -13.5;
	198
	199	x[0] = 0;
	200	y[0] = expected_c0 + dis(gen);
	201	for(size_t i = 1; i < x.size(); ++i) {
	202	Real t = dis(gen);
	203	x[i] = x[i-1] + t*t;
	204	y[i] = expected_c0 + expected_c1*x[i] + dis(gen);
	205	}
	206
1e59de90 TL	207	std::tuple<Real, Real, Real> temp = simple_ordinary_least_squares_with_R_squared(x, y);
	208	Real c0 = std::get<0>(temp);
	209	Real c1 = std::get<1>(temp);
	210	Real Rsquared = std::get<2>(temp);
f67539c2 TL	211	CHECK_MOLLIFIED_CLOSE(expected_c0, c0, 0.005);
	212	CHECK_MOLLIFIED_CLOSE(expected_c1, c1, 0.005);
	213
	214	// If y -> lambda y, then c0 -> lambda c0 and c1 -> lambda c1.
	215	Real lambda = 6;
	216
	217	for (auto& s : y) {
	218	s *= lambda;
	219	}
	220
1e59de90 TL	221	temp = simple_ordinary_least_squares_with_R_squared(x, y);
	222	Real c0_lambda = std::get<0>(temp);
	223	Real c1_lambda = std::get<1>(temp);
	224	Real Rsquared_lambda = std::get<2>(temp);
f67539c2 TL	225
	226	CHECK_ULP_CLOSE(lambda*c0, c0_lambda, 30);
	227	CHECK_ULP_CLOSE(lambda*c1, c1_lambda, 30);
	228	CHECK_ULP_CLOSE(Rsquared, Rsquared_lambda, 3);
	229
	230	// If x -> lambda x, then c0 -> c0 and c1 -> c1/lambda
	231	for (auto& s : x) {
	232	s *= lambda;
	233	}
	234	// Put y back into it's original state:
	235	for (auto& s : y) {
	236	s /= lambda;
	237	}
f67539c2	238
1e59de90 TL	239	temp = simple_ordinary_least_squares_with_R_squared(x, y);
	240	Real c0_ = std::get<0>(temp);
	241	Real c1_ = std::get<1>(temp);
	242	Real Rsquared_ = std::get<2>(temp);
	243
	244	CHECK_ULP_CLOSE(c0, c0_, 70);
f67539c2 TL	245	CHECK_ULP_CLOSE(c1, c1_*lambda, 50);
	246	CHECK_ULP_CLOSE(Rsquared, Rsquared_, 50);
	247
	248	}
	249
	250
	251	int main()
	252	{
	253	test_line<float>();
	254	test_line<double>();
1e59de90	255	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
f67539c2	256	test_line<long double>();
1e59de90 TL	257	#endif
	258	test_integer_line<int>();
	259	test_integer_line<int32_t>();
	260	test_integer_line<int64_t>();
	261	test_integer_line<uint32_t>();
f67539c2 TL	262
	263	test_constant<float>();
	264	test_constant<double>();
1e59de90	265	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
f67539c2	266	test_constant<long double>();
1e59de90 TL	267	#endif
	268	test_integer_constant<int>();
	269	test_integer_constant<int32_t>();
	270	test_integer_constant<int64_t>();
	271	test_integer_constant<uint32_t>();
f67539c2 TL	272
	273	test_permutation_invariance<float>();
	274	test_permutation_invariance<double>();
1e59de90	275	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
f67539c2	276	test_permutation_invariance<long double>();
1e59de90	277	#endif
f67539c2 TL	278
	279	test_scaling_relations<float>();
	280	test_scaling_relations<double>();
1e59de90	281	#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
f67539c2	282	test_scaling_relations<long double>();
1e59de90	283	#endif
f67539c2 TL	284	return boost::math::test::report_errors();
f67539c2 TL	285	}