JKQTPlotter trunk/v5.0.0
an extensive Qt5+Qt6 Plotter framework (including a feature-richt plotter widget, a speed-optimized, but limited variant and a LaTeX equation renderer!), written fully in C/C++ and without external dependencies
Loading...
Searching...
No Matches
jkqtpstatpoly.h
1/*
2 Copyright (c) 2008-2024 Jan W. Krieger (<jan@jkrieger.de>)
3
4 last modification: $LastChangedDate$ (revision $Rev$)
5
6 This software is free software: you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License (LGPL) as published by
8 the Free Software Foundation, either version 2.1 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU Lesser General Public License (LGPL) for more details.
15
16 You should have received a copy of the GNU Lesser General Public License (LGPL)
17 along with this program. If not, see <http://www.gnu.org/licenses/>.
18*/
19
20
21#ifndef JKQTPSTATPOLY_H_INCLUDED
22#define JKQTPSTATPOLY_H_INCLUDED
23
24#include <stdint.h>
25#include <cmath>
26#include <stdlib.h>
27#include <string.h>
28#include <iostream>
29#include <stdio.h>
30#include <limits>
31#include <vector>
32#include <utility>
33#include <cfloat>
34#include <ostream>
35#include <iomanip>
36#include <sstream>
37#include "jkqtmath/jkqtmath_imexport.h"
38#include "jkqtmath/jkqtplinalgtools.h"
39#include "jkqtmath/jkqtparraytools.h"
40#include "jkqtcommon/jkqtpdebuggingtools.h"
41#include <stdexcept>
42
43
44
45
46
47/*! \brief fits (in a least-squares sense) a polynomial \f$ f(x)=\sum\limits_{i=0}^Pp_ix^i \f$ of order P to a set of N data pairs \f$ (x_i,y_i) \f$
48 \ingroup jkqtptools_math_statistics_poly
49
50 \tparam InputItX standard iterator type of \a firstX and \a lastX.
51 \tparam InputItY standard iterator type of \a firstY and \a lastY.
52 \tparam OutputItP output iterator for the polynomial coefficients
53 \param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
54 \param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
55 \param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
56 \param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
57 \param P degree of the polynomial (P>=N !!!)
58 \param[out] firstRes Iterator (of type \a OutputItP ), which receives the (P+1)-entry vector with the polynomial coefficients \f$ p_i \f$
59
60 This function uses jkqtpstatLinSolve() to solve the system of equations
61 \f[ \begin{bmatrix} y_1\\ y_2\\ y_3 \\ \vdots \\ y_n \end{bmatrix}= \begin{bmatrix} 1 & x_1 & x_1^2 & \dots & x_1^P \\ 1 & x_2 & x_2^2 & \dots & x_2^P\\ 1 & x_3 & x_3^2 & \dots & x_3^P \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \dots & x_n^P \end{bmatrix} \begin{bmatrix} p_0\\ p_1\\ p_2\\ \vdots \\ p_P \end{bmatrix} \f]
62 \f[ \vec{y}=V\vec{p}\ \ \ \ \ \Rightarrow\ \ \ \ \ \vec{p}=(V^TV)^{-1}V^T\vec{y} \f]
63
64 \image html datastore_regression_polynom.png
65
66 \see https://en.wikipedia.org/wiki/Polynomial_regression
67*/
68template <class InputItX, class InputItY, class OutputItP>
69inline void jkqtpstatPolyFit(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, size_t P, OutputItP firstRes) {
70 {
71 const int Nx=std::distance(firstX,lastX);
72 const int Ny=std::distance(firstY,lastY);
73 JKQTPASSERT(Nx>1 && Ny>1);
74 }
75
76 size_t N=0;
77
78 std::vector<double> X,Y;
79 auto itX=firstX;
80 auto itY=firstY;
81 for (; itX!=lastX && itY!=lastY; ++itX, ++itY) {
82 const double fit_x=jkqtp_todouble(*itX);
83 const double fit_y=jkqtp_todouble(*itY);
84 if (JKQTPIsOKFloat(fit_x) && JKQTPIsOKFloat(fit_y)) {
85 X.push_back(fit_x);
86 Y.push_back(fit_y);
87 N++;
88 }
89 }
90
91 // build Vandermonde matrix V
92 std::vector<double> V;
93 V.resize(N*(P+1));
94 for (size_t l=0; l<N; l++) {
95 V[jkqtplinalgMatIndex(l,0,P+1)]=1.0;
96 double x=X[l];
97 const double xx=x;
98 for (size_t c=1; c<P+1; c++) {
99 V[jkqtplinalgMatIndex(l,c,P+1)]=x;
100 x=x*xx;
101 }
102 }
103#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
104 std::cout<<"V = \n";
105 jkqtplinalgPrintMatrix(V.data(),N,P+1);
106 std::cout<<"\n";
107#endif
108
109 // calculate V^T
110 std::vector<double> VT=V;
111 jkqtplinalgTransposeMatrix(VT.data(), static_cast<long>(N), static_cast<long>(P+1));
112
113#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
114 std::cout<<"V^T = \n";
115 jkqtplinalgPrintMatrix(VT.data(),P+1,N);
116 std::cout<<"\n";
117#endif
118
119 // calculate V^T*V
120 std::vector<double> VTV;
121 VTV.resize((P+1)*(P+1));
122 jkqtplinalgMatrixProduct(VT.data(), static_cast<long>(P+1), static_cast<long>(N), V.data(), static_cast<long>(N), static_cast<long>(P+1), VTV.data());
123
124#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
125 std::cout<<"V^T*V = \n";
126 jkqtplinalgPrintMatrix(VTV.data(),P+1,P+1);
127 std::cout<<"\n";
128#endif
129
130 // calculate V^T*y
131 std::vector<double> VTY;
132 VTY.resize(P+1);
133 jkqtplinalgMatrixProduct(VT.data(), static_cast<long>(P+1), static_cast<long>(N), Y.data(), static_cast<long>(N), 1, VTY.data());
134
135#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
136 std::cout<<"V^T*y = \n";
137 jkqtplinalgPrintMatrix(VTY.data(),P+1,1);
138 std::cout<<"\n";
139#endif
140
141 // solve V^T*y = V^T*V*p
142 const bool ok=jkqtplinalgLinSolve(VTV.data(), VTY.data(), static_cast<long>(P+1));
143
144 if (ok) {
145 auto itR=firstRes;
146 for (size_t p=0; p<P+1; p++) {
147 *itR=VTY[p];
148 ++itR;
149 }
150 } else {
151 throw std::runtime_error("jkqtplinalgLinSolve() didn't return a result!");
152 }
153
154#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
155 std::cout<<"result_out = \n";
156 jkqtplinalgPrintMatrix(result_out,P+1,1);
157 std::cout<<"\n";
158#endif
159
160}
161
162
163
164
165
166#endif // JKQTPSTATPOLY_H_INCLUDED
167
168
#define JKQTPASSERT(condition)
dynamic assertion, throws an exception with the given message, when the given condition condition eva...
Definition jkqtpdebuggingtools.h:77
constexpr double jkqtp_todouble(const T &d)
converts a boolean to a double, is used to convert boolean to double by JKQTPDatastore
Definition jkqtpmathtools.h:113
bool JKQTPIsOKFloat(T v)
check whether the dlotaing point number is OK (i.e. non-inf, non-NAN)
Definition jkqtpmathtools.h:496
#define jkqtplinalgMatIndex(l, c, C)
calculate the index of the entry in line l and column c in a row-major matrix with C columns
Definition jkqtplinalgtools.h:79
void jkqtplinalgTransposeMatrix(T *matrix, long N)
transpose the given NxN matrix
Definition jkqtplinalgtools.h:319
void jkqtplinalgPrintMatrix(const T *matrix, long L, long C, int width=9, int precision=3, char mode='f')
print the given LxC matrix to std::cout
Definition jkqtplinalgtools.h:94
bool jkqtplinalgLinSolve(const T *A, const T *B, long N, long C, T *result_out)
solve a system of N linear equations simultaneously for C columns in B
Definition jkqtplinalgtools.h:631
void jkqtplinalgMatrixProduct(const T *M1, long L1, long C1, const T *M2, long L2, long C2, T *M)
matrix-matrix product
Definition jkqtplinalgtools.h:378
void jkqtpstatPolyFit(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, size_t P, OutputItP firstRes)
fits (in a least-squares sense) a polynomial of order P to a set of N data pairs
Definition jkqtpstatpoly.h:69