added regression, IRLS robust regression, weighted regression and polynomial fitting/evaluation to statistics library

added example for regression, IRLS robust regression, weighted regression and polynomial fitting
This commit is contained in:
jkriege2 2019-06-01 18:15:04 +02:00
parent 199e8493fe
commit 2d0b1e7935
25 changed files with 2127 additions and 15 deletions

View File

@ -82,6 +82,7 @@ addSimpleTest(imageplot_nodatastore)
addSimpleTest(datastore)
addSimpleTest(datastore_iterators)
addSimpleTest(datastore_statistics)
addSimpleTest(datastore_regression)
addSimpleTest(contourplot)
#addSimpleTest(rgbimageplot_opencv)
#addSimpleTest(imageplot_opencv)

View File

@ -159,6 +159,9 @@ All test-projects are Qt-projects that use qmake to build. You can load them int
<tr><td> \image html jkqtplotter_simpletest_datastore_statistics_small.png
<td> \subpage JKQTPlotterBasicJKQTPDatastoreStatistics
<td> Advanced 1-Dimensional Statistical Computation with JKQTPDatastore<br>using the internal statistics library (see \ref jkqtptools_math_statistics )<br>basic statistics (mean, standard deviation, ...)<br>boxplots<br>histograms<br>kernel density estimates (KDE)
<tr><td> \image html jkqtplotter_simpletest_datastore_regression_small.png
<td> \subpage JKQTPlotterBasicJKQTPDatastoreRegression
<td> Advanced 1-Dimensional Statistical Computation with JKQTPDatastore<br>using the internal statistics library<br>Regression Analysis (with the Statistics Library)<br>robust regression (IRLS)<br>weighted regression<br>non-linear regression<br>polynomial fitting
</table>

View File

@ -65,6 +65,7 @@ All test-projects are Qt-projects that use qmake to build. You can load them int
| [![](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/simpletest_datastore_small.png)](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore) | [Tutorial: Basic Usage of JKQTPDatastore](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore) | Basic Data Management with JKQTPDatastore <br/> Copying data into a JKQTPDatastore <br/> Editing data inside a JKQTPDatastore <br/> Editing Image Data in a JKQTPDatastore |
| [![](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/simpletest_datastore_iterators_small.png)](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_iterators) | [Tutorial: Iterator-based access to JKQTPDatastore](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_iterators) | Iterator-based Data Management with JKQTPDatastore |
| [![](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_statistics_small.png)](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_statistics) | [Tutorial: Advanced 1-Dimensional Statistics with JKQTPDatastore](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_statistics) | Advanced 1-Dimensional Statistical Computation with JKQTPDatastore<br>using the internal statistics library<br>basic statistics (mean, standard deviation, ...)<br>boxplots<br>histograms<br>kernel density estimates (KDE) |
| [![](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_small.png)](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_regression) | [Tutorial: Regression Analysis (with the Statistics Library)](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_regression) | Advanced 1-Dimensional Statistical Computation with JKQTPDatastore<br>using the internal statistics library<br>Regression Analysis (with the Statistics Library)<br>robust regression (IRLS)<br>weighted regression<br>non-linear regression<br>polynomial fitting |
## More Complex Examples

View File

@ -1,8 +1,10 @@
# Tutorial (JKQTPDatastore): Basic Usage of JKQTPDatastore {#JKQTPlotterBasicJKQTPDatastore}
[JKQTPlotterBasicJKQTPDatastore]: @ref JKQTPlotterBasicJKQTPDatastore "Basic Usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreIterators]: @ref JKQTPlotterBasicJKQTPDatastoreIterators "Iterator-Based usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreStatistics]: @ref JKQTPlotterBasicJKQTPDatastoreStatistics "Advanced 1-Dimensional Statistics with JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreRegression]: @ref JKQTPlotterBasicJKQTPDatastoreRegression "Regression Analysis (with the Statistics Library)"
[statisticslibrary]: @ref jkqtptools_math_statistics "JKQTPlotter Statistics Library"
This tutorial project (see `./examples/simpletest_datastore/`) explains several options of JKQTPDatastore, which is the class used to centrally store the data for (most) graphs on a JKQTPlotter widget.
@ -11,6 +13,8 @@ This tutorial project (see `./examples/simpletest_datastore/`) explains several
- [JKQTPlotterBasicJKQTPDatastore]
- [JKQTPlotterBasicJKQTPDatastoreIterators]
- [JKQTPlotterBasicJKQTPDatastoreStatistics]
- [JKQTPlotterBasicJKQTPDatastoreRegression]
[TOC]

View File

@ -1,8 +1,10 @@
# Tutorial (JKQTPDatastore): Iterator-Based usage of JKQTPDatastore {#JKQTPlotterBasicJKQTPDatastoreIterators}
[JKQTPlotterBasicJKQTPDatastore]: @ref JKQTPlotterBasicJKQTPDatastore "Basic Usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreIterators]: @ref JKQTPlotterBasicJKQTPDatastoreIterators "Iterator-Based usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreStatistics]: @ref JKQTPlotterBasicJKQTPDatastoreStatistics "Advanced 1-Dimensional Statistics with JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreRegression]: @ref JKQTPlotterBasicJKQTPDatastoreRegression "Regression Analysis (with the Statistics Library)"
[statisticslibrary]: @ref jkqtptools_math_statistics "JKQTPlotter Statistics Library"
This tutorial project (see `./examples/simpletest_datastore_iterators/`) explains how to use the iterator-based interface to JKQTPDatastore.
@ -11,6 +13,8 @@ This tutorial project (see `./examples/simpletest_datastore_iterators/`) explain
- [JKQTPlotterBasicJKQTPDatastore]
- [JKQTPlotterBasicJKQTPDatastoreIterators]
- [JKQTPlotterBasicJKQTPDatastoreStatistics]
- [JKQTPlotterBasicJKQTPDatastoreRegression]
[TOC]

View File

