使用GSL C ++的Logistic S曲线的Jacobian矩阵

问题描述

我尝试使用Jacobian矩阵通过GSL c ++库求解非线性回归的逻辑参数y = k /（1 + exp（-rt + b）），但未成功，可能分别衍生了三个变量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,3,2,4,14,5,28,10,6,18,12,20,9,39,41,190,125,120,117,110,130,153,123,212,106,172,235,159,150,156,140,142,208,217,179,131,170,109,118,184,134,85,69,54,84,36,57,50,71,88,51,38,40,31,94,105,122,55,30,45,68,67,70,16,37,17,22,47,78,48,60,187,15,103,93,277,19,7,33,43,8,11,21,13,23,25,26,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,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,1);
        gsl_matrix_set(J,-t[i]*e);
        gsl_matrix_set(J,-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,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);

}

解决方法

在针对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,1/u);
        gsl_matrix_set(J,1,(k*t[i]*e)/(u*u));
        gsl_matrix_set(J,2,(-k*e)/(u*u));

    }