gsl - 使用 GSL C++ 的 Logistic S 曲线的雅可比矩阵

Question

我尝试使用雅可比矩阵来解决逻辑参数 y=k/(1+exp(-rt+b)) 用于使用 GSL c++ 库的非线性回归，但没有成功，可能是三个变量 k、r 和 b 的相应导数错了，下面是我的 C++ 代码片段。如果我设置 fdf.df=NULL; 我可以获得与已知 Fortran 求解器进行比较的未优化参数，k 值似乎未优化。有人可以建议如何正确编码 int expb_df 下的雅可比矩阵吗？

#include<gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include <gsl/gsl_randist.h>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_multifit_nlinear.h>



vector<int> vec_indata = {3, 1, 0, 0, 3, 1, 0, 0, 0, 0, 2, 2, 0, 0, 4, 0, 2, 0, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
                        2, 4, 3, 4, 14, 5, 28, 10, 6, 18, 12, 20, 9, 39, 41, 190, 125, 120, 117, 110, 130, 153, 123, 212, 106,
                        172, 235, 130, 159, 150, 156, 140, 142, 208, 217, 150, 179, 131, 170, 156, 109, 118, 184, 153, 134, 170,
                        85, 110, 69, 54, 84, 36, 57, 50, 71, 88, 51, 38, 40, 31, 94, 57, 69, 105, 122, 55, 30, 45, 39, 68, 54,
                        67, 70, 16, 37, 40, 36, 17, 22, 47, 37, 31, 50, 78, 48, 60, 172, 187, 15, 10, 103, 30, 57, 38, 20, 93,
                        277, 19, 37, 19, 7, 7, 2, 31, 33, 43, 8, 41, 11, 10, 14, 6, 21, 16, 15, 3, 6, 4, 6, 10, 18, 3, 2, 1, 3,
                         5, 10, 5, 5, 6,3,6,13,8,14,7,4,5,3,18,9,15,21,15,16,9,21,23,13,7,39,13,8,12,9,14,2,1,21,15,25,7,13,
                         11,9,11,15,20,26,25,12,7,16,5,9,8,10,7,11,6,5,10,11,17,6,14,6,14,11,6,6,62,100,24};
 


int  expb_f(const gsl_vector *x, void *indata, gsl_vector *f) {

    size_t n = ((struct idata*)indata)->n;
    double *t = ((struct idata*)indata)->t;
    double *y = ((struct idata*)indata)->y;
    double k = gsl_vector_get(x, 0);
    double r = gsl_vector_get(x, 1);
    double b = gsl_vector_get(x, 2);

    size_t i;

    for (int i = 0; i < n; i++) {

        /* logistic Model y=k/(1+exp(-r*t+b)) */
        double y_hat = k / (1 + std::exp(-r * t[i] + b));
        gsl_vector_set(f, i, y_hat - y[i]);

    }

    return GSL_SUCCESS;

}



int  expb_df(const gsl_vector * x, void *data, gsl_matrix * J) {

    size_t n = ((struct idata*)data)->n;
    double *t = ((struct idata*)data)->t;
    double *y = ((struct idata*)data)->y;
    double k = gsl_vector_get(x, 0);
    double r = gsl_vector_get(x, 1);
    double b = gsl_vector_get(x, 2);
    size_t i;
    for (int i = 0; i < n; i++) {

        double dv = -1 * std::pow((1 + std::exp(-r * t[i])), -2);
        double e = std::exp(-r * t[i] + b);
        gsl_matrix_set(J, i, 0, 1);
        gsl_matrix_set(J, i, 1, -t[i]*e);
        gsl_matrix_set(J, i, 2, -e);

    }

    return GSL_SUCCESS;

}