@ -0,0 +1,362 @@
# Tutorial (JKQTPDatastore): Regression Analysis (with the Statistics Library) {#JKQTPlotterBasicJKQTPDatastoreRegression}
[JKQTPlotterBasicJKQTPDatastore]: @ref JKQTPlotterBasicJKQTPDatastore "Basic Usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreIterators]: @ref JKQTPlotterBasicJKQTPDatastoreIterators "Iterator-Based usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreStatistics]: @ref JKQTPlotterBasicJKQTPDatastoreStatistics "Advanced 1-Dimensional Statistics with JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreRegression]: @ref JKQTPlotterBasicJKQTPDatastoreRegression "Regression Analysis (with the Statistics Library)"
[statisticslibrary]: @ref jkqtptools_math_statistics "JKQTPlotter Statistics Library"
This tutorial project (see `./examples/simpletest_datastore_statistics/`) explains several advanced functions of JKQTPDatastore in combination with the [[statisticslibrary]] conatined in JKQTPlotter.
***Note*** that there are additional tutorial explaining other aspects of data mangement in JKQTPDatastore:
- [JKQTPlotterBasicJKQTPDatastore]
- [JKQTPlotterBasicJKQTPDatastoreIterators]
- [JKQTPlotterBasicJKQTPDatastoreStatistics]
- [JKQTPlotterBasicJKQTPDatastoreRegression]
[TOC]
The source code of the main application can be found in [`jkqtplotter_simpletest_datastore_regression.cpp`](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_statistics/jkqtplotter_simpletest_datastore_regression.cpp).
This tutorial cites only parts of this code to demonstrate different ways of performing regression analysis.
# Simple Linear Regression
First we generate a set of datapoints (x,y), which scatter randomly around a linear function.
```.cpp
std::random_device rd; // random number generators:
std::mt19937 gen{rd()};
std::normal_distribution<> d1{0,1};
double a0=-5;
double b0=2;
size_t colLinX=datastore1->addColumn("lin data, x");
size_t colLinY=datastore1->addColumn("lin data, y");
for (double x=-5; x<=10; x++) {
datastore1->appendToColumn(colLinX, x);
datastore1->appendToColumn(colLinY, a0+b0*x+d1(gen));
}
```
... and we visualize this data with a simple scatter graph:
```.cpp
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
graphD->setDrawLine(false);
graphD->setTitle("data");
```
Now we can caluate the regression line (i.e. the two regression coefficients a and b of the function \c f(x)=a+b*x) using the function `jkqtpstatLinearRegression()` from the [statisticslibrary]:
```.cpp
double coeffA=0, coeffB=0;
jkqtpstatLinearRegression(datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY), coeffA, coeffB);
```
... and add a `JKQTPXFunctionLineGraph` to draw the resulting linear function:
```.cpp
JKQTPXFunctionLineGraph *graphRegLine=new JKQTPXFunctionLineGraph(plot1);
graphRegLine->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
graphRegLine->setParamsV(coeffA, coeffB);
graphRegLine->setTitle(QString("regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(coeffA)).arg(jkqtp_floattolatexqstr(coeffB)));
plot1->addGraph(graphRegLine);
```
These two steps can be simplified using an "adaptor":
```.cpp
jkqtpstatAddLinearRegression(plot1->getPlotter(), datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY));
```
... or even shorter:
```.cpp
jkqtpstatAddLinearRegression(graphD);
```
Here the x- and y-columns from the `JKQTPXYGraph`-based graph `graphD` (see above) are used as datasources for the plot.
The plot resulting from any of the variants above looks like this:
![jkqtplotter_simpletest_datastore_regression_lin](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_lin.png)
# Robust Linear Regression
Sometimes data contains outliers that can render the results of a regression analysis inaccurate. For such cases the [statisticslibrary] offers the function `jkqtpstatRobustIRLSLinearRegression()`, which is a drop-in replacement for `jkqtpstatLinearRegression()` and solves the optimization problem a) in the Lp-norm (which is more robust to outliers) and b) uses the [iteratively reweighted least-squares algorithm (IRLS)](https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares), which performs a series of regressions, where in each instance the data-points are weighted differently. The method assigns a lower weight to those points that are far from the current best-fit (typically the outliers) and thus slowly comes nearer to an estimate that is not distorted by the outliers.
To demonstrate this method, we use the same dataset as above, but add a few outliers:
```.cpp
std::random_device rd; // random number generators:
std::mt19937 gen{rd()};
std::normal_distribution<> d1{0,1};
double a0=-5;
double b0=2;
size_t colLinX=datastore1->addColumn("lin data, x");
size_t colLinY=datastore1->addColumn("lin data, y");
for (double x=-5; x<=10; x++) {
datastore1->appendToColumn(colLinX, x);
if (jkqtp_approximatelyEqual(x, -5)||jkqtp_approximatelyEqual(x, -3)) datastore1->appendToColumn(colRLinY, a0+b0*x+d1(gen)+12);
else datastore1->appendToColumn(colRLinY, a0+b0*x+d1(gen));
}
```
Note the outliers ar x=-5 and x=-3!
With this dataset we can use the same code as above, but with `jkqtpstatRobustIRLSLinearRegression()` instead of `jkqtpstatLinearRegression()`:
```.cpp
double coeffA=0, coeffB=0;
jkqtpstatRobustIRLSLinearRegression(datastore1->begin(colRLinX), datastore1->end(colRLinX), datastore1->begin(colRLinY), datastore1->end(colRLinY), coeffA, coeffB);
JKQTPXFunctionLineGraph *graphRegLine=new JKQTPXFunctionLineGraph(plot3);
graphRegLine->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
graphRegLine->setParamsV(coeffA, coeffB);
graphRegLine->setTitle(QString("robust regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(coeffA)).arg(jkqtp_floattolatexqstr(coeffB)));
plot3->addGraph(graphRegLine);
```
Also for the robust regression, there are two shortcuts in the form of "adaptors":
```.cpp
jkqtpstatAddRobustIRLSLinearRegression(plot3->getPlotter(), datastore1->begin(colRLinX), datastore1->end(colRLinX), datastore1->begin(colRLinY), datastore1->end(colRLinY));
```
and
```.cpp
jkqtpstatAddRobustIRLSLinearRegression(graphD);
```
The following screenshot shows the result of the IRLS regression analysis and for comparison the normal linear regression for the same dataset (plotted using `jkqtpstatAddLinearRegression(graphD);`):
![jkqtplotter_simpletest_datastore_regression_linrobust](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_linrobust.png)
The following screenshot shows the influence of the regularization parameter p (default value 1.1) onto the fit result:
- the closer `p` is to 1, the more robust the fit is (it is closer to the L1-norm)
- the closer `p` is to 2, the closer the fit is to the least squares solution (i.e. the normal regression obtained with the L2 norm)
![jkqtplotter_simpletest_datastore_regression_linrobust_p](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_linrobust_p.png)
# Weighted Linear Regression
Another option to react to measurement errors is to take these into account when calculating the regression. To do so, you can use weighted regression that uses the measurement errors as inverse weights. This algorithm is implemented in the function `jkqtpstatLinearWeightedRegression()`.
First we generate again a set of datapoints (x,y), which scatter randomly around a linear function. In addition we calculate an "error" `err` for each datapoint:
```.cpp
std::random_device rd; // random number generators:
std::mt19937 gen{rd()};
std::uniform_real_distribution<> de{0.5,1.5};
std::uniform_int_distribution<> ddecide{0,4};
double a0=-5;
double b0=2;
size_t colWLinX=datastore1->addColumn("wlin data, x");
size_t colWLinY=datastore1->addColumn("wlin data, y");
size_t colWLinE=datastore1->addColumn("wlin data, errors");
for (double x=-5; x<=10; x++) {
double factor=1;
if (ddecide(gen)==4) {
factor=4;
}
const double err=de(gen)*factor;
datastore1->appendToColumn(colWLinX, x);
datastore1->appendToColumn(colWLinY, a0+b0*x+err);
datastore1->appendToColumn(colWLinE, 1.0/err);
}
```
We use distribution `de` to draw deviations from the ideal linear function from the range 0.5...1.5. then - for good measure - we use a second distribution `ddecide` (dice tossing) to select a few datapoints to have a 4-fold increased error.
Finally we visualize this data with a simple scatter graph with error indicators:
```.cpp
JKQTPXYLineErrorGraph* graphE;
plot2->addGraph(graphE=new JKQTPXYLineErrorGraph(plot2));
graphE->setXYColumns(colWLinX, colWLinY);
graphE->setYErrorColumn(static_cast<int>(colWLinE));
graphE->setDrawLine(false);
graphE->setTitle("data");
```
Now we can caluate the regression line (i.e. the two regression coefficients a and b of the function \c f(x)=a+b*x) using the function `jkqtpstatLinearWeightedRegression()` from the [statisticslibrary]:
```.cpp
double coeffA=0, coeffB=0;
jkqtpstatLinearWeightedRegression(datastore1->begin(colWLinX), datastore1->end(colWLinX),
datastore1->begin(colWLinY), datastore1->end(colWLinY),
datastore1->begin(colWLinE), datastore1->end(colWLinE),
coeffA, coeffB, false, false,
&jkqtp_inversePropSaveDefault<double>);
```
***Note*** that in addition to the three data-columns we also provided a C++ functor `jkqtp_inversePropSaveDefault()`, which calculates 1/error. This is done, because the function `jkqtpstatLinearWeightedRegression()` uses the data from the range `datastore1->begin(colWLinE)` ... `datastore1->end(colWLinE)` directly as weights, but we calculated errors, which are inversely proportional to the weight of each data point when solving the least squares problem, as data points with larger errors should be weighted less than thos with smaller errors (outliers).
Again these two steps can be simplified using an "adaptor":
```.cpp
jkqtpstatAddLinearWeightedRegression(plot1->getPlotter(),
datastore1->begin(colLinX), datastore1->end(colLinX),
datastore1->begin(colLinY), datastore1->end(colLinY),
datastore1->begin(colWLinE), datastore1->end(colWLinE),
&coeffA, &coeffB, false, false,
&jkqtp_inversePropSaveDefault<double>);
```... or even shorter:
```.cpp
jkqtpstatAddLinearWeightedRegression(graphD);
```
Here the x- and y-columns from the `JKQTPXYGraph`-based graph `graphE` (see above) and the weights from the error column of `graphE` are used as datasources for the plot. This function implicitly uses the function `jkqtp_inversePropSaveDefault()` to convert plot errors to weights, as it is already clear that we are dealing with errors rather than direct weights.
The plot resulting from any of the variants above looks like this:
![jkqtplotter_simpletest_datastore_regression_linweight](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_linweight.png)
For this plot we also added a call to
```.cpp
jkqtpstatAddLinearRegression(graphE);
```
which performs a simple non-weighted regression. The difference between the two resulting linear functions (blue: simple regression, green: weighted regression) demonstrates the influence of the weighting.
# Linearizable Regression Models
In addition to the simple linear regression model `f(x)=a+b*x`, it is also possible to fit a few non-linear models by transforming the data:
- power-law function`f(x)=a*x^b`, which is a linear function in a log-log-plot
- exponential function `f(x)=a*exp(b*x)`, which is a linear function in a semi-log-plot
The available models are defined in the enum `JKQTPStatRegressionModelType`. And there exists a function `jkqtpStatGenerateRegressionModel()`, which returns a C++-functor representing the function.
To demonstrate these fitting options, we first generate data from an exponential and a power-law model. Note that we also add normally distributed errors, but in order to ensure that we do not obtain y-values <0, we use loops that draw normally distributed random numbers, until this condition is met:
```.cpp
std::random_device rd; // random number generators:
std::mt19937 gen{rd()};
std::normal_distribution<> d1{0,1};
double a0_powerlaw=20;
double b0_powerlaw=0.25;
double a0_exp=5;
double b0_exp=0.5;
size_t colNLLinX=datastore1->addColumn("non-lin data, x");
size_t colNLLinYExp=datastore1->addColumn("non-lin data, y, exponential model");
size_t colNLLinYPow=datastore1->addColumn("non-lin data, y, power-law model");
auto model_powerlaw=jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::PowerLaw);
auto model_exp=jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::Exponential);
for (double x=0.1; x<=10; x+=0.5) {
datastore1->appendToColumn(colNLLinX, x);
double ypow=model_powerlaw(x, a0_powerlaw, b0_powerlaw)+d1(gen);
while (ypow<0) {
ypow=model_powerlaw(x, a0_powerlaw, b0_powerlaw)+d1(gen);
}
datastore1->appendToColumn(colNLLinYPow, ypow);
double yexp=model_exp(x, a0_powerlaw, b0_powerlaw)+d1(gen);
while (yexp<0) {
yexp=model_exp(x, a0_powerlaw, b0_powerlaw)+d1(gen);
}
datastore1->appendToColumn(colNLLinYExp, model_exp(x, a0_exp, b0_exp)+d1(gen));
}
```
The generated data is visualized with scatter-plots:
```.cpp
JKQTPXYLineGraph* graphD_powerlaw;
plot5->addGraph(graphD_powerlaw=new JKQTPXYLineGraph(plot5));
graphD_powerlaw->setXYColumns(colNLLinX, colNLLinYPow);
graphD_powerlaw->setDrawLine(false);
graphD_powerlaw->setTitle(QString("data $%1+\\mathcal{N}(0,1)$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::PowerLaw, a0_powerlaw, b0_powerlaw)));
JKQTPXYLineGraph* graphD_exp;
plot4->addGraph(graphD_exp=new JKQTPXYLineGraph(plot4));
graphD_exp->setXYColumns(colNLLinX, colNLLinYExp);
graphD_exp->setDrawLine(false);
graphD_exp->setTitle(QString("data $%1+\\mathcal{N}(0,1)$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::Exponential, a0_exp, b0_exp)));
```
Now we can fit the regression models using `jkqtpstatRegression()`, which receives the model type as first parameter:
```.cpp
double cA=0, cB=0;
JKQTPXFunctionLineGraph* gFunc;
jkqtpstatRegression(JKQTPStatRegressionModelType::Exponential, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYExp), datastore1->end(colNLLinYExp), cA, cB);
plot4->addGraph(gFunc=new JKQTPXFunctionLineGraph(plot4));
gFunc->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::Exponential, cA, cB));
gFunc->setTitle(QString("regression: $%1$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::Exponential, cA, cB)));
cA=0; cB=0;
jkqtpstatRegression(JKQTPStatRegressionModelType::PowerLaw, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYPow), datastore1->end(colNLLinYPow), cA, cB);
plot5->addGraph(gFunc=new JKQTPXFunctionLineGraph(plot5));
gFunc->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::PowerLaw, cA, cB));
gFunc->setTitle(QString("regression: $%1$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::PowerLaw, cA, cB)));
```
The regression models can be plotted using a `JKQTPXFunctionLineGraph`. the fucntion to plot is again generated by calling `jkqtpStatGenerateRegressionModel()`, but now with the parameters determined above the respective lines. Note that `jkqtpstatRegressionModel2Latex()` outputs the model as LaTeX string, which can be used as plot label.
The resulting plot looks like this:
![jkqtplotter_simpletest_datastore_regression_nonlinreg](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_nonlinreg.png)
Of course also "adaptors" exist that allow to perform the steps above in a single function call:
```.cpp
jkqtpstatAddRegression(plot4->getPlotter(), JKQTPStatRegressionModelType::Exponential, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYExp), datastore1->end(colNLLinYExp));
jkqtpstatAddRegression(plot5->getPlotter(), JKQTPStatRegressionModelType::PowerLaw, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYPow), datastore1->end(colNLLinYPow));
```
... or even shorter:
```.cpp
jkqtpstatAddRegression(graphD_exp, JKQTPStatRegressionModelType::Exponential);
jkqtpstatAddRegression(graphD_powerlaw, JKQTPStatRegressionModelType::PowerLaw);
```
Also note that we used the function `jkqtpstatRegression()` above, which performs a linear regression (internally uses `jkqtpstatLinearRegression()`). But there also exist variants for robust IRLS regression adn weighted regression:
- `jkqtpstatRobustIRLSRegression()` / `jkqtpstatAddRobustIRLSRegression()`
- `jkqtpstatWeightedRegression()` / `jkqtpstatAddWeightedRegression()`
# Polynom Fitting
Finally the [statisticslibrary] also supports one option for non-linear model fitting, namely fitting of polynomial models. This is implemented in the function `jkqtpstatPolyFit()`.
To demonstrate this function we first generate data from a poylnomial model (with gaussian noise):
```.cpp
std::random_device rd; // random number generators:
std::mt19937 gen{rd()};
std::normal_distribution<> d1{0,50};
std::vector<double> pPoly {1,2,-2,0.5};
size_t colPolyX=datastore1->addColumn("polynomial data, x");
size_t colPolyY=datastore1->addColumn("polynomial data, y");
for (double x=-10; x<=10; x++) {
datastore1->appendToColumn(colPolyX, x);
datastore1->appendToColumn(colPolyY, jkqtpstatPolyEval(x, pPoly.begin(), pPoly.end())+d1(gen));
}
```
The function `jkqtpstatPolyEval()` is used to evaluate a given polynomial (coefficients in `pPoly`) at a position `x`.
The generated data is visualized with scatter-plots:
```.cpp
JKQTPXYLineGraph* graphP;
plot6->addGraph(graphP=new JKQTPXYLineGraph(plot6));
graphP->setXYColumns(colPolyX, colPolyY);
graphP->setDrawLine(false);
graphP->setTitle(QString("data $%1+\\mathcal{N}(0,50)$").arg(jkqtpstatPolynomialModel2Latex(pPoly.begin(), pPoly.end())));
```
Here the function `jkqtpstatPolynomialModel2Latex()` generates a string from a polynomial model.
Now we can call `jkqtpstatPolyFit()` to fit different polynomial regression models to the data:
```.cpp
for (size_t p=0; p<=5; p++) {
std::vector<double> pFit;
JKQTPXFunctionLineGraph* gPoly;
jkqtpstatPolyFit(datastore1->begin(colPolyX), datastore1->end(colPolyX), datastore1->begin(colPolyY), datastore1->end(colPolyY), p, std::back_inserter(pFit));
plot6->addGraph(gPoly=new JKQTPXFunctionLineGraph(plot6));
gPoly->setPlotFunctionFunctor(jkqtpstatGeneratePolynomialModel(pFit.begin(), pFit.end()));
gPoly->setTitle(QString("regression: $%1$").arg(jkqtpstatPolynomialModel2Latex(pFit.begin(), pFit.end())));
}
```
Each model is also ploted using a `JKQTPXFunctionLineGraph`. The plot function assigned to these `JKQTPXFunctionLineGraph` is generated by calling `jkqtpstatGeneratePolynomialModel()`, which returns a C++-functor for a polynomial.
The resulting plots look like this (without added gaussian noise):
![jkqtplotter_simpletest_datastore_regression_polynom](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_polynom.png)
... and with added gaussian noise:
![jkqtplotter_simpletest_datastore_regression_polynom_errros](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression_polynom_errros.png)
Of course also the "adaptor" shortcuts are available:
```.cpp
for (size_t p=0; p<=5; p++) {
jkqtpstatAddPolyFit(plot6->getPlotter(), datastore1->begin(colPolyX), datastore1->end(colPolyX), datastore1->begin(colPolyY), datastore1->end(colPolyY), p);
jkqtpstatAddPolyFit(graphP, p);
}
```
# Screenshot of the full Program
The output of the full test program [`jkqtplotter_simpletest_datastore_regression.cpp`](https://github.com/jkriege2/JKQtPlotter/tree/master/examples/simpletest_datastore_statistics/jkqtplotter_simpletest_datastore_regression.cpp) looks like this:
![jkqtplotter_simpletest_datastore_regression](https://raw.githubusercontent.com/jkriege2/JKQtPlotter/master/screenshots/jkqtplotter_simpletest_datastore_regression.png)

View File

@ -0,0 +1,327 @@
/** \example jkqtplotter_simpletest_datastore_regression.cpp
* Explains how to use the internal statistics library (see \ref jkqtptools_statistics ) together with JKQTPDatastore to perform different types of regression and polynomial fitting.
*
* \ref JKQTPlotterBasicJKQTPDatastoreRegression
*/
#include <QApplication>
#include "jkqtplotter/jkqtplotter.h"
#include "jkqtplotter/jkqtpgraphspeakstream.h"
#include "jkqtplotter/jkqtpgraphsboxplot.h"
#include "jkqtplotter/jkqtpgraphsstatisticsadaptors.h"
#include "jkqtplotter/jkqtpgraphsevaluatedfunction.h"
#include "jkqtcommon/jkqtpstatisticstools.h"
#include "jkqtcommon/jkqtpstringtools.h"
#include <random>
#include <cmath>
int main(int argc, char* argv[])
{
QApplication app(argc, argv);
// 1. create a window with several plotters and get a pointer to the internal datastores (for convenience)
QWidget mainWidget;
QGridLayout* lay;
mainWidget.setLayout(lay=new QGridLayout);
JKQTPlotter* plot1=new JKQTPlotter(&mainWidget);
plot1->getPlotter()->setPlotLabel("Simple Linear Regression");
JKQTPDatastore* datastore1=plot1->getDatastore();
lay->addWidget(plot1,0,0);
JKQTPlotter *plot2=new JKQTPlotter(datastore1, &mainWidget);
plot2->getPlotter()->setPlotLabel("Weighted Linear Regression");
lay->addWidget(plot2,1,0);
JKQTPlotter* plot3=new JKQTPlotter(datastore1, &mainWidget);
plot3->getPlotter()->setPlotLabel("Robust Linear Regression");
lay->addWidget(plot3,0,1);
JKQTPlotter *plot6=new JKQTPlotter(datastore1, &mainWidget);
plot6->getPlotter()->setPlotLabel("Polynomial Fitting");
lay->addWidget(plot6,1,1);
JKQTPlotter* plot4=new JKQTPlotter(datastore1, &mainWidget);
plot4->getPlotter()->setPlotLabel("Exponential Regression");
lay->addWidget(plot4,0,2);
JKQTPlotter* plot5=new JKQTPlotter(datastore1, &mainWidget);
plot5->getPlotter()->setPlotLabel("Power-Law Regression");
lay->addWidget(plot5,1,2);
// 2.1. To demonstrate linear regression, we create a dataset with a linear dependence between two
// columns and added gaussian noise
std::random_device rd; // random number generators:
std::mt19937 gen{rd()};
std::normal_distribution<> d1{0,1};
double a0=-5;
double b0=2;
size_t colLinX=datastore1->addColumn("lin data, x");
size_t colLinY=datastore1->addColumn("lin data, y");
for (double x=-5; x<=10; x++) {
datastore1->appendToColumn(colLinX, x);
datastore1->appendToColumn(colLinY, a0+b0*x+d1(gen));
}
// we visualize this data with a simple scatter graph:
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
graphD->setDrawLine(false);
graphD->setTitle(QString("data $f(x)=%1+%2\\cdot x+\\mathcal{N}(0,1)$").arg(jkqtp_floattolatexqstr(a0,1)).arg(jkqtp_floattolatexqstr(b0,1)));
// 2.2. Now we calculate the regression line and add a plot to the graph:
/*double coeffA=0, coeffB=0;
jkqtpstatLinearRegression(datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY), coeffA, coeffB);
JKQTPXFunctionLineGraph *graphRegLine=new JKQTPXFunctionLineGraph(plot1);
graphRegLine->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
graphRegLine->setParamsV(coeffA, coeffB);
graphRegLine->setTitle(QString("regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(coeffA)).arg(jkqtp_floattolatexqstr(coeffB)));
plot1->addGraph(graphRegLine);*/
// this code can also be written with one function call, using the "adaptor" jkqtpstatAddLinearRegression():
//jkqtpstatAddLinearRegression(plot1->getPlotter(), datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY));
// or even shorter:
jkqtpstatAddLinearRegression(graphD);
// 3.1. We extend the example above by
//std::random_device rd; // random number generators:
//std::mt19937 gen{rd()};
std::uniform_real_distribution<> de{0.5,1.5};
std::uniform_int_distribution<> ddecide{0,4};
//double a0=-5;
//double b0=2;
size_t colWLinX=datastore1->addColumn("wlin data, x");
size_t colWLinY=datastore1->addColumn("wlin data, y");
size_t colWLinE=datastore1->addColumn("wlin data, errors");
for (double x=-5; x<=10; x++) {
double factor=1;
if (ddecide(gen)==4) {
factor=4;
}
const double err=de(gen)*factor;
datastore1->appendToColumn(colWLinX, x);
datastore1->appendToColumn(colWLinY, a0+b0*x+err);
datastore1->appendToColumn(colWLinE, err);
}
// we visualize this data with a simple scatter graph:
JKQTPXYLineErrorGraph* graphE;
plot2->addGraph(graphE=new JKQTPXYLineErrorGraph(plot2));
graphE->setXYColumns(colWLinX, colWLinY);
graphE->setYErrorColumn(static_cast<int>(colWLinE));
graphE->setDrawLine(false);
graphE->setTitle(QString("data $f(x)=%1+%2\\cdot x+\\mbox{Noise}$").arg(jkqtp_floattolatexqstr(a0,1)).arg(jkqtp_floattolatexqstr(b0,1)));
// 2.2. Now we calculate the regression line and add a plot to the graph:
/*double coeffA=0, coeffB=0;
jkqtpstatLinearWeightedRegression(datastore1->begin(colWLinX), datastore1->end(colWLinX),
datastore1->begin(colWLinY), datastore1->end(colWLinY),
datastore1->begin(colWLinE), datastore1->end(colWLinE),
coeffA, coeffB, false, false,
&jkqtp_inversePropSaveDefault<double>);
// note that in addition to the three data-columns we also provided a C++ functor
// jkqtp_inversePropSaveDefault(), which calculates 1/error. This is done, because the function
// jkqtpstatLinearWeightedRegression() uses the data from the range datastore1->begin(colWLinE) ... datastore1->end(colWLinE)
// directly as weights, but we calculated errors, which are inversely proportional to the
// weight of each data point when solving the least squares problem, as data points with
// larger errors should be weighted less than thos with smaller errors
//
// Now we can plot the resulting linear function:
JKQTPXFunctionLineGraph *graphRegLine=new JKQTPXFunctionLineGraph(plot2);
graphRegLine->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
graphRegLine->setParamsV(coeffA, coeffB);
graphRegLine->setTitle(QString("weighted regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(coeffA)).arg(jkqtp_floattolatexqstr(coeffB)));
plot2->addGraph(graphRegLine);*/
// this code can also be written with one function call, using the "adaptor" jkqtpstatAddLinearRegression():
//jkqtpstatAddLinearWeightedRegression(plot2->getPlotter(),
// datastore1->begin(colLinX), datastore1->end(colLinX),
// datastore1->begin(colLinY), datastore1->end(colLinY),
// datastore1->begin(colWLinE), datastore1->end(colWLinE),
// &coeffA, &coeffB, false, false,
// &jkqtp_inversePropSaveDefault<double>);
// or even shorter:
jkqtpstatAddLinearWeightedRegression(graphE);
// to demonstrate the effect of the weighting, we also add a simple linear regression that
// does not take into account the errors:
jkqtpstatAddLinearRegression(graphE);
// 4.1. To demonstrate IRLS linear regression, we create a dataset with a linear dependence between two
// columns and added gaussian noise and some outliers
//std::random_device rd; // random number generators:
//std::mt19937 gen{rd()};
//std::normal_distribution<> d1{0,1};
//double a0=-5;
//double b0=2;
size_t colRLinX=datastore1->addColumn("lin data, x");
size_t colRLinY=datastore1->addColumn("lin data, y");
for (double x=-5; x<=10; x++) {
datastore1->appendToColumn(colRLinX, x);
if (jkqtp_approximatelyEqual(x, -5)||jkqtp_approximatelyEqual(x, -3)) datastore1->appendToColumn(colRLinY, a0+b0*x+d1(gen)+12);
else datastore1->appendToColumn(colRLinY, a0+b0*x+d1(gen));
}
// we visualize this data with a simple scatter graph:
//JKQTPXYLineGraph* graphD;
plot3->addGraph(graphD=new JKQTPXYLineGraph(plot3));
graphD->setXYColumns(colRLinX, colRLinY);
graphD->setDrawLine(false);
graphD->setTitle(QString("data $f(x)=%1+%2\\cdot x+\\mathcal{N}(0,1)$").arg(jkqtp_floattolatexqstr(a0,1)).arg(jkqtp_floattolatexqstr(b0,1)));
// 4.2. Now we calculate the regression line and add a plot to the graph:
double coeffA=0, coeffB=0;
jkqtpstatRobustIRLSLinearRegression(datastore1->begin(colRLinX), datastore1->end(colRLinX), datastore1->begin(colRLinY), datastore1->end(colRLinY), coeffA, coeffB);
JKQTPXFunctionLineGraph *graphRegLine=new JKQTPXFunctionLineGraph(plot3);
graphRegLine->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
graphRegLine->setParamsV(coeffA, coeffB);
graphRegLine->setTitle(QString("robust regression: $f(x) = %1 + %2 \\cdot x$, $p=1.1$").arg(jkqtp_floattolatexqstr(coeffA)).arg(jkqtp_floattolatexqstr(coeffB)));
plot3->addGraph(graphRegLine);
// this code can also be written with one function call, using the "adaptor" jkqtpstatAddLinearRegression():
//jkqtpstatAddRobustIRLSLinearRegression(plot3->getPlotter(), datastore1->begin(colRLinX), datastore1->end(colRLinX), datastore1->begin(colRLinY), datastore1->end(colRLinY));
// or even shorter:
//jkqtpstatAddRobustIRLSLinearRegression(graphD);
// as a comparison, we also add the result of the normal/non-robust linear regression:
jkqtpstatAddLinearRegression(graphD);
// the following code demonstrates the influence of the rgularization parameter p:
// - the closer it is to 1, the more robust the fit is (it is closer to the L1-norm)
// - the closer it is to 2, the closer the fit is to the least squares solution (i.e. the normal regression)
double p;
p=1.1;
auto g=jkqtpstatAddRobustIRLSLinearRegression(graphD, nullptr, nullptr, false, false, p);
g->setTitle(g->getTitle()+", $p="+jkqtp_floattolatexqstr(p)+"$");
p=1.5;
g=jkqtpstatAddRobustIRLSLinearRegression(graphD, nullptr, nullptr, false, false, p);
g->setTitle(g->getTitle()+", $p="+jkqtp_floattolatexqstr(p)+"$");
p=1.7;
g=jkqtpstatAddRobustIRLSLinearRegression(graphD, nullptr, nullptr, false, false, p);
g->setTitle(g->getTitle()+", $p="+jkqtp_floattolatexqstr(p)+"$");
p=2;
g=jkqtpstatAddRobustIRLSLinearRegression(graphD, nullptr, nullptr, false, false, p);
g->setTitle(g->getTitle()+", $p="+jkqtp_floattolatexqstr(p)+"$");
// 5.1. The functions for linear regression can also be used to calculate some non-linear models by transforming the input data.
// This is also supported by the statistics library. the supported models are defined in JKQTPStatRegressionModelType
//std::random_device rd; // random number generators:
//std::mt19937 gen{rd()};
//std::normal_distribution<> d1{0,1};
double a0_powerlaw=20;
double b0_powerlaw=0.25;
double a0_exp=5;
double b0_exp=0.5;
size_t colNLLinX=datastore1->addColumn("non-lin data, x");
size_t colNLLinYExp=datastore1->addColumn("non-lin data, y, exponential model");
size_t colNLLinYPow=datastore1->addColumn("non-lin data, y, power-law model");
auto model_powerlaw=jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::PowerLaw);
auto model_exp=jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::Exponential);
for (double x=0.1; x<=10; x+=0.5) {
datastore1->appendToColumn(colNLLinX, x);
double ypow=model_powerlaw(x, a0_powerlaw, b0_powerlaw)+d1(gen);
while (ypow<0) {
ypow=model_powerlaw(x, a0_powerlaw, b0_powerlaw)+d1(gen);
}
datastore1->appendToColumn(colNLLinYPow, ypow);
double yexp=model_exp(x, a0_powerlaw, b0_powerlaw)+d1(gen);
while (yexp<0) {
yexp=model_exp(x, a0_powerlaw, b0_powerlaw)+d1(gen);
}
datastore1->appendToColumn(colNLLinYExp, model_exp(x, a0_exp, b0_exp)+d1(gen));
}
// we visualize this data with a simple scatter graphs:
JKQTPXYLineGraph* graphD_powerlaw;
plot5->addGraph(graphD_powerlaw=new JKQTPXYLineGraph(plot5));
graphD_powerlaw->setXYColumns(colNLLinX, colNLLinYPow);
graphD_powerlaw->setDrawLine(false);
graphD_powerlaw->setTitle(QString("data $%1+\\mathcal{N}(0,1)$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::PowerLaw, a0_powerlaw, b0_powerlaw)));
JKQTPXYLineGraph* graphD_exp;
plot4->addGraph(graphD_exp=new JKQTPXYLineGraph(plot4));
graphD_exp->setXYColumns(colNLLinX, colNLLinYExp);
graphD_exp->setDrawLine(false);
graphD_exp->setTitle(QString("data $%1+\\mathcal{N}(0,1)$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::Exponential, a0_exp, b0_exp)));
// 5.2. Now we calculate the regression models and add a plot to the graph:
double cA=0, cB=0;
JKQTPXFunctionLineGraph* gFunc;
jkqtpstatRegression(JKQTPStatRegressionModelType::Exponential, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYExp), datastore1->end(colNLLinYExp), cA, cB);
plot4->addGraph(gFunc=new JKQTPXFunctionLineGraph(plot4));
gFunc->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::Exponential, cA, cB));
gFunc->setTitle(QString("regression: $%1$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::Exponential, cA, cB)));
cA=0; cB=0;
jkqtpstatRegression(JKQTPStatRegressionModelType::PowerLaw, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYPow), datastore1->end(colNLLinYPow), cA, cB);
plot5->addGraph(gFunc=new JKQTPXFunctionLineGraph(plot5));
gFunc->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType::PowerLaw, cA, cB));
gFunc->setTitle(QString("regression: $%1$").arg(jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType::PowerLaw, cA, cB)));
// Note: Here we used the normal linear regression functions, but variants for IRLS and weighted regression are also available!
// 5.3. Of course also adaptors exist:
//jkqtpstatAddRegression(plot4->getPlotter(), JKQTPStatRegressionModelType::Exponential, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYExp), datastore1->end(colNLLinYExp));
//jkqtpstatAddRegression(plot5->getPlotter(), JKQTPStatRegressionModelType::PowerLaw, datastore1->begin(colNLLinX), datastore1->end(colNLLinX), datastore1->begin(colNLLinYPow), datastore1->end(colNLLinYPow));
//jkqtpstatAddRegression(graphD_exp, JKQTPStatRegressionModelType::Exponential);
//jkqtpstatAddRegression(graphD_powerlaw, JKQTPStatRegressionModelType::PowerLaw);
// 6.1. To demonstrate polynomial fitting, we generate data for a polynomial model
std::vector<double> pPoly {1,2,-2,0.5};
size_t colPolyX=datastore1->addColumn("polynomial data, x");
size_t colPolyY=datastore1->addColumn("polynomial data, y");
for (double x=-10; x<=10; x++) {
datastore1->appendToColumn(colPolyX, x);
datastore1->appendToColumn(colPolyY, jkqtpstatPolyEval(x, pPoly.begin(), pPoly.end())+d1(gen)*50.0);
}
// we visualize this data with a simple scatter graph:
JKQTPXYLineGraph* graphP;
plot6->addGraph(graphP=new JKQTPXYLineGraph(plot6));
graphP->setXYColumns(colPolyX, colPolyY);
graphP->setDrawLine(false);
graphP->setTitle(QString("data $%1+\\mathcal{N}(0,50)$").arg(jkqtpstatPolynomialModel2Latex(pPoly.begin(), pPoly.end())));
// 6.2. now we can fit polynomials with different number of coefficients:
for (size_t p=0; p<=5; p++) {
std::vector<double> pFit;
JKQTPXFunctionLineGraph* gPoly;
jkqtpstatPolyFit(datastore1->begin(colPolyX), datastore1->end(colPolyX), datastore1->begin(colPolyY), datastore1->end(colPolyY), p, std::back_inserter(pFit));
plot6->addGraph(gPoly=new JKQTPXFunctionLineGraph(plot6));
gPoly->setPlotFunctionFunctor(jkqtpstatGeneratePolynomialModel(pFit.begin(), pFit.end()));
gPoly->setTitle(QString("regression: $%1$").arg(jkqtpstatPolynomialModel2Latex(pFit.begin(), pFit.end())));
}
// 6.3. of course also the "adaptor" shortcuts are available:
//for (size_t p=0; p<=5; p++) {
// jkqtpstatAddPolyFit(plot6->getPlotter(), datastore1->begin(colPolyX), datastore1->end(colPolyX), datastore1->begin(colPolyY), datastore1->end(colPolyY), p);
// jkqtpstatAddPolyFit(graphP, p);
//}
// autoscale the plot so the graph is contained
plot1->zoomToFit();
plot1->getXAxis()->setShowZeroAxis(false);
plot1->getYAxis()->setShowZeroAxis(false);
plot1->getPlotter()->setKeyPosition(JKQTPKeyPosition::JKQTPKeyInsideTopLeft);
plot2->zoomToFit();
plot2->getXAxis()->setShowZeroAxis(false);
plot2->getYAxis()->setShowZeroAxis(false);
plot2->getPlotter()->setKeyPosition(JKQTPKeyPosition::JKQTPKeyInsideTopLeft);
plot3->zoomToFit();
plot3->getXAxis()->setShowZeroAxis(false);
plot3->getYAxis()->setShowZeroAxis(false);
plot3->getPlotter()->setKeyPosition(JKQTPKeyPosition::JKQTPKeyInsideTopLeft);
plot4->zoomToFit();
plot4->getXAxis()->setShowZeroAxis(false);
plot4->getYAxis()->setShowZeroAxis(false);
plot4->getPlotter()->setKeyPosition(JKQTPKeyPosition::JKQTPKeyInsideTopLeft);
plot4->setAbsoluteX(0.05, plot4->getXMax());
plot4->zoomToFit();
plot5->getXAxis()->setShowZeroAxis(false);
plot5->getYAxis()->setShowZeroAxis(false);
plot5->getPlotter()->setKeyPosition(JKQTPKeyPosition::JKQTPKeyInsideTopLeft);
plot5->setAbsoluteX(0.05, plot5->getXMax());
plot5->zoomToFit();
plot6->getXAxis()->setShowZeroAxis(false);
plot6->getYAxis()->setShowZeroAxis(false);
plot6->getPlotter()->setKeyPosition(JKQTPKeyPosition::JKQTPKeyInsideBottomRight);
plot6->zoomToFit();
// show plotter and make it a decent size
mainWidget.show();
mainWidget.resize(1600,800);
return app.exec();
}

