使用 GSL C++ 的逻辑 S 曲线的雅可比矩阵
Jacobian matrix for Logistic S-Curve using GSL C++
我尝试使用雅可比矩阵求解逻辑参数 y=k/(1+exp(-rt+b)) 以使用 GSL c++ 库进行非线性回归，但没有成功，可能是三个变量 k 的相应导数,r 和 b 是错误的，下面是我的 C++ 代码片段。
如果我设置 fdf.df=NULL;
我可以获得与已知 Fortran 求解器进行比较的未优化参数，k 值似乎未优化。有人可以建议如何正确编码下的雅可比矩阵
请输入 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);
}
在调整了关于 k,r,b 的偏导数误差后，我设法获得了与上述代码完全相似的输出结果，该代码是在没有分配雅可比矩阵或将 fdf.df=空；我无法确定解决方案是否最佳，也许有人可以发表评论，谢谢
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));
}
我尝试使用雅可比矩阵求解逻辑参数 y=k/(1+exp(-rt+b)) 以使用 GSL c++ 库进行非线性回归，但没有成功，可能是三个变量 k 的相应导数,r 和 b 是错误的，下面是我的 C++ 代码片段。 如果我设置 fdf.df=NULL; 我可以获得与已知 Fortran 求解器进行比较的未优化参数，k 值似乎未优化。有人可以建议如何正确编码下的雅可比矩阵 请输入 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);
}
在调整了关于 k,r,b 的偏导数误差后，我设法获得了与上述代码完全相似的输出结果，该代码是在没有分配雅可比矩阵或将 fdf.df=空；我无法确定解决方案是否最佳，也许有人可以发表评论，谢谢
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));
}