使用 GNU 科学库进行曲线拟合

问题描述 投票:0回答:1

我正在尝试使用 GSL 将曲线拟合到一组点。我对 GSL 不太有经验,因此从 chatGPT 获得了帮助。但代码无法编译,过了一会儿我感觉就像在对门把手说话。我想知道有人可以帮助我。

我正在尝试拟合岭回归,但附加的约束是曲线参数之间的差异也最小化。

#include <iostream>
#include <vector>
#include <gsl/gsl_multifit.h>

// Ridge loss function with regularization term and parameter difference constraints
int ridge_loss_with_constraints(const gsl_vector* x, void* params, gsl_vector* f) {
    gsl_vector* c = (gsl_vector*)params; // Coefficients (a, b, c) at time t2
    gsl_vector* c_t1 = (gsl_vector*)params + 3; // Coefficients (a_t1, b_t1, c_t1) at time t1

    double a = gsl_vector_get(c, 0);
    double b = gsl_vector_get(c, 1);
    double c_val = gsl_vector_get(c, 2);

    double a_t1 = gsl_vector_get(c_t1, 0);
    double b_t1 = gsl_vector_get(c_t1, 1);
    double c_t1_val = gsl_vector_get(c_t1, 2);

    double lambda = 0.1; // Regularization parameter
    double regularization_term = lambda * (a * a + b * b + c_val * c_val);

    // Calculate squared differences between parameters at t1 and t2
    double diff_a = a - a_t1;
    double diff_b = b - b_t1;
    double diff_c = c_val - c_t1_val;
    double parameter_difference_term = diff_a * diff_a + diff_b * diff_b + diff_c * diff_c;

    // Calculate the fitted values using the current parameters
    for (size_t i = 0; i < x->size; ++i) {
        double xi = gsl_vector_get(x, i);
        double fitted_y = a * xi * xi + b * xi + c_val;
        gsl_vector_set(f, i, fitted_y + regularization_term + parameter_difference_term);
    }

    return GSL_SUCCESS;
}

int main() {
    // Sample data points (x, y)
    std::vector<double> x_data = {1.0, 2.0, 3.0, 4.0, 5.0};
    std::vector<double> y_data = {2.5, 3.7, 4.8, 5.9, 7.1};

    // Number of data points
    size_t n = x_data.size();

    // Define the dependent variable vector 'y'
    gsl_vector* y = gsl_vector_alloc(n);
    for (size_t i = 0; i < n; ++i) {
        gsl_vector_set(y, i, y_data[i]);
    }

    // Initialize the coefficient vectors for time t2 and t1
    gsl_vector* c = gsl_vector_alloc(3); // Coefficients (a, b, c) at time t2
    gsl_vector* c_t1 = gsl_vector_alloc(3); // Coefficients (a_t1, b_t1, c_t1) at time t1

    // Set coefficients at time t1 (replace with actual values)
    gsl_vector_set(c_t1, 0, 1/* a_t1 value */);
    gsl_vector_set(c_t1, 1, 1/* b_t1 value */);
    gsl_vector_set(c_t1, 2, 1/* c_t1 value */);

    // Combine the coefficient vectors for use in the ridge_loss_with_constraints function
    gsl_vector* combined_params = gsl_vector_alloc(6); // [a, b, c, a_t1, b_t1, c_t1]
    for (size_t i = 0; i < 3; ++i) {
        gsl_vector_set(combined_params, i, gsl_vector_get(c, i));
        gsl_vector_set(combined_params, i + 3, gsl_vector_get(c_t1, i));
    }

    // Initialize the fitting function
    gsl_multifit_function f;
    f.f = ridge_loss_with_constraints;
    f.n = n;
    f.p = 3;
    f.params = combined_params;

    // Perform the fit using the Levenberg-Marquardt optimization algorithm
    gsl_multifit_fdfsolver* solver = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, n, 3);
    gsl_multifit_fdfsolver_set(solver, &f, c);

    // Print the results
    std::cout << "Fitted coefficients with L2 regularization and constraints: Intercept = "
              << gsl_vector_get(c, 0) << ", Slope = " << gsl_vector_get(c, 1) << ", c-term = " << gsl_vector_get(c, 2) << std::endl;

    // Free GSL objects
    gsl_vector_free(y);
    gsl_vector_free(c);
    gsl_vector_free(c_t1);
    gsl_vector_free(combined_params);
    gsl_multifit_fdfsolver_free(solver);

    return 0;
}

我正在编译以下内容

g++ -o my_program ridge_regression3.cpp -lgsl -lgslcblas -lm -O3

我遇到的一些错误是:

ridge_regression3.cpp: In function ‘int main()’:
ridge_regression3.cpp:68:5: error: ‘gsl_multifit_function’ was not declared in this scope; did you mean ‘gsl_multifit_linear’?
   68 |     gsl_multifit_function f;
      |     ^~~~~~~~~~~~~~~~~~~~~
      |     gsl_multifit_linear
ridge_regression3.cpp:69:5: error: ‘f’ was not declared in this scope
   69 |     f.f = ridge_loss_with_constraints;
      |     ^
ridge_regression3.cpp:75:5: error: ‘gsl_multifit_fdfsolver’ was not declared in this scope; did you mean ‘gsl_multifit_robust’?
   75 |     gsl_multifit_fdfsolver* solver = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, n, 3);
      |     ^~~~~~~~~~~~~~~~~~~~~~
      |     gsl_multifit_robust
ridge_regression3.cpp:75:29: error: ‘solver’ was not declared in this scope
   75 |     gsl_multifit_fdfsolver* solver = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, n, 3);
      |                             ^~~~~~
ridge_regression3.cpp:75:67: error: ‘gsl_multifit_fdfsolver_lmsder’ was not declared in this scope
   75 |     gsl_multifit_fdfsolver* solver = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, n, 3);
      |                                                                   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ridge_regression3.cpp:75:38: error: ‘gsl_multifit_fdfsolver_alloc’ was not declared in this scope; did you mean ‘gsl_multifit_robust_alloc’?
   75 |     gsl_multifit_fdfsolver* solver = gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, n, 3);
      |                                      ^~~~~~~~~~~~~~~~~~~~~~~~~~~~
      |                                      gsl_multifit_robust_alloc
ridge_regression3.cpp:76:5: error: ‘gsl_multifit_fdfsolver_set’ was not declared in this scope; did you mean ‘gsl_multifit_robust_est’?
   76 |     gsl_multifit_fdfsolver_set(solver, &f, c);
      |     ^~~~~~~~~~~~~~~~~~~~~~~~~~
      |     gsl_multifit_robust_est
ridge_regression3.cpp:87:5: error: ‘gsl_multifit_fdfsolver_free’ was not declared in this scope; did you mean ‘gsl_multifit_robust_free’?
   87 |     gsl_multifit_fdfsolver_free(solver);
      |     ^~~~~~~~~~~~~~~~~~~~~~~~~~~
      |     gsl_multifit_robust_free

我尝试阅读文档here,但有很多我不明白。

提前致谢。

c++ curve-fitting gnu gsl
1个回答
0
投票

要编译代码片段:

  1. 包含gsl/gsl_multifit_nlin.h:

    #include <gsl/gsl_multifit_nlin.h>

  2. gsl_multifit_fdfsolver_set
    的第二个参数不匹配,所以将
    gsl_multifit_function f;
    更改为
    gsl_multifit_function_fdf f;

程序仍然存在段错误,但现在可以编译了。

© www.soinside.com 2019 - 2024. All rights reserved.