View File

@ -0,0 +1,27 @@
# source code for this simple demo
SOURCES = jkqtplotter_simpletest_datastore_regression.cpp
# configure Qt
CONFIG += link_prl qt
QT += core gui xml svg
greaterThan(QT_MAJOR_VERSION, 4): QT += widgets printsupport
# output executable name
TARGET = jkqtplotter_simpletest_datastore_regression
# include JKQTPlotter source headers and link against library
DEPENDPATH += ../../lib ../../staticlib/jkqtplotterlib
INCLUDEPATH += ../../lib
CONFIG (debug, debug|release) {
LIBS += -L../../staticlib/jkqtplotterlib/debug -ljkqtplotterlib_debug
} else {
LIBS += -L../../staticlib/jkqtplotterlib/release -ljkqtplotterlib
}
message("LIBS = $$LIBS")
win32-msvc*: DEFINES += _USE_MATH_DEFINES
win32-msvc*: DEFINES += NOMINMAX

View File

@ -0,0 +1,8 @@
TEMPLATE = subdirs
SUBDIRS += jkqtplotterlib jkqtplotter_simpletest_datastore_regression
jkqtplotterlib.file = ../../staticlib/jkqtplotterlib/jkqtplotterlib.pro
jkqtplotter_simpletest_datastore_regression.file=$$PWD/jkqtplotter_simpletest_datastore_regression.pro
jkqtplotter_simpletest_datastore_regression.depends = jkqtplotterlib

