logistic

logistic4p.cpp

// The contents of this file are in the public domain. See LICENSE_FOR_EXAMPLE_PROGRAMS.txt
/*

    This is an example illustrating the use the general purpose non-linear
    least squares optimization routines from the dlib C++ Library.

    This example program will demonstrate how these routines can be used for data fitting.
    In particular, we will generate a set of data and then use the least squares
    routines to infer the parameters of the model which generated the data.
*/


#include <dlib/optimization.h>
#include <iostream>
#include <vector>


using namespace std;
using namespace dlib;

// ----------------------------------------------------------------------------------------

// typedef matrix<double,1,1> input_vector;
typedef matrix<double,4,1> parameter_vector;

// ----------------------------------------------------------------------------------------

// We will use this function to generate data.  It represents a function of 2 variables
// and 4 parameters.   The least squares procedure will be used to infer the values of
// the 4 parameters based on a set of input/output pairs.
double model (
    const double input,
    const parameter_vector& params
)
{
    double denominal;
    double nominal;

    const double a = params(0);
    const double b = params(1);
    const double c = params(2);
    const double d = params(3);

    const double x = input;

    denominal = a - d;
    nominal = 1 + pow(x/c, b);

    return denominal/nominal + d;
}

// ----------------------------------------------------------------------------------------

// This function is the "residual" for a least squares problem.   It takes an input/output
// pair and compares it to the output of our model and returns the amount of error.  The idea
// is to find the set of parameters which makes the residual small on all the data pairs.
double residual (
    const std::pair<double, double>& data,
    const parameter_vector& params
)
{
    return model(data.first, params) - data.second;
}

// ----------------------------------------------------------------------------------------
void logistic4PL(
    const double xs[16],
    const double ys[16],
    const int iDots,
    const int iTrend,
    double ymodle[16],
    double params[]) {

    // Now let's generate a bunch of input/output pairs according to our model.
    std::vector<std::pair<double, double>> data_samples;
    for(int i=0;i<iDots;i++) {
        data_samples.push_back(make_pair(xs[i], ys[i]));
    }

    int antilog = 1;
    if(iTrend < 1) {
        antilog = -1;
    }

    // Now let's use the solve_least_squares_lm() routine to figure out what the
    // parameters are based on just the data_samples.
    double ymin = data_samples[0].second;
    double ymax = data_samples[iDots-1].second;
    double cinit = data_samples[iDots-1].first/2;
    parameter_vector pvec;
    pvec = ymin,
        antilog,
        cinit,
        ymax;
    cout << "Use Levenberg-Marquardt, approximate derivatives" << endl;
    // If we didn't create the residual_derivative function then we could
    // have used this method which numerically approximates the derivatives for you.
    // .be_verbose here to show the result of iteration.
    solve_least_squares_lm(objective_delta_stop_strategy(1e-9).be_verbose(),
                         residual,
                         derivative(residual),
                         data_samples,
                         pvec);

    for(int i=0;i<iDots;i++) {
        ymodle[i] = model(xs[i], pvec);
    }

    for(int r=0;r<4;r++) {
        params[r] = pvec(r);
    }
}


int main()
{
    try
    {
        double xs[16] = {30.0, 300.0, 1000.0, 3000.0, 15000.0, 30000.0};
        double ys[16] = {0.0171, 0.0699, 0.1738, 0.4831, 1.8927, 3.1621};
        int iDots = 6;
        double yfit[16] = {0};
        double params[] = {0};

        logistic4PL(xs, ys, iDots, 1, yfit, params);
        cout << "inferred parameters: ";
        for(int i=0; i<4; i++){
            cout << params[i] << "   ";
        }
        cout << endl;

        for(int i=0; i<iDots; i++){
            cout << "y = " << ys[i] << "; " << "yfit = " << yfit[i] << endl;
        }
    }
    catch (std::exception& e)
    {
        cout << e.what() << endl;
    }
}