21#ifndef JKQTPSTATREGRESSION_H_INCLUDED
22#define JKQTPSTATREGRESSION_H_INCLUDED
37#include "jkqtmath/jkqtmath_imexport.h"
38#include "jkqtmath/jkqtplinalgtools.h"
39#include "jkqtmath/jkqtparraytools.h"
40#include "jkqtcommon/jkqtpdebuggingtools.h"
41#include "jkqtmath/jkqtpstatbasics.h"
42#include "jkqtmath/jkqtpstatpoly.h"
70template <
class InputItX,
class InputItY>
71inline void jkqtpstatLinearRegression(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY,
double& coeffA,
double& coeffB,
bool fixA=
false,
bool fixB=
false) {
72 if (fixA&&fixB)
return;
73 const int Nx=std::distance(firstX,lastX);
74 const int Ny=std::distance(firstY,lastY);
78 double sumx=0, sumy=0, sumxy=0, sumx2=0;
82 for (; itX!=lastX && itY!=lastY; ++itX, ++itY) {
88 sumxy=sumxy+fit_x*fit_y;
89 sumx2=sumx2+fit_x*fit_x;
93 const double NN=
static_cast<double>(N);
96 coeffB=(double(sumxy)-double(sumx)*double(sumy)/NN)/(
double(sumx2)-double(sumx)*double(sumx)/NN);;
97 coeffA=double(sumy)/NN-coeffB*double(sumx)/NN;
98 }
else if (fixA && !fixB) {
99 coeffB=(double(sumy)/NN-coeffA)/(
double(sumx)/NN);
100 }
else if (!fixA && fixB) {
101 coeffA=double(sumy)/NN-coeffB*double(sumx)/NN;
143template <
class InputItX,
class InputItY,
class InputItW>
144inline void jkqtpstatLinearWeightedRegression(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, InputItW firstW, InputItW lastW,
double& coeffA,
double& coeffB,
bool fixA=
false,
bool fixB=
false, std::function<
double(
double)> fWeightDataToWi=&jkqtp_identity<double>) {
145 if (fixA&&fixB)
return;
146 const int Nx=std::distance(firstX,lastX);
147 const int Ny=std::distance(firstY,lastY);
148 const int Nw=std::distance(firstW,lastW);
152 double sumx=0, sumy=0, sumxy=0, sumx2=0, sumw2=0;
157 for (; itX!=lastX && itY!=lastY && itW!=lastW; ++itX, ++itY, ++itW) {
162 sumx=sumx+fit_w2*fit_x;
163 sumy=sumy+fit_w2*fit_y;
164 sumxy=sumxy+fit_w2*fit_x*fit_y;
165 sumx2=sumx2+fit_w2*fit_x*fit_x;
170 const double NN=
static_cast<double>(N);
172 if (!fixA && !fixB) {
173 coeffB=(double(sumxy)*double(sumw2)-double(sumx)*double(sumy))/(
double(sumx2)*double(sumw2)-double(sumx)*double(sumx));
174 coeffA=(double(sumy)-coeffB*double(sumx))/
double(sumw2);
175 }
else if (fixA && !fixB) {
176 coeffB=(double(sumy)-coeffA*double(sumw2))/
double(sumx);
177 }
else if (!fixA && fixB) {
178 coeffA=(double(sumy)-coeffB*double(sumx))/
double(sumw2);
228template <
class InputItX,
class InputItY>
229inline void jkqtpstatRobustIRLSLinearRegression(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY,
double& coeffA,
double& coeffB,
bool fixA=
false,
bool fixB=
false,
double p=1.1,
int iterations=100) {
230 if (fixA&&fixB)
return;
231 const int Nx=std::distance(firstX,lastX);
232 const int Ny=std::distance(firstY,lastY);
233 const int N=std::min(Nx,Ny);
237 std::vector<double> weights;
238 std::fill_n(std::back_inserter(weights), N, 1.0);
240 double alast=coeffA, blast=coeffB;
241 jkqtpstatLinearWeightedRegression(firstX, lastX, firstY, lastY, weights.begin(), weights.end(), alast, blast, fixA, fixB, &jkqtp_identity<double>);
242 for (
int it=0; it<iterations-1; it++) {
246 for (
double& w: weights) {
247 const double fit_x=*itX;
248 const double fit_y=*itY;
249 const double e=alast+blast*fit_x-fit_y;
250 w=pow(std::max<double>(
JKQTP_EPSILON*100.0, fabs(e)), (p-2.0)/2.0);
255 jkqtpstatLinearWeightedRegression(firstX, lastX, firstY, lastY, weights.begin(), weights.end(), alast, blast, fixA, fixB, &jkqtp_identity<double>);
337template <
class InputItX,
class InputItY>
339 std::vector<double> x, y;
344 std::transform(firstX, lastX, std::back_inserter(x), trafo.first);
345 std::transform(firstY, lastY, std::back_inserter(y), trafo.second);
347 double a=aTrafo.first(coeffA);
348 double b=bTrafo.first(coeffB);
352 coeffA=aTrafo.second(a);
353 coeffB=bTrafo.second(b);
383template <
class InputItX,
class InputItY>
384inline void jkqtpstatRobustIRLSRegression(
JKQTPStatRegressionModelType type, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY,
double& coeffA,
double& coeffB,
bool fixA=
false,
bool fixB=
false,
double p=1.1,
int iterations=100) {
385 std::vector<double> x, y;
390 std::transform(firstX, lastX, std::back_inserter(x), trafo.first);
391 std::transform(firstY, lastY, std::back_inserter(y), trafo.second);
393 double a=aTrafo.first(coeffA);
394 double b=bTrafo.first(coeffB);
398 coeffA=aTrafo.second(a);
399 coeffB=bTrafo.second(b);
435template <
class InputItX,
class InputItY,
class InputItW>
436inline void jkqtpstatWeightedRegression(
JKQTPStatRegressionModelType type, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, InputItW firstW, InputItW lastW,
double& coeffA,
double& coeffB,
bool fixA=
false,
bool fixB=
false, std::function<
double(
double)> fWeightDataToWi=&jkqtp_identity<double>) {
437 std::vector<double> x, y;
442 std::transform(firstX, lastX, std::back_inserter(x), trafo.first);
443 std::transform(firstY, lastY, std::back_inserter(y), trafo.second);
445 double a=aTrafo.first(coeffA);
446 double b=bTrafo.first(coeffB);
448 jkqtpstatLinearWeightedRegression(x.begin(), x.end(), y.begin(), y.end(), firstW, lastW, a, b, fixA, fixB, fWeightDataToWi);
450 coeffA=aTrafo.second(a);
451 coeffB=bTrafo.second(b);
473template <
class InputItX,
class InputItY>
482 for (; itX!=lastX && itY!=lastY; ++itX, ++itY) {
491 return 1.0-SSres/SStot;
519template <
class InputItX,
class InputItY,
class InputItW>
520inline double jkqtpstatWeightedCoefficientOfDetermination(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, InputItW firstW, InputItW lastW, std::function<
double(
double)> f, std::function<
double(
double)> fWeightDataToWi=&jkqtp_identity<double>) {
529 for (; itX!=lastX && itY!=lastY && itW!=lastW; ++itX, ++itY, ++itW) {
534 SSres+=(fit_w2*
jkqtp_sqr(fit_y-f(fit_x)));
539 return 1.0-SSres/SStot;
561template <
class InputItX,
class InputItY>
562inline double jkqtpstatSumOfDeviations(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, std::function<
double(
double)> f) {
568 for (; itX!=lastX && itY!=lastY; ++itX, ++itY) {
604template <
class InputItX,
class InputItY,
class InputItW>
605inline double jkqtpstatWeightedSumOfDeviations(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, InputItW firstW, InputItW lastW, std::function<
double(
double)> f, std::function<
double(
double)> fWeightDataToWi=&jkqtp_identity<double>) {
612 for (; itX!=lastX && itY!=lastY && itW!=lastW; ++itX, ++itY, ++itW) {
#define jkqtmath_LIB_EXPORT
Definition jkqtmath_imexport.h:87