View File

@ -1,8 +1,10 @@
# Tutorial (JKQTPDatastore): Advanced 1-Dimensional Statistics with JKQTPDatastore {#JKQTPlotterBasicJKQTPDatastoreStatistics}
[JKQTPlotterBasicJKQTPDatastore]: @ref JKQTPlotterBasicJKQTPDatastore "Basic Usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreIterators]: @ref JKQTPlotterBasicJKQTPDatastoreIterators "Iterator-Based usage of JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreStatistics]: @ref JKQTPlotterBasicJKQTPDatastoreStatistics "Advanced 1-Dimensional Statistics with JKQTPDatastore"
[JKQTPlotterBasicJKQTPDatastoreRegression]: @ref JKQTPlotterBasicJKQTPDatastoreRegression "Regression Analysis (with the Statistics Library)"
[statisticslibrary]: @ref jkqtptools_math_statistics "JKQTPlotter Statistics Library"
@ -12,6 +14,7 @@ This tutorial project (see `./examples/simpletest_datastore_statistics/`) explai
- [JKQTPlotterBasicJKQTPDatastore]
- [JKQTPlotterBasicJKQTPDatastoreIterators]
- [JKQTPlotterBasicJKQTPDatastoreStatistics]
- [JKQTPlotterBasicJKQTPDatastoreRegression]
[TOC]

View File

@ -389,13 +389,13 @@ inline void jkqtplinalgMatrixProduct(const T* M1, long L1, long C1, const T* M2,
}
} else if (M1==M && M2!=M) {
JKQTPArrayScopedPointer<T> MM(jkqtpArrayDuplicate(M1, L1*C1));
jkqtplinalgMatrixProduct(MM,L1,C1,M2,L2,C2,M);
jkqtplinalgMatrixProduct(MM.data(),L1,C1,M2,L2,C2,M);
} else if (M1!=M && M2==M) {
JKQTPArrayScopedPointer<T> MM(jkqtpArrayDuplicate(M1, L1*C1));
jkqtplinalgMatrixProduct(M1,L1,C1,MM,L2,C2,M);
jkqtplinalgMatrixProduct(M1,L1,C1,MM.data(),L2,C2,M);
} else if (M1==M && M2==M) {
JKQTPArrayScopedPointer<T> MM(jkqtpArrayDuplicate(M1, L1*C1));
jkqtplinalgMatrixProduct(MM,L1,C1,MM,L2,C2,M);
jkqtplinalgMatrixProduct(MM.data(),L1,C1,MM.data(),L2,C2,M);
}
}
@ -581,7 +581,7 @@ inline bool jkqtplinalgMatrixInversion(const T* matrix, T* matrix_out, long N) {
}
bool ok=linalgGaussJordanV(m, N, 2*N);
bool ok=jkqtplinalgGaussJordan(m.data(), N, 2*N);
if (ok) {
// finally we copy the result to matrix_out
@ -714,7 +714,7 @@ inline bool jkqtplinalgLinSolve(const T* A, const T* B, long N, long C, T* resul
}
bool ok=linalgGaussJordan(m.data(), N, N+C);
bool ok=jkqtplinalgGaussJordan(m.data(), N, N+C);
if (ok) {
for (long k=0; k<N; k++) {

View File

@ -213,16 +213,39 @@ inline T jkqtp_cube(T x) {
}
/*! \brief calculates the sign of number \a x
/*! \brief calculates the sign of number \a x (-1 for x<0 and +1 for x>=0)
\ingroup jkqtptools_math_basic
*/
template <class T>
inline T jkqtp_sign(T x) {
if (x<0) return -1;
//else if (x==0) return 0;
else return 1;
}
/** \brief returns the inversely proportional value 1/\a v of \a v
* \ingroup jkqtptools_math_basic */
template<typename T>
inline T jkqtp_inverseProp(const T& v) {
return T(1.0)/v;
}
/** \brief returns the inversely proportional value 1/\a v of \a v and ensures that \f$ |v|\geq \mbox{absMinV} \f$
* \ingroup jkqtptools_math_basic */
template<typename T>
inline T jkqtp_inversePropSave(const T& v, const T& absMinV) {
T vv=v;
if (fabs(vv)<absMinV) vv=jkqtp_sign(v)*absMinV;
return T(1.0)/vv;
}
/** \brief returns the inversely proportional value 1/\a v of \a v and ensures that \f$ |v|\geq \mbox{absMinV} \f$, uses \c absMinV=std::numeric_limits<T>::epsilon()*100.0
* \ingroup jkqtptools_math_basic */
template<typename T>
inline T jkqtp_inversePropSaveDefault(const T& v) {
return jkqtp_inversePropSave<T>(v, std::numeric_limits<T>::epsilon()*100.0);
}
/** \brief calculate the distance between two QPointF points
* \ingroup jkqtptools_math_basic
*

View File

@ -104,3 +104,61 @@ double JKQTPStat5NumberStatistics::IQR() const {
double JKQTPStat5NumberStatistics::IQRSignificanceEstimate() const {
return 2.0*(1.58*(IQR()))/sqrt(static_cast<double>(N));
}
std::function<double (double, double, double)> jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType type) {
switch(type) {
case JKQTPStatRegressionModelType::Linear: return [](double x, double a, double b)->double { return a+b*x; };
case JKQTPStatRegressionModelType::PowerLaw: return [](double x, double a, double b)->double { return a*pow(x,b); };
case JKQTPStatRegressionModelType::Exponential: return [](double x, double a, double b)->double { return a*exp(b*x); };
}
throw std::runtime_error("unknown JKQTPStatRegressionModelType in jkqtpStatGenerateRegressionModel()");
}
QString jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType type, double a, double b) {
switch(type) {
case JKQTPStatRegressionModelType::Linear: return QString("f(x)=%1+%2{\\cdot}x").arg(jkqtp_floattolatexqstr(a, 3)).arg(jkqtp_floattolatexqstr(b, 3));
case JKQTPStatRegressionModelType::PowerLaw: return QString("f(x)=%1{\\cdot}x^{%2}").arg(jkqtp_floattolatexqstr(a, 3)).arg(jkqtp_floattolatexqstr(b, 3));
case JKQTPStatRegressionModelType::Exponential: return QString("f(x)=%1{\\cdot}\\exp(%2{\\cdot}x)").arg(jkqtp_floattolatexqstr(a, 3)).arg(jkqtp_floattolatexqstr(b, 3));
}
throw std::runtime_error("unknown JKQTPStatRegressionModelType in jkqtpstatRegressionModel2Latex()");
}
std::function<double (double)> jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType type, double a, double b) {
auto res=jkqtpStatGenerateRegressionModel(type);
return std::bind(res, std::placeholders::_1, a, b);
}
std::pair<std::function<double (double)>, std::function<double (double)> > jkqtpStatGenerateTransformation(JKQTPStatRegressionModelType type) {
auto logF=[](double x)->double { return log(x); };
auto idF=&jkqtp_identity<double>;
switch(type) {
case JKQTPStatRegressionModelType::Linear: return std::pair<std::function<double(double)>,std::function<double(double)> >(idF, idF);
case JKQTPStatRegressionModelType::PowerLaw: return std::pair<std::function<double(double)>,std::function<double(double)> >(logF, logF);
case JKQTPStatRegressionModelType::Exponential: return std::pair<std::function<double(double)>,std::function<double(double)> >(idF, logF);
}
throw std::runtime_error("unknown JKQTPStatRegressionModelType in jkqtpStatGenerateTransformation()");
}
std::pair<std::function<double (double)>, std::function<double (double)> > jkqtpStatGenerateParameterATransformation(JKQTPStatRegressionModelType type) {
auto logF=[](double x)->double { return log(x); };
auto expF=[](double x)->double { return exp(x); };
auto idF=&jkqtp_identity<double>;
switch(type) {
case JKQTPStatRegressionModelType::Linear: return std::pair<std::function<double(double)>,std::function<double(double)> >(idF, idF);
case JKQTPStatRegressionModelType::PowerLaw: return std::pair<std::function<double(double)>,std::function<double(double)> >(logF, expF);
case JKQTPStatRegressionModelType::Exponential: return std::pair<std::function<double(double)>,std::function<double(double)> >(logF, expF);
}
throw std::runtime_error("unknown JKQTPStatRegressionModelType in jkqtpStatGenerateTransformation()");
}
std::pair<std::function<double (double)>, std::function<double (double)> > jkqtpStatGenerateParameterBTransformation(JKQTPStatRegressionModelType type) {
//auto logF=[](double x)->double { return log(x); };
//auto expF=[](double x)->double { return exp(x); };
auto idF=&jkqtp_identity<double>;
switch(type) {
case JKQTPStatRegressionModelType::Linear: return std::pair<std::function<double(double)>,std::function<double(double)> >(idF, idF);
case JKQTPStatRegressionModelType::PowerLaw: return std::pair<std::function<double(double)>,std::function<double(double)> >(idF, idF);
case JKQTPStatRegressionModelType::Exponential: return std::pair<std::function<double(double)>,std::function<double(double)> >(idF, idF);
}
throw std::runtime_error("unknown JKQTPStatRegressionModelType in jkqtpStatGenerateTransformation()");
}

View File

@ -37,13 +37,7 @@
#include "jkqtcommon/jkqtp_imexport.h"
#include "jkqtcommon/jkqtplinalgtools.h"
#include "jkqtcommon/jkqtparraytools.h"
#ifdef _OPENMP
# include <omp.h>
#endif
#include "jkqtcommon/jkqtptoolsdebugging.h"
@ -1801,6 +1795,612 @@ inline void jkqtpstatKDE1D(InputIt first, InputIt last, double binXLeft, double
/*! \brief calculate the linear regression coefficients for a given data range \a firstX / \a firstY ... \a lastX / \a lastY where the model is \f$ f(x)=a+b\cdot x \f$
So this function solves the least-squares optimization problem: \f[ (a^\ast, b^\ast)=\argmin\limits_{a,b}\sum\limits_i\left(y_i-(a+b\cdot x_i)\right)^2 \f]
\ingroup jkqtptools_math_statistics_regression
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used
This function computes internally:
\f[ a=\overline{y}-b\cdot\overline{x} \f]
\f[ b=\frac{\sum x_iy_i-N\cdot\overline{x}\cdot\overline{y}}{\sum x_i^2-N\cdot(\overline{x})^2} \f]
\image html jkqtplotter_simpletest_datastore_regression_lin.png
*/
template <class InputItX, class InputItY>
inline void jkqtpstatLinearRegression(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, double& coeffA, double& coeffB, bool fixA=false, bool fixB=false) {
if (fixA&&fixB) return;
const int Nx=std::distance(firstX,lastX);
const int Ny=std::distance(firstY,lastY);
JKQTPASSERT(Nx>1 && Ny>1);
double sumx=0, sumy=0, sumxy=0, sumx2=0;
size_t N=0;
auto itX=firstX;
auto itY=firstY;
for (; itX!=lastX && itY!=lastY; ++itX, ++itY) {
const double fit_x=jkqtp_todouble(*itX);
const double fit_y=jkqtp_todouble(*itY);
if (JKQTPIsOKFloat(fit_x) && JKQTPIsOKFloat(fit_y)) {
sumx=sumx+fit_x;
sumy=sumy+fit_y;
sumxy=sumxy+fit_x*fit_y;
sumx2=sumx2+fit_x*fit_x;
N++;
}
}
const double NN=static_cast<double>(N);
JKQTPASSERT_M(NN>1, "too few datapoints");
if (!fixA && !fixB) {
coeffB=(double(sumxy)-double(sumx)*double(sumy)/NN)/(double(sumx2)-double(sumx)*double(sumx)/NN);;
coeffA=double(sumy)/NN-coeffB*double(sumx)/NN;
} else if (fixA && !fixB) {
coeffB=(double(sumy)/NN-coeffA)/(double(sumx)/NN);
} else if (!fixA && fixB) {
coeffA=double(sumy)/NN-coeffB*double(sumx)/NN;
}
}
/*! \brief calculate the weighted linear regression coefficients for a given for a given data range \a firstX / \a firstY / \a firstW ... \a lastX / \a lastY / \a lastW where the model is \f$ f(x)=a+b\cdot x \f$
So this function solves the least-squares optimization problem: \f[ (a^\ast, b^\ast)=\argmin\limits_{a,b}\sum\limits_iw_i^2\cdot\left(y_i-(a+b\cdot x_i)\right)^2 \f]
\ingroup jkqtptools_math_statistics_regression
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\tparam InputItW standard iterator type of \a firstW and \a lastW.
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param firstW iterator pointing to the first item in the weight-dataset to use \f$ w_1 \f$
\param lastW iterator pointing behind the last item in the weight-dataset to use \f$ w_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used
\param fWeightDataToWi an optional function, which is applied to the data from \a firstW ... \a lastW to convert them to weight, i.e. \c wi=fWeightDataToWi(*itW)
e.g. if you use data used to draw error bars, you can use jkqtp_inversePropSaveDefault(). The default is jkqtp_identity(), which just returns the values.
In the case of jkqtp_inversePropSaveDefault(), a datapoint x,y, has a large weight, if it's error is small and in the case if jkqtp_identity() it's weight
is directly proportional to the given value.
This function internally computes:
\f[ a=\frac{\overline{y}-b\cdot\overline{x}}{\overline{w^2}} \f]
\f[ b=\frac{\overline{w^2}\cdot\overline{x\cdot y}-\overline{x}\cdot\overline{y}}{\overline{x^2}\cdot\overline{w^2}-\overline{x}^2} \f]
Here the averages are defined in terms of a weight vector \f$ w_i\f$:
\f[ \overline{x}=\sum\limits_iw_i^2\cdot x_i \f]
\f[ \overline{y}=\sum\limits_iw_i^2\cdot y_i \f]
\f[ \overline{x\cdot y}=\sum\limits_iw_i^2\cdot x_i\cdot y_i \f]
\f[ \overline{x^2}=\sum\limits_iw_i^2\cdot x_i^2 \f]
\f[ \overline{w^2}=\sum\limits_iw_i^2 \f]
\image html jkqtplotter_simpletest_datastore_regression_linweight.png
*/
template <class InputItX, class InputItY, class InputItW>
inline 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>) {
if (fixA&&fixB) return;
const int Nx=std::distance(firstX,lastX);
const int Ny=std::distance(firstY,lastY);
const int Nw=std::distance(firstW,lastW);
JKQTPASSERT(Nx>1 && Ny>1 && Nw>1);
double sumx=0, sumy=0, sumxy=0, sumx2=0, sumw2=0;
size_t N=0;
auto itX=firstX;
auto itY=firstY;
auto itW=firstW;
for (; itX!=lastX && itY!=lastY && itW!=lastW; ++itX, ++itY, ++itW) {
const double fit_x=jkqtp_todouble(*itX);
const double fit_y=jkqtp_todouble(*itY);
const double fit_w2=jkqtp_sqr(fWeightDataToWi(jkqtp_todouble(*itW)));
if (JKQTPIsOKFloat(fit_x)&&JKQTPIsOKFloat(fit_y)&&JKQTPIsOKFloat(fit_w2)) {
sumx=sumx+fit_w2*fit_x;
sumy=sumy+fit_w2*fit_y;
sumxy=sumxy+fit_w2*fit_x*fit_y;
sumx2=sumx2+fit_w2*fit_x*fit_x;
sumw2=sumw2+fit_w2;
N++;
}
}
const double NN=static_cast<double>(N);
JKQTPASSERT_M(NN>1, "too few datapoints");
if (!fixA && !fixB) {
coeffB=(double(sumxy)*double(sumw2)-double(sumx)*double(sumy))/(double(sumx2)*double(sumw2)-double(sumx)*double(sumx));
coeffA=(double(sumy)-coeffB*double(sumx))/double(sumw2);
} else if (fixA && !fixB) {
coeffB=(double(sumy)-coeffA*double(sumw2))/double(sumx);
} else if (!fixA && fixB) {
coeffA=(double(sumy)-coeffB*double(sumx))/double(sumw2);
}
}
/*! \brief calculate the (robust) iteratively reweighted least-squares (IRLS) estimate for the parameters of the model \f$ f(x)=a+b\cdot x \f$
for a given data range \a firstX / \a firstY ... \a lastX / \a lastY
So this function finds an outlier-robust solution to the optimization problem:
\f[ (a^\ast,b^\ast)=\argmin\limits_{a,b}\sum\limits_i|a+b\cdot x_i-y_i|^p \f]
\ingroup jkqtptools_math_statistics_regression
\ingroup jkqtptools_math_statistics_regression
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used
\param p regularization parameter, the optimization problem is formulated in the \f$ L_p \f$ norm, using this \a p (see image below for an example)
\param iterations the number of iterations the IRLS algorithm performs
This is a simple form of the IRLS algorithm to estimate the parameters a and b in a linear model \f$ f(x)=a+b\cdot x \f$.
This algorithm solves the optimization problem for a \f$ L_p\f$-norm:
\f[ (a^\ast,b^\ast)=\argmin\limits_{a,b}\sum\limits_i|a+b\cdot x_i-y_i|^p \f]
by iteratively optimization weights \f$ \vec{w} \f$ and solving a weighted least squares problem in each iteration:
\f[ (a_n,b_n)=\argmin\limits_{a,b}\sum\limits_i|a+b\cdot x_i-y_i|^{(p-2)}\cdot|a+b\cdot x_i-y_i|^2 \f]
The IRLS-algorithm works as follows:
- calculate initial \f$ a_0\f$ and \f$ b_0\f$ with unweighted regression from x and y
- perform a number of iterations (parameter \a iterations ). In each iteration \f$ n\f$:
- calculate the error vector \f$\vec{e}\f$: \f[ e_i = a+b\cdot x_i -y_i \f]
- estimate new weights \f$\vec{w}\f$: \[ w_i=|e_i|^{(p-2)/2} \]
- calculate new estimates \f$ a_n\f$ and \f$ b_n\f$ with weighted regression from \f$ \vec{x}\f$ and \f$ \vec{y}\f$ and \f$ \vec{w}\f$
.
- return the last estimates \f$ a_n\f$ and \f$ b_n\f$
.
\image html irls.png
\image html jkqtplotter_simpletest_datastore_regression_linrobust_p.png
\see https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares, C. Sidney Burrus: "Iterative Reweighted Least Squares", <a href="http://cnx.org/content/m45285/latest/">http://cnx.org/content/m45285/latest/</a>
*/
template <class InputItX, class InputItY>
inline 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) {
if (fixA&&fixB) return;
const int Nx=std::distance(firstX,lastX);
const int Ny=std::distance(firstY,lastY);
const int N=std::min(Nx,Ny);
JKQTPASSERT(Nx>1 && Ny>1);
std::vector<double> weights;
std::fill_n(std::back_inserter(weights), N, 1.0);
double alast=coeffA, blast=coeffB;
jkqtpstatLinearWeightedRegression(firstX, lastX, firstY, lastY, weights.begin(), weights.end(), alast, blast, fixA, fixB, &jkqtp_identity<double>);
for (int it=0; it<iterations-1; it++) {
// calculate weights
auto itX=firstX;
auto itY=firstY;
for (double& w: weights) {
const double fit_x=*itX;
const double fit_y=*itY;
const double e=alast+blast*fit_x-fit_y;
w=pow(std::max<double>(JKQTP_EPSILON*100.0, fabs(e)), (p-2.0)/2.0);
++itX;
++itY;
}
// solve weighted linear least squares
jkqtpstatLinearWeightedRegression(firstX, lastX, firstY, lastY, weights.begin(), weights.end(), alast, blast, fixA, fixB, &jkqtp_identity<double>);
}
coeffA=alast;
coeffB=blast;
}
/*! \brief when performing linear regression, different target functions can be fitted, if the input data is transformed accordingly. This library provides the options in this enum by default.
\ingroup jkqtptools_math_statistics_regression
*/
enum class JKQTPStatRegressionModelType {
Linear, /*!< \brief linear model \f$ f(x)=a+b\cdot x \f$ */
PowerLaw, /*!< \brief power law model \f$ f(x)=a\cdot x^b \f$ */
Exponential, /*!< \brief exponential model \f$ f(x)=a\cdot \exp(b\cdot x) \f$ */
};
/*! \brief Generates functors \c f(x,a,b) for the models from JKQTPStatRegressionModelType in \a type
\ingroup jkqtptools_math_statistics_regression
*/
JKQTP_LIB_EXPORT std::function<double(double, double, double)> jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType type);
/*! \brief Generates a LaTeX string for the models from JKQTPStatRegressionModelType in \a type
\ingroup jkqtptools_math_statistics_regression
*/
JKQTP_LIB_EXPORT QString jkqtpstatRegressionModel2Latex(JKQTPStatRegressionModelType type, double a, double b);
/*! \brief Generates functors \c f(x) for the models from JKQTPStatRegressionModelType in \a type and binds the parameter values \a and \a b to the returned function
\ingroup jkqtptools_math_statistics_regression
*/
JKQTP_LIB_EXPORT std::function<double(double)> jkqtpStatGenerateRegressionModel(JKQTPStatRegressionModelType type, double a, double b);
/*! \brief Generates the transformation function for x-data (\c result.first ) and y-data (\c result.second ) for each regression model in JKQTPStatRegressionModelType in \a type
\ingroup jkqtptools_math_statistics_regression
\internal
*/
JKQTP_LIB_EXPORT std::pair<std::function<double(double)>,std::function<double(double)> > jkqtpStatGenerateTransformation(JKQTPStatRegressionModelType type);
/*! \brief Generates the transformation function for a-parameter (offset, \c result.first : transform, \c result.second : back-transform) for each regression model in JKQTPStatRegressionModelType in \a type
\ingroup jkqtptools_math_statistics_regression
\internal
*/
JKQTP_LIB_EXPORT std::pair<std::function<double(double)>,std::function<double(double)> > jkqtpStatGenerateParameterATransformation(JKQTPStatRegressionModelType type);
/*! \brief Generates the transformation function for b-parameter (slope, \c result.first : transform, \c result.second : back-transform) for each regression model in JKQTPStatRegressionModelType in \a type
\ingroup jkqtptools_math_statistics_regression
\internal
*/
JKQTP_LIB_EXPORT std::pair<std::function<double(double)>,std::function<double(double)> > jkqtpStatGenerateParameterBTransformation(JKQTPStatRegressionModelType type);
/*! \brief calculate the linear regression coefficients for a given data range \a firstX / \a firstY ... \a lastX / \a lastY where the model is defined by \a type
So this function solves the least-squares optimization problem: \f[ (a^\ast, b^\ast)=\argmin\limits_{a,b}\sum\limits_i\left(y_i-f_\mbox{type}(x_i,a,b)\right)^2 \f]
by reducing it to a linear fit by transforming x- and/or y-data
\ingroup jkqtptools_math_statistics_regression
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used
This function computes internally first transforms the data, as appropriate to fit the model defined by \a type and then calls jkqtpstatLinearRegression()
to obtain the parameters. The output parameters are transformed, so they can be used with jkqtpStatGenerateRegressionModel() to generate a functor
that evaluates the model
\see JKQTPStatRegressionModelType, jkqtpStatGenerateRegressionModel(), jkqtpstatLinearRegression(), jkqtpStatGenerateTransformation()
*/
template <class InputItX, class InputItY>
inline void jkqtpstatRegression(JKQTPStatRegressionModelType type, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, double& coeffA, double& coeffB, bool fixA=false, bool fixB=false) {
std::vector<double> x, y;
auto trafo=jkqtpStatGenerateTransformation(type);
auto aTrafo =jkqtpStatGenerateParameterATransformation(type);
auto bTrafo =jkqtpStatGenerateParameterBTransformation(type);
std::transform(firstX, lastX, std::back_inserter(x), trafo.first);
std::transform(firstY, lastY, std::back_inserter(y), trafo.second);
double a=aTrafo.first(coeffA);
double b=bTrafo.first(coeffB);
jkqtpstatLinearRegression(x.begin(), x.end(), y.begin(), y.end(), a, b, fixA, fixB);
coeffA=aTrafo.second(a);
coeffB=bTrafo.second(b);
}
/*! \brief calculate the robust linear regression coefficients for a given data range \a firstX / \a firstY ... \a lastX / \a lastY where the model is defined by \a type
So this function solves the Lp-norm optimization problem: \f[ (a^\ast, b^\ast)=\argmin\limits_{a,b}\sum\limits_i\left(y_i-f_\mbox{type}(x_i,a,b)\right)^p \f]
by reducing it to a linear fit by transforming x- and/or y-data
\ingroup jkqtptools_math_statistics_regression
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used
\param p regularization parameter, the optimization problem is formulated in the \f$ L_p \f$ norm, using this \a p (see image below for an example)
\param iterations the number of iterations the IRLS algorithm performs
This function computes internally first transforms the data, as appropriate to fit the model defined by \a type and then calls jkqtpstatRobustIRLSLinearRegression()
to obtain the parameters. The output parameters are transformed, so they can be used with jkqtpStatGenerateRegressionModel() to generate a functor
that evaluates the model
\see JKQTPStatRegressionModelType, jkqtpStatGenerateRegressionModel(), jkqtpstatRobustIRLSLinearRegression(), jkqtpStatGenerateTransformation()
*/
template <class InputItX, class InputItY>
inline 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) {
std::vector<double> x, y;
auto trafo=jkqtpStatGenerateTransformation(type);
auto aTrafo =jkqtpStatGenerateParameterATransformation(type);
auto bTrafo =jkqtpStatGenerateParameterBTransformation(type);
std::transform(firstX, lastX, std::back_inserter(x), trafo.first);
std::transform(firstY, lastY, std::back_inserter(y), trafo.second);
double a=aTrafo.first(coeffA);
double b=bTrafo.first(coeffB);
jkqtpstatRobustIRLSLinearRegression(x.begin(), x.end(), y.begin(), y.end(), a, b, fixA, fixB, p, iterations);
coeffA=aTrafo.second(a);
coeffB=bTrafo.second(b);
}
/*! \brief calculate the robust linear regression coefficients for a given data range \a firstX / \a firstY ... \a lastX / \a lastY where the model is defined by \a type
So this function solves the Lp-norm optimization problem: \f[ (a^\ast, b^\ast)=\argmin\limits_{a,b}\sum\limits_i\left(y_i-f_\mbox{type}(x_i,a,b)\right)^p \f]
by reducing it to a linear fit by transforming x- and/or y-data
\ingroup jkqtptools_math_statistics_regression
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\tparam InputItW standard iterator type of \a firstW and \a lastW.
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param firstW iterator pointing to the first item in the weight-dataset to use \f$ w_1 \f$
\param lastW iterator pointing behind the last item in the weight-dataset to use \f$ w_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used
\param fWeightDataToWi an optional function, which is applied to the data from \a firstW ... \a lastW to convert them to weight, i.e. \c wi=fWeightDataToWi(*itW)
e.g. if you use data used to draw error bars, you can use jkqtp_inversePropSaveDefault(). The default is jkqtp_identity(), which just returns the values.
In the case of jkqtp_inversePropSaveDefault(), a datapoint x,y, has a large weight, if it's error is small and in the case if jkqtp_identity() it's weight
is directly proportional to the given value.
This function computes internally first transforms the data, as appropriate to fit the model defined by \a type and then calls jkqtpstatLinearWeightedRegression()
to obtain the parameters. The output parameters are transformed, so they can be used with jkqtpStatGenerateRegressionModel() to generate a functor
that evaluates the model
\see JKQTPStatRegressionModelType, jkqtpStatGenerateRegressionModel(), jkqtpstatLinearWeightedRegression(), jkqtpStatGenerateTransformation()
*/
template <class InputItX, class InputItY, class InputItW>
inline 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>) {
std::vector<double> x, y;
auto trafo=jkqtpStatGenerateTransformation(type);
auto aTrafo =jkqtpStatGenerateParameterATransformation(type);
auto bTrafo =jkqtpStatGenerateParameterBTransformation(type);
std::transform(firstX, lastX, std::back_inserter(x), trafo.first);
std::transform(firstY, lastY, std::back_inserter(y), trafo.second);
double a=aTrafo.first(coeffA);
double b=bTrafo.first(coeffB);
jkqtpstatLinearWeightedRegression(x.begin(), x.end(), y.begin(), y.end(), firstW, lastW, a, b, fixA, fixB, fWeightDataToWi);
coeffA=aTrafo.second(a);
coeffB=bTrafo.second(b);
}
/*! \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$
\ingroup jkqtptools_math_statistics_poly
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\tparam OutputItP output iterator for the polynomial coefficients
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param P degree of the polynomial (P>=N !!!)
\param[out] firstRes Iterator (of type \a OutputItP ), which receives the (P+1)-entry vector with the polynomial coefficients \f$ p_i \f$
This function uses jkqtpstatLinSolve() to solve the system of equations
\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]
\f[ \vec{y}=V\vec{p}\ \ \ \ \ \Rightarrow\ \ \ \ \ \vec{p}=(V^TV)^{-1}V^T\vec{y} \]
\image html jkqtplotter_simpletest_datastore_regression_polynom.png
\see https://en.wikipedia.org/wiki/Polynomial_regression
*/
template <class InputItX, class InputItY, class OutputItP>
inline void jkqtpstatPolyFit(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, size_t P, OutputItP firstRes) {
{
const int Nx=std::distance(firstX,lastX);
const int Ny=std::distance(firstY,lastY);
JKQTPASSERT(Nx>1 && Ny>1);
}
size_t N=0;
std::vector<double> X,Y;
auto itX=firstX;
auto itY=firstY;
for (; itX!=lastX && itY!=lastY; ++itX, ++itY) {
const double fit_x=jkqtp_todouble(*itX);
const double fit_y=jkqtp_todouble(*itY);
if (JKQTPIsOKFloat(fit_x) && JKQTPIsOKFloat(fit_y)) {
X.push_back(fit_x);
Y.push_back(fit_y);
N++;
}
}
// build Vandermonde matrix V
std::vector<double> V;
V.resize(N*(P+1));
for (size_t l=0; l<N; l++) {
V[jkqtplinalgMatIndex(l,0,P+1)]=1.0;
double x=X[l];
const double xx=x;
for (size_t c=1; c<P+1; c++) {
V[jkqtplinalgMatIndex(l,c,P+1)]=x;
x=x*xx;
}
}
#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
std::cout<<"V = \n";
jkqtplinalgPrintMatrix(V.data(),N,P+1);
std::cout<<"\n";
#endif
// calculate V^T
std::vector<double> VT=V;
jkqtplinalgTransposeMatrix(VT.data(), static_cast<long>(N), static_cast<long>(P+1));
#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
std::cout<<"V^T = \n";
jkqtplinalgPrintMatrix(VT.data(),P+1,N);
std::cout<<"\n";
#endif
// calculate V^T*V
std::vector<double> VTV;
VTV.resize((P+1)*(P+1));
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());
#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
std::cout<<"V^T*V = \n";
jkqtplinalgPrintMatrix(VTV.data(),P+1,P+1);
std::cout<<"\n";
#endif
// calculate V^T*y
std::vector<double> VTY;
VTY.resize(P+1);
jkqtplinalgMatrixProduct(VT.data(), static_cast<long>(P+1), static_cast<long>(N), Y.data(), static_cast<long>(N), 1, VTY.data());
#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
std::cout<<"V^T*y = \n";
jkqtplinalgPrintMatrix(VTY.data(),P+1,1);
std::cout<<"\n";
#endif
// solve V^T*y = V^T*V*p
const bool ok=jkqtplinalgLinSolve(VTV.data(), VTY.data(), static_cast<long>(P+1));
if (ok) {
auto itR=firstRes;
for (size_t p=0; p<P+1; p++) {
*++itR=VTY[p];
}
} else {
throw std::runtime_error("jkqtplinalgLinSolve() didn't return a result!");
}
#ifdef STATISTICS_TOOLS_DEBUG_statisticsPolyFit
std::cout<<"result_out = \n";
jkqtplinalgPrintMatrix(result_out,P+1,1);
std::cout<<"\n";
#endif
}
/*! \brief evaluate a polynomial \f$ f(x)=\sum\limits_{i=0}^Pp_ix^i \f$ with \f$ p_i \f$ taken from the range \a firstP ... \a lastP
\ingroup jkqtptools_math_statistics_poly
\tparam PolyItP iterator for the polynomial coefficients
\param x where to evaluate
\param firstP points to the first polynomial coefficient \f$ p_1 \f$ (i.e. the offset with \f$ x^0 \f$ )
\param lastP points behind the last polynomial coefficient \f$ p_P \f$
\return value of polynomial \f$ f(x)=\sum\limits_{i=0}^Pp_ix^i \f$ at location \a x
*/
template <class PolyItP>
inline double jkqtpstatPolyEval(double x, PolyItP firstP, PolyItP lastP) {
double v=0.0;
double xx=1.0;
for (auto itP=firstP; itP!=lastP; ++itP) {
v=v+(*itP)*xx;
xx=xx*x;
}
return v;
}
/*! \brief a C++-functor, which evaluates a polynomial
\ingroup jkqtptools_math_statistics_poly
*/
struct JKQTPStatPolynomialFunctor {
std::vector<double> P;
template <class PolyItP>
inline JKQTPStatPolynomialFunctor(PolyItP firstP, PolyItP lastP) {
for (auto itP=firstP; itP!=lastP; ++itP) {
P.push_back(*itP);
}
}
inline double operator()(double x) const { return jkqtpstatPolyEval(x, P.begin(), P.end()); }
};
/*! \brief returns a C++-functor, which evaluates a polynomial
\ingroup jkqtptools_math_statistics_poly
\tparam PolyItP iterator for the polynomial coefficients
\param firstP points to the first polynomial coefficient \f$ p_1 \f$ (i.e. the offset with \f$ x^0 \f$ )
\param lastP points behind the last polynomial coefficient \f$ p_P \f$
*/
template <class PolyItP>
inline std::function<double(double)> jkqtpstatGeneratePolynomialModel(PolyItP firstP, PolyItP lastP) {
return JKQTPStatPolynomialFunctor(firstP, lastP);
}
/*! \brief Generates a LaTeX string for the polynomial model with the coefficients \a firstP ... \a lastP
\ingroup jkqtptools_math_statistics_regression
\tparam PolyItP iterator for the polynomial coefficients
\param firstP points to the first polynomial coefficient \f$ p_1 \f$ (i.e. the offset with \f$ x^0 \f$ )
\param lastP points behind the last polynomial coefficient \f$ p_P \f$
*/
template <class PolyItP>
QString jkqtpstatPolynomialModel2Latex(PolyItP firstP, PolyItP lastP) {
QString str="f(x)=";
size_t p=0;
for (auto itP=firstP; itP!=lastP; ++itP) {
if (p==0) str+=jkqtp_floattolatexqstr(*itP, 3);
else {
if (*itP>=0) str+="+";
str+=QString("%2{\\cdot}x^{%1}").arg(p).arg(jkqtp_floattolatexqstr(*itP, 3));
}
p++;
}
return str;
}
#endif // JKQTPSTATISTICSTOOLS_H_INCLUDED