void callback(const size_t iter, void *params,
    const gsl_multifit_nlinear_workspace *w) {
    String out_txt;
    gsl_vector *f = gsl_multifit_nlinear_residual(w);
    gsl_vector *x = gsl_multifit_nlinear_position(w);
    double rcond;
    /* compute reciprocal condition number of J(x) */
    // gsl_multifit_nlinear_rcond(&rcond, w);
    out_txt = "Iter " + String(iter) + " k=" + String(gsl_vector_get(x, 0)) +
        " r=" + String(gsl_vector_get(x, 1)) + " b=" +
        String(gsl_vector_get(x, 2)) + "   ," + String(gsl_blas_dnrm2(f)) +
        "\r\n" + out_txt;
    
}

void FitLogistic() {
    const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
    gsl_multifit_nlinear_workspace *w;
    gsl_multifit_nlinear_fdf fdf;
    gsl_multifit_nlinear_parameters fdf_params =
        gsl_multifit_nlinear_default_parameters();
    size_t N = size(vec_indata);
    const size_t n = N;
    const size_t p = 3;
    gsl_vector *f;
    gsl_matrix *J;
    gsl_matrix *covar = gsl_matrix_alloc(p, p);
    double t[N], y[N], weights[N];

    struct idata d = {
        n, t, y
    };


    double x_init[3] = {20000, 0.01, 1};
    gsl_vector_view x = gsl_vector_view_array(x_init, p);
    gsl_vector_view wts = gsl_vector_view_array(weights, n);
    double chisq, chisq0;
    int status, info;
    size_t i;
    const double xtol = 1e-8;
    const double gtol = 1e-8;
    const double ftol = 0.0;

    /* define the function to be minimized */
    fdf.f = expb_f;
   

// jacobian matrix 
 fdf.df=NULL;
    //fdf.df = expb_df;
    fdf.fvv = NULL;
    fdf.n = n;
    fdf.p = p;
    fdf.params = &d;
    int cum_actual_cnt = 0;
    // ShowMessage(String(N));

    for (i = 0; i < N; i++) {
        // t[i] = i*N/(n-1);
        t[i] = i;
        y[i] = vec_indata.at(i) + cum_actual_cnt;
        cum_actual_cnt += vec_indata.at(i);
        weights[i] = 1 / ((0.1 * y[i]) * (0.1 * y[i]));
    }

    // ShowMessage(String(y[N-1]));
    /* allocate workspace with default parameters */
    w = gsl_multifit_nlinear_alloc(T, &fdf_params, n, p);
    /* initialize solver with starting point and weights */
    gsl_multifit_nlinear_init(&x.vector, &fdf, w);
    // gsl_multifit_nlinear_winit(&x.vector, &wts.vector, &fdf, w);
    /* solve the system with a maximum of 1000 iterations */
    status = gsl_multifit_nlinear_driver(1000, xtol, gtol, ftol, callback, NULL,
        &info, w);
    ShowMessage("Fit Status " + String(gsl_strerror(status)) + " " +
        gsl_multifit_nlinear_name(w) + "/" + gsl_multifit_nlinear_trs_name(w));
    /* compute covariance of best fit parameters */
    J = gsl_multifit_nlinear_jac(w);
    gsl_multifit_nlinear_covar(J, 0.0, covar);
    gsl_multifit_nlinear_free(w);

}

标签： gsl

Answer 1

在针对 k、r 和 b 调整偏导数的误差后，我设法获得了与上述代码完全相同的输出结果，该代码在没有分配雅可比矩阵或放置 fdf.df=NULL 的情况下进行编码；我无法确定解决方案是否最佳，也许有人能够发表评论，谢谢

 fdf.df = expb_df;

for (int i = 0; i < n; i++) {

        double dv = std::pow((1 + std::exp(-r * t[i] + b)), 2);
        double e = std::exp(-r * t[i] + b);
            double u =1+ std::exp(-r * t[i] + b);
        sigma = i - 1;
        gsl_matrix_set(J, i, 0, 1/u);
        gsl_matrix_set(J, i, 1, (k*t[i]*e)/(u*u));
        gsl_matrix_set(J, i, 2, (-k*e)/(u*u));

    }