View File

@ -21,4 +21,58 @@
#include "jkqtplotter/jkqtpgraphsstatisticsadaptors.h"
#include "jkqtplotter/jkqtplotter.h"
#include "jkqtplotter/jkqtpgraphsbase.h"
JKQTPXFunctionLineGraph *jkqtpstatAddLinearRegression(JKQTPXYGraph *datagraph, double *coeffA, double *coeffB, bool fixA, bool fixB) {
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
return jkqtpstatAddLinearRegression(plt, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()), coeffA, coeffB, fixA, fixB);
}
JKQTPXFunctionLineGraph *jkqtpstatAddLinearWeightedRegression(JKQTPXYGraph *datagraph, double *coeffA, double *coeffB, bool fixA, bool fixB)
{
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
JKQTPYGraphErrorData* ge=dynamic_cast<JKQTPYGraphErrorData*>(datagraph);
JKQTPASSERT_M(ge!=nullptr, "datagraph needs to be convertible to JKQTPYGraphErrorData with a dynamic_cast!");
return jkqtpstatAddLinearWeightedRegression(plt, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()), ds->begin(ge->getYErrorColumn()), ds->end(ge->getYErrorColumn()), coeffA, coeffB, fixA, fixB, &jkqtp_inversePropSaveDefault<double>);
}
JKQTPXFunctionLineGraph *jkqtpstatAddRobustIRLSLinearRegression(JKQTPXYGraph *datagraph, double *coeffA, double *coeffB, bool fixA, bool fixB, double p, int iterations)
{
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
return jkqtpstatAddRobustIRLSLinearRegression(plt, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()), coeffA, coeffB, fixA, fixB, p, iterations);
}
JKQTPXFunctionLineGraph *jkqtpstatAddRegression(JKQTPXYGraph *datagraph, JKQTPStatRegressionModelType type, double *coeffA, double *coeffB, bool fixA, bool fixB) {
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
return jkqtpstatAddRegression(plt, type, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()), coeffA, coeffB, fixA, fixB);
}
JKQTPXFunctionLineGraph *jkqtpstatAddWeightedRegression(JKQTPXYGraph *datagraph, JKQTPStatRegressionModelType type, double *coeffA, double *coeffB, bool fixA, bool fixB)
{
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
JKQTPYGraphErrorData* ge=dynamic_cast<JKQTPYGraphErrorData*>(datagraph);
JKQTPASSERT_M(ge!=nullptr, "datagraph needs to be convertible to JKQTPYGraphErrorData with a dynamic_cast!");
return jkqtpstatAddWeightedRegression(plt, type, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()), ds->begin(ge->getYErrorColumn()), ds->end(ge->getYErrorColumn()), coeffA, coeffB, fixA, fixB, &jkqtp_inversePropSaveDefault<double>);
}
JKQTPXFunctionLineGraph *jkqtpstatAddRobustIRLSRegression(JKQTPXYGraph *datagraph, JKQTPStatRegressionModelType type, double *coeffA, double *coeffB, bool fixA, bool fixB, double p, int iterations)
{
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
return jkqtpstatAddRobustIRLSRegression(plt, type, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()), coeffA, coeffB, fixA, fixB, p, iterations);
}
JKQTPXFunctionLineGraph *jkqtpstatAddPolyFit(JKQTPXYGraph *datagraph, size_t P) {
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
return jkqtpstatAddPolyFit(plt, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()),P);
}

View File

@ -20,10 +20,14 @@
#include "jkqtcommon/jkqtp_imexport.h"
#include "jkqtcommon/jkqtpstatisticstools.h"
#include "jkqtcommon/jkqtptoolsdebugging.h"
#include "jkqtplotter/jkqtpgraphsbase.h"
#include "jkqtplotter/jkqtpgraphsbaseerrors.h"
#include "jkqtplotter/jkqtpgraphsboxplot.h"
#include "jkqtplotter/jkqtpgraphsscatter.h"
#include "jkqtplotter/jkqtpgraphssinglecolumnsymbols.h"
#include "jkqtplotter/jkqtpgraphsbarchart.h"
#include "jkqtplotter/jkqtpgraphsevaluatedfunction.h"
#ifndef JKQTPGRAPHSSTATISTICSADAPTORS_H_INCLUDED
#define JKQTPGRAPHSSTATISTICSADAPTORS_H_INCLUDED
@ -48,6 +52,9 @@
jkqtpstatAddHBoxplot(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), -0.3);
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_boxplots_simple.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstat5NumberStatistics()
*/
template <class InputIt>
@ -143,6 +150,8 @@ inline JKQTPBoxplotVerticalElement* jkqtpstatAddVBoxplot(JKQTBasePlotter* plotte
0.05, 0.95 // Quantiles for the boxplot box whiskers' ends
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_boxplots_outliers.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstat5NumberStatistics()
*/
template <class InputIt>
@ -201,7 +210,7 @@ inline std::pair<JKQTPBoxplotHorizontalElement*,JKQTPSingleColumnSymbolsGraph*>
jkqtpstatAddVBoxplotAndOutliers(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), -0.3,
0.25, 0.75, // 1. and 3. Quartile for the boxplot box
0.05, 0.95 // Quantiles for the boxplot box whiskers' ends
\endcode
\endcode
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstat5NumberStatistics()
*/
@ -257,6 +266,8 @@ inline std::pair<JKQTPBoxplotVerticalElement*,JKQTPSingleColumnSymbolsGraph*> jk
jkqtpstatAddHHistogram1DAutoranged(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), 11);
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_hist.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatHistogram1DAutoranged(), JKQTPBarVerticalGraph
*/
template <class InputIt>
@ -293,6 +304,8 @@ inline JKQTPBarVerticalGraph* jkqtpstatAddHHistogram1DAutoranged(JKQTBasePlotter
jkqtpstatAddHHistogram1DAutoranged(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), 0.5);
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_hist.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatHistogram1DAutoranged(), JKQTPBarVerticalGraph
*/
template <class InputIt>
@ -331,6 +344,8 @@ inline JKQTPBarVerticalGraph* jkqtpstatAddHHistogram1DAutoranged(JKQTBasePlotter
jkqtpstatAddHHistogram1D(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), bins.begin(), bins.end());
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_hist.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatHistogram1D(), JKQTPBarVerticalGraph
*/
template <class InputIt, class BinsInputIt>
@ -486,6 +501,9 @@ inline JKQTPBarHorizontalGraph* jkqtpstatAddVHistogram1D(JKQTBasePlotter* plotte
jkqtpstatAddHKDE1DAutoranged(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), 200);
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_kde.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatKDE1DAutoranged(), JKQTPXYLineGraph
*/
template <class InputIt>
@ -525,6 +543,9 @@ inline JKQTPXYLineGraph* jkqtpstatAddHKDE1DAutoranged(JKQTBasePlotter* plotter,
jkqtpstatAddHKDE1DAutoranged(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), 0.01);
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_kde.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatKDE1DAutoranged(), JKQTPXYLineGraph
*/
template <class InputIt>
@ -566,6 +587,9 @@ inline JKQTPXYLineGraph* jkqtpstatAddHKDE1DAutoranged(JKQTBasePlotter* plotter,
jkqtpstatAddHKDE1D(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), bins.begin(), bins.end());
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_kde.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatKDE1D(), JKQTPXYLineGraph
*/
template <class InputIt, class BinsInputIt>
@ -606,6 +630,8 @@ inline JKQTPXYLineGraph* jkqtpstatAddHKDE1D(JKQTBasePlotter* plotter, InputIt fi
jkqtpstatAddHKDE1D(plot1->getPlotter(), datastore1->begin(randomdatacol1), datastore1->end(randomdatacol1), bins.begin(), bins.end());
\endcode
\image html jkqtplotter_simpletest_datastore_statistics_kde.png
\see \ref JKQTPlotterBasicJKQTPDatastoreStatistics, jkqtpstatKDE1D(), JKQTPXYLineGraph
*/
template <class InputIt>
@ -625,6 +651,617 @@ inline JKQTPXYLineGraph* jkqtpstatAddHKDE1D(JKQTBasePlotter* plotter, InputIt fi
/*! \brief calculate the linear regression coefficients for a given data range \a firstX / \a firstY ... \a lastX / \a lastY where the model is \f$ f(x)=a+b\cdot x \f$
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param plotter the plotter to which to add the resulting graph
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
Example:
\code
jkqtpstatAddLinearRegression(plot1->getPlotter(), datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY));
\endcode
\image html jkqtplotter_simpletest_datastore_regression_lin.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatLinearRegression()
*/
template <class InputItX, class InputItY>
inline JKQTPXFunctionLineGraph* jkqtpstatAddLinearRegression(JKQTBasePlotter* plotter, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false) {
double cA;
if (coeffA) cA=*coeffA;
double cB;
if (coeffB) cB=*coeffB;
JKQTPASSERT_M(!fixA || (fixA && coeffA!=nullptr), "if fixA=true, coeffA needs to be provided");
JKQTPASSERT_M(!fixB || (fixB && coeffB!=nullptr), "if fixB=true, coeffB needs to be provided");
jkqtpstatLinearRegression(firstX, lastX, firstY, lastY, cA, cB, fixA, fixB);
JKQTPXFunctionLineGraph* g=new JKQTPXFunctionLineGraph(plotter);
g->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
g->setParamsV(cA, cB);
g->setTitle(QString("regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(cA)).arg(jkqtp_floattolatexqstr(cB)));
plotter->addGraph(g);
if (coeffA) *coeffA=cA;
if (coeffB) *coeffB=cB;
return g;
}
/*! \brief calculate the linear regression coefficients for a given data data used to draw any `JKQTPXYGraph` \a datagraph where the model is \f$ f(x)=a+b\cdot x \f$
\ingroup jkqtptools_math_statistics_adaptors
\param datagraph graph representing the (x,y) datapairs to which to fit the regression line
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
Example:
\code
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
jkqtpstatAddLinearRegression(graphD);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_lin.png
\note The line graph is added to the same plotter that is the parent of \a datagraph !
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatLinearRegression()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph* jkqtpstatAddLinearRegression(JKQTPXYGraph *datagraph, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false);
/*! \brief calculate the (robust) iteratively reweighted least-squares (IRLS) estimate for the parameters of the model \f$ f(x)=a+b\cdot x \f$
for a given data range \a firstX / \a firstY ... \a lastX / \a lastY
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param plotter the plotter to which to add the resulting graph
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
\param p regularization parameter, the optimization problem is formulated in the \f$ L_p \f$ norm, using this \a p (see image below for an example)
\param iterations the number of iterations the IRLS algorithm performs
Example:
\code
jkqtpstatAddRobustIRLSLinearRegression(plot1->getPlotter(), datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY));
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linrobust.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatRobustIRLSLinearRegression()
*/
template <class InputItX, class InputItY>
inline JKQTPXFunctionLineGraph* jkqtpstatAddRobustIRLSLinearRegression(JKQTBasePlotter* plotter, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false, double p=1.1, int iterations=100) {
double cA;
if (coeffA) cA=*coeffA;
double cB;
if (coeffB) cB=*coeffB;
JKQTPASSERT_M(!fixA || (fixA && coeffA!=nullptr), "if fixA=true, coeffA needs to be provided");
JKQTPASSERT_M(!fixB || (fixB && coeffB!=nullptr), "if fixB=true, coeffB needs to be provided");
jkqtpstatRobustIRLSLinearRegression(firstX, lastX, firstY, lastY, cA, cB, fixA, fixB, p, iterations);
JKQTPXFunctionLineGraph* g=new JKQTPXFunctionLineGraph(plotter);
g->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
g->setParamsV(cA, cB);
g->setTitle(QString("robust regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(cA)).arg(jkqtp_floattolatexqstr(cB)));
plotter->addGraph(g);
if (coeffA) *coeffA=cA;
if (coeffB) *coeffB=cB;
return g;
}
/*! \brief calculate the (robust) iteratively reweighted least-squares (IRLS) estimate for the parameters of the model \f$ f(x)=a+b\cdot x \f$
for a given data range \a firstX / \a firstY ... \a lastX / \a lastY
\ingroup jkqtptools_math_statistics_adaptors
\param datagraph graph representing the (x,y) datapairs to which to fit the regression line
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
\param p regularization parameter, the optimization problem is formulated in the \f$ L_p \f$ norm, using this \a p (see image below for an example)
\param iterations the number of iterations the IRLS algorithm performs
Example:
\code
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
jkqtpstatAddRobustIRLSLinearRegression(graphD);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linrobust.png
\note The line graph is added to the same plotter that is the parent of \a datagraph !
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatRobustIRLSLinearRegression()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph *jkqtpstatAddRobustIRLSLinearRegression(JKQTPXYGraph *datagraph, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false, double p=1.1, int iterations=100);
/*! \brief calculate the weighted linear regression coefficients for a given for a given data range \a firstX / \a firstY / \a firstW ... \a lastX / \a lastY / \a lastW where the model is \f$ f(x)=a+b\cdot x \f$
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\tparam InputItW standard iterator type of \a firstW and \a lastW.
\param plotter the plotter to which to add the resulting graph
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param firstW iterator pointing to the first item in the weight-dataset to use \f$ w_1 \f$
\param lastW iterator pointing behind the last item in the weight-dataset to use \f$ w_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
\param fWeightDataToWi an optional function, which is applied to the data from \a firstW ... \a lastW to convert them to weight, i.e. \c wi=fWeightDataToWi(*itW)
e.g. if you use data used to draw error bars, you can use jkqtp_inversePropSaveDefault(). The default is jkqtp_identity(), which just returns the values.
In the case of jkqtp_inversePropSaveDefault(), a datapoint x,y, has a large weight, if it's error is small and in the case if jkqtp_identity() it's weight
is directly proportional to the given value.
Example:
\code
double coeffA=0, coeffB=0;
jkqtpstatLinearWeightedRegression(datastore1->begin(colWLinX), datastore1->end(colWLinX),
datastore1->begin(colWLinY), datastore1->end(colWLinY),
datastore1->begin(colWLinE), datastore1->end(colWLinE),
coeffA, coeffB, false, false,
&jkqtp_inversePropSaveDefault<double>);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linweight.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatLinearRegression()
*/
template <class InputItX, class InputItY, class InputItW>
inline JKQTPXFunctionLineGraph* jkqtpstatAddLinearWeightedRegression(JKQTBasePlotter* plotter, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, InputItW firstW, InputItW lastW, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false, std::function<double(double)> fWeightDataToWi=&jkqtp_identity<double>) {
double cA;
if (coeffA) cA=*coeffA;
double cB;
if (coeffB) cB=*coeffB;
JKQTPASSERT_M(!fixA || (fixA && coeffA!=nullptr), "if fixA=true, coeffA needs to be provided");
JKQTPASSERT_M(!fixB || (fixB && coeffB!=nullptr), "if fixB=true, coeffB needs to be provided");
jkqtpstatLinearWeightedRegression(firstX, lastX, firstY, lastY, firstW, lastW, cA, cB, fixA, fixB, fWeightDataToWi);
JKQTPXFunctionLineGraph* g=new JKQTPXFunctionLineGraph(plotter);
g->setSpecialFunction(JKQTPXFunctionLineGraph::SpecialFunction::Line);
g->setParamsV(cA, cB);
g->setTitle(QString("weighted regression: $f(x) = %1 + %2 \\cdot x$").arg(jkqtp_floattolatexqstr(cA)).arg(jkqtp_floattolatexqstr(cB)));
plotter->addGraph(g);
if (coeffA) *coeffA=cA;
if (coeffB) *coeffB=cB;
return g;
}
/*! \brief calculate the linear weighted regression coefficients for a given data data used to draw any `JKQTPXYGraph` \a datagraph , which also implements JKQTPYGraphErrorData and where the model is \f$ f(x)=a+b\cdot x \f$
\ingroup jkqtptools_math_statistics_adaptors
\param datagraph graph representing the (x,y,error) data triples to which to fit the regression line
The errors are used as iverse weights!
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
Example:
\code
JKQTPXYLineErrorGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineErrorGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
graphD->setYErrorColumn(static_cast<int>(colWLinE));
jkqtpstatAddLinearWeightedRegression(graphD);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linweight.png
\note The line graph is added to the same plotter that is the parent of \a datagraph !
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatLinearRegression()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph* jkqtpstatAddLinearWeightedRegression(JKQTPXYGraph *datagraph, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false);
/*! \brief calculate the linear regression coefficients for a given data range \a firstX / \a firstY ... \a lastX / \a lastY where the model is defined by \a type
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param plotter the plotter to which to add the resulting graph
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
Example:
\code
jkqtpstatRegression(plot1->getPlotter(), JKQTPStatRegressionModelType::Exponential, datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY));
\endcode
\image html jkqtplotter_simpletest_datastore_regression_nonlinreg.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatRegression()
*/
template <class InputItX, class InputItY>
inline JKQTPXFunctionLineGraph* jkqtpstatAddRegression(JKQTBasePlotter* plotter, JKQTPStatRegressionModelType type, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false) {
double cA;
if (coeffA) cA=*coeffA;
double cB;
if (coeffB) cB=*coeffB;
JKQTPASSERT_M(!fixA || (fixA && coeffA!=nullptr), "if fixA=true, coeffA needs to be provided");
JKQTPASSERT_M(!fixB || (fixB && coeffB!=nullptr), "if fixB=true, coeffB needs to be provided");
jkqtpstatRegression(type, firstX, lastX, firstY, lastY, cA, cB, fixA, fixB);
JKQTPXFunctionLineGraph* g=new JKQTPXFunctionLineGraph(plotter);
g->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(type, cA, cB));
g->setTitle(QString("regression: $%1$").arg(jkqtpstatRegressionModel2Latex(type, cA, cB)));
plotter->addGraph(g);
if (coeffA) *coeffA=cA;
if (coeffB) *coeffB=cB;
return g;
}
/*! \brief calculate the linear regression coefficients for a given data data used to draw any `JKQTPXYGraph` \a datagraph where the model is defined by \a type
\param datagraph graph representing the (x,y) datapairs to which to fit the regression line
\param type model to be fitted
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
Example:
\code
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
jkqtpstatRegression(graphD, JKQTPStatRegressionModelType::Exponential);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_nonlinreg.png
\note The line graph is added to the same plotter that is the parent of \a datagraph !
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatRegression()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph* jkqtpstatAddRegression(JKQTPXYGraph *datagraph, JKQTPStatRegressionModelType type, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false);
/*! \brief calculate the (robust) iteratively reweighted least-squares (IRLS) estimate for the parameters where the model is defined by \a type
for a given data range \a firstX / \a firstY ... \a lastX / \a lastY
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param plotter the plotter to which to add the resulting graph
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
\param p regularization parameter, the optimization problem is formulated in the \f$ L_p \f$ norm, using this \a p (see image below for an example)
\param iterations the number of iterations the IRLS algorithm performs
Example:
\code
jkqtpstatAddRobustIRLSRegression(plot1->getPlotter(), JKQTPStatRegressionModelType::Exponential, datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY));
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linrobust.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatRobustIRLSRegression()
*/
template <class InputItX, class InputItY>
inline JKQTPXFunctionLineGraph* jkqtpstatAddRobustIRLSRegression(JKQTBasePlotter* plotter, JKQTPStatRegressionModelType type, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false, double p=1.1, int iterations=100) {
double cA;
if (coeffA) cA=*coeffA;
double cB;
if (coeffB) cB=*coeffB;
JKQTPASSERT_M(!fixA || (fixA && coeffA!=nullptr), "if fixA=true, coeffA needs to be provided");
JKQTPASSERT_M(!fixB || (fixB && coeffB!=nullptr), "if fixB=true, coeffB needs to be provided");
jkqtpstatRobustIRLSRegression(type, firstX, lastX, firstY, lastY, cA, cB, fixA, fixB, p, iterations);
JKQTPXFunctionLineGraph* g=new JKQTPXFunctionLineGraph(plotter);
g->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(type, cA, cB));
g->setTitle(QString("robust regression: $%1$").arg(jkqtpstatRegressionModel2Latex(type, cA, cB)));
plotter->addGraph(g);
if (coeffA) *coeffA=cA;
if (coeffB) *coeffB=cB;
return g;
}
/*! \brief calculate the (robust) iteratively reweighted least-squares (IRLS) estimate for the parameters where the model is defined by \a type
for a given data range \a firstX / \a firstY ... \a lastX / \a lastY
\ingroup jkqtptools_math_statistics_adaptors
\param datagraph graph representing the (x,y) datapairs to which to fit the regression line
\param type model to be fitted
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
\param p regularization parameter, the optimization problem is formulated in the \f$ L_p \f$ norm, using this \a p (see image below for an example)
\param iterations the number of iterations the IRLS algorithm performs
Example:
\code
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
jkqtpstatAddRobustIRLSRegression(graphD, JKQTPStatRegressionModelType::Exponential);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linrobust.png
\note The line graph is added to the same plotter that is the parent of \a datagraph !
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatRobustIRLSRegression()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph *jkqtpstatAddRobustIRLSRegression(JKQTPXYGraph *datagraph, JKQTPStatRegressionModelType type, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false, double p=1.1, int iterations=100);
/*! \brief calculate the weighted linear regression coefficients for a given for a given data range \a firstX / \a firstY / \a firstW ... \a lastX / \a lastY / \a lastW where the model is defined by \a type
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\tparam InputItW standard iterator type of \a firstW and \a lastW.
\param plotter the plotter to which to add the resulting graph
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param firstW iterator pointing to the first item in the weight-dataset to use \f$ w_1 \f$
\param lastW iterator pointing behind the last item in the weight-dataset to use \f$ w_N \f$
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
\param fWeightDataToWi an optional function, which is applied to the data from \a firstW ... \a lastW to convert them to weight, i.e. \c wi=fWeightDataToWi(*itW)
e.g. if you use data used to draw error bars, you can use jkqtp_inversePropSaveDefault(). The default is jkqtp_identity(), which just returns the values.
In the case of jkqtp_inversePropSaveDefault(), a datapoint x,y, has a large weight, if it's error is small and in the case if jkqtp_identity() it's weight
is directly proportional to the given value.
Example:
\code
double coeffA=0, coeffB=0;
jkqtpstatAddWeightedRegression(plotter, JKQTPStatRegressionModelType::Exponential,
datastore1->begin(colWLinX), datastore1->end(colWLinX),
datastore1->begin(colWLinY), datastore1->end(colWLinY),
datastore1->begin(colWLinE), datastore1->end(colWLinE),
coeffA, coeffB, false, false,
&jkqtp_inversePropSaveDefault<double>);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linweight.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatWeightedRegression()
*/
template <class InputItX, class InputItY, class InputItW>
inline JKQTPXFunctionLineGraph* jkqtpstatAddWeightedRegression(JKQTBasePlotter* plotter, JKQTPStatRegressionModelType type, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, InputItW firstW, InputItW lastW, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false, std::function<double(double)> fWeightDataToWi=&jkqtp_identity<double>) {
double cA;
if (coeffA) cA=*coeffA;
double cB;
if (coeffB) cB=*coeffB;
JKQTPASSERT_M(!fixA || (fixA && coeffA!=nullptr), "if fixA=true, coeffA needs to be provided");
JKQTPASSERT_M(!fixB || (fixB && coeffB!=nullptr), "if fixB=true, coeffB needs to be provided");
jkqtpstatWeightedRegression(type, firstX, lastX, firstY, lastY, firstW, lastW, cA, cB, fixA, fixB, fWeightDataToWi);
JKQTPXFunctionLineGraph* g=new JKQTPXFunctionLineGraph(plotter);
g->setPlotFunctionFunctor(jkqtpStatGenerateRegressionModel(type, cA, cB));
g->setTitle(QString("weighted regression: $%1$").arg(jkqtpstatRegressionModel2Latex(type, cA, cB)));
plotter->addGraph(g);
if (coeffA) *coeffA=cA;
if (coeffB) *coeffB=cB;
return g;
}
/*! \brief calculate the linear weighted regression coefficients for a given data data used to draw any `JKQTPXYGraph` \a datagraph , which also implements JKQTPYGraphErrorData and where the model is defined by \a type
\ingroup jkqtptools_math_statistics_adaptors
\param datagraph graph representing the (x,y,error) data triples to which to fit the regression line
The errors are used as iverse weights!
\param type model to be fitted
\param[in,out] coeffA returns the offset of the linear model
\param[in,out] coeffB returns the slope of the linear model
\param fixA if \c true, the offset coefficient \f$ a \f$ is not determined by the fit, but the value provided in \a coeffA is used \note If \a fixA \c ==true, You need to provide a value for A in \a coeffA
\param fixB if \c true, the slope coefficient \f$ b \f$ is not determined by the fit, but the value provided in \a coeffB is used \note If \a fixB \c ==true, You need to provide a value for B in \a coeffB
Example:
\code
JKQTPXYLineErrorGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineErrorGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
graphD->setYErrorColumn(static_cast<int>(colWLinE));
jkqtpstatAddWeightedRegression(graphD, JKQTPStatRegressionModelType::Exponential);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_linweight.png
\note The line graph is added to the same plotter that is the parent of \a datagraph !
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatWeightedRegression()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph* jkqtpstatAddWeightedRegression(JKQTPXYGraph *datagraph, JKQTPStatRegressionModelType type, double* coeffA=nullptr, double* coeffB=nullptr, bool fixA=false, bool fixB=false);
/*! \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$ from a given data range \a firstX / \a firstY ... \a lastX / \a lastY
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\tparam OutputItP output iterator for the polynomial coefficients
\param plotter the plotter to which to add the resulting graph
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param P degree of the polynomial (P>=N !!!)
\param[out] firstRes Iterator (of type \a OutputItP ), which receives the (P+1)-entry vector with the polynomial coefficients \f$ p_i \f$
Example:
\code
std::vector<double> pFit;
jkqtpstatAddPolyFit(plot1->getPlotter(), JKQTPStatRegressionModelType::Exponential, datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY), 3, std::back_inserter(pFit));
\endcode
\image html jkqtplotter_simpletest_datastore_regression_polynom.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatPolyFit()
*/
template <class InputItX, class InputItY, class OutputItP>
inline JKQTPXFunctionLineGraph* jkqtpstatAddPolyFit(JKQTBasePlotter* plotter, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, size_t P, OutputItP firstRes) {
std::vector<double> pFit;
JKQTPXFunctionLineGraph* gPoly=new JKQTPXFunctionLineGraph(plotter);
jkqtpstatPolyFit(firstX,lastX,firstY,lastY,P,std::back_inserter(pFit));
gPoly->setPlotFunctionFunctor(jkqtpstatGeneratePolynomialModel(pFit.begin(), pFit.end()));
gPoly->setTitle(QString("regression: $%1$").arg(jkqtpstatPolynomialModel2Latex(pFit.begin(), pFit.end())));
std::copy(pFit.begin(), pFit.end(), firstRes);
plotter->addGraph(gPoly);
return gPoly;
}
/*! \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$ from a given data range \a firstX / \a firstY ... \a lastX / \a lastY
\ingroup jkqtptools_math_statistics_adaptors
\tparam InputItX standard iterator type of \a firstX and \a lastX.
\tparam InputItY standard iterator type of \a firstY and \a lastY.
\param plotter the plotter to which to add the resulting graph
\param type model to be fitted
\param firstX iterator pointing to the first item in the x-dataset to use \f$ x_1 \f$
\param lastX iterator pointing behind the last item in the x-dataset to use \f$ x_N \f$
\param firstY iterator pointing to the first item in the y-dataset to use \f$ y_1 \f$
\param lastY iterator pointing behind the last item in the y-dataset to use \f$ y_N \f$
\param P degree of the polynomial (P>=N !!!)
Example:
\code
jkqtpstatAddPolyFit(plot1->getPlotter(), JKQTPStatRegressionModelType::Exponential, datastore1->begin(colLinX), datastore1->end(colLinX), datastore1->begin(colLinY), datastore1->end(colLinY), 3);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_polynom.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatPolyFit()
*/
template <class InputItX, class InputItY>
inline JKQTPXFunctionLineGraph* jkqtpstatAddPolyFit(JKQTBasePlotter* plotter, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, size_t P) {
std::vector<double> pFit;
return jkqtpstatAddPolyFit(plotter, firstX,lastX,firstY,lastY,P,std::back_inserter(pFit));
}
/*! \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$ from a given JKQTPXYGraph \a datagraph
\ingroup jkqtptools_math_statistics_adaptors
\tparam OutputItP output iterator for the polynomial coefficients
\param datagraph graph representing the (x,y) datapairs to which to fit the regression line
\param P degree of the polynomial (P>=N !!!)
\param[out] firstRes Iterator (of type \a OutputItP ), which receives the (P+1)-entry vector with the polynomial coefficients \f$ p_i \f$
Example:
\code
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
std::vector<double> pFit;
jkqtpstatAddPolyFit(graphD, 3,std::back_inserter(pFit));
\endcode
\image html jkqtplotter_simpletest_datastore_regression_polynom.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatPolyFit()
*/
template <class OutputItP>
inline JKQTPXFunctionLineGraph* jkqtpstatAddPolyFit(JKQTPXYGraph *datagraph, size_t P, OutputItP firstRes) {
JKQTBasePlotter* plt=datagraph->getParent();
JKQTPDatastore* ds=plt->getDatastore();
return jkqtpstatAddPolyFit(plt, ds->begin(datagraph->getXColumn()), ds->end(datagraph->getXColumn()), ds->begin(datagraph->getYColumn()), ds->end(datagraph->getYColumn()),P,firstRes);
}
/*! \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$ from a given JKQTPXYGraph \a datagraph
\ingroup jkqtptools_math_statistics_adaptors
\param datagraph graph representing the (x,y) datapairs to which to fit the regression line
\param P degree of the polynomial (P>=N !!!)
Example:
\code
JKQTPXYLineGraph* graphD;
plot1->addGraph(graphD=new JKQTPXYLineGraph(plot1));
graphD->setXYColumns(colLinX, colLinY);
jkqtpstatAddPolyFit(graphD, 3);
\endcode
\image html jkqtplotter_simpletest_datastore_regression_polynom.png
\see \ref JKQTPlotterBasicJKQTPDatastoreRegression, jkqtpstatPolyFit()
*/
JKQTP_LIB_EXPORT JKQTPXFunctionLineGraph* jkqtpstatAddPolyFit(JKQTPXYGraph *datagraph, size_t P);
#endif // JKQTPGRAPHSSTATISTICSADAPTORS_H_INCLUDED

Binary file not shown.

After

Width:  |  Height:  |  Size: 137 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 23 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 32 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 24 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 33 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 33 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 36 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB