四参数拟合之C代码(基于GSL)

前言

前面已经通过3篇介绍了关于四参数拟合算法，主要用到拟牛顿和LMF方法，都是用matlab实现的。本想用C自己实现一遍，一想到其中涉及不少的矩阵运算，考虑到用C实现的工作量，立马打退堂鼓。几年前用过GSL，感觉调用它的线性代数运算还是挺方便的，于是果断去网上找GSL。目前，GSL已经更新到2.6，它真是一个宝库，里面居然已经包含了LMF算法。尽管在编码，让人觉得有些诡异，而且在windows下编译让人觉得心累，但是确实好用。花了一天多的时间终于将它调试通过。

过程

调试环境：Win10+VC6，是的，正是古老的VC6，不是不想用最新的IDE，实在是太大太慢，不想安装。
下载最新的GSL2.6，解压，将头文件包含到VC6的目录里面。新建vc工程，在GSL官网上找了LMF的调用例程，编译，问题一大堆，足以令人望而生畏。好在有过GSL老版本的移植经验，没有被吓倒。接下来本着缺啥补啥，错啥改啥。花了一天功夫，将例程调试通过，再花一天的功夫修改例程到L4P算法，大功告成。
LM算法可以用到许多无约束的问题当中去，只要提供拟合函数，雅可比矩阵生成函数和必要的设定参数即可。
四参数的B可以有正有负，其他为正数，真是啰嗦，其导数形式复杂，可以用matlab的diff函数，避免手动推导出错。省事起见，这里，只保证A>0,D>0
$$\begin{gathered} f(w_1;x_i)=D_1^2+\frac{A_1^2-D_1^2}{1+(x_i/C{)}^B} \\ {w}_1=[A_1,B,C,D_1] \end{gathered}$$
用L4P的起始点非常重要，起始点选的不好，很可能导致结果错误。这里就参照L4P的起始点设置算法

slope=(y(end)-y(1))/(x(end)-x(1));
if isempty(st)
    [~,Idx]=min(abs((y-((max(y)-min(y))/2))));
    st=[min(y)^2 sign(slope) x(Idx) max(y)^2];
end

参数C的选择意义不太直观，也不去理会

设置好拟合函数，求雅可比矩阵的函数和运行参数即可

代码

// testgsl.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_sort.h>


struct data {
	size_t n;
	double * t;
	double * y;
};

int
L4P_fit (const gsl_vector * x, void *data,
		 gsl_vector * f)
{
	size_t n = ((struct data *)data)->n;
	double *t = ((struct data *)data)->t;
	double *y = ((struct data *)data)->y;
	double A = gsl_vector_get (x, 0);
	double B = gsl_vector_get (x, 1);
	double C = gsl_vector_get (x, 2);
	double D = gsl_vector_get (x, 3);
	size_t i;
	for (i = 0; i < n; i++)
    {
		double Yi = D*D+(A*A-D*D)/(1+pow(t[i]/(C),B));
		gsl_vector_set (f, i, Yi - y[i]);
    }
	
	return GSL_SUCCESS;
}
int
L4P_df (const gsl_vector * x, void *data,
		gsl_matrix * J)
{
	size_t n = ((struct data *)data)->n;
	double *t = ((struct data *)data)->t;
	double A = gsl_vector_get (x, 0);
	double B = gsl_vector_get (x, 1);
	double C = gsl_vector_get (x, 2);
	double D = gsl_vector_get (x, 3);
	size_t i;
	
	for (i = 0; i < n; i++)
    {
		/* Jacobian matrix J(i,j) = dfi / dxj, */
		double x = t[i];
		//py_A = (2*A)/((x/C)^B + 1)
		//py_B = -(log(x/C)*(A^2 - D^2)*(x/C)^B)/((x/C)^B + 1)^2
		//py_C = (B*x*(A^2 - D^2)*(x/C)^(B - 1))/(C^2*((x/C)^B + 1)^2)
		//py_D = 2*D - (2*D)/((x/C)^B + 1)
		double DD = pow(D,2);
		double AA = pow(A,2);
		double CC = pow(C,2);
		double pow_xc_b = pow((x/C),B);
		double py_A =  (2*A)/(pow_xc_b + 1); 
		double py_B = -log(x/C)*pow_xc_b*(AA-DD)/pow((pow_xc_b + 1),2); 
		double py_C = (B*x*pow((x/C),(B - 1))*( AA -DD))/(CC*pow((pow_xc_b + 1),2)); 
		double py_D = 2*D - 2*D/(pow_xc_b + 1);
		gsl_matrix_set (J, i, 0, py_A);
		gsl_matrix_set (J, i, 1, py_B);
		gsl_matrix_set (J, i, 2, py_C);
		gsl_matrix_set (J, i, 3, py_D);
    }
	
	return GSL_SUCCESS;
}

void
callback(const size_t iter, void *params,
         const gsl_multifit_nlinear_workspace *w)
{
	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);
	double a,b,c,d,e,fx;
	a = gsl_vector_get(x, 0);
	b = gsl_vector_get(x, 1);
	c = gsl_vector_get(x, 2);
	d = gsl_vector_get(x, 3);
	e = 1.0 / rcond;
	fx = gsl_blas_dnrm2(f);  
	fprintf(stderr, "iter %d: A = %.4f, B = %.4f, C = %.4f, D = %.4f,cond(J) = %8.4f, |f(x)| = %.4f\n",
		iter,
		a*a,
		b,
		c,
		d*d,
		1.0 / rcond,
		gsl_blas_dnrm2(f));
}


/**/
static void initL4P(double *x,double *y,int n,double *a,double *b,double *c,double *d)
{
	int i = 0,idx = 0;
	   double max_y=-1,min_y=0x7fff,mind = 0x7fff,distance=0;
	   for(i=0;i<n;i++){
		   max_y = max_y>y[i]?max_y:y[i];
		   min_y = min_y<y[i]?min_y:y[i];
	   }
	   *a = sqrt(min_y);
	   *b = ((y[0]-y[n-1])/(x[0]-x[n-1])<0)?-1:1;
       *d = sqrt(max_y);
       distance = (max_y-min_y)/2;
       mind = max_y;
	   for(i=0;i<n;i++){
		   if (abs(y[i]-distance)<mind){
			   mind = abs(y[i]-distance);
			   idx = i;
		   }
	   }
	   *c = x[idx];
}
void  L4P(double *xDat,double *yDat,int n)
{
	
	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();
	const size_t p = 4;
	gsl_vector *f;
	gsl_matrix *J;
	gsl_matrix *covar = gsl_matrix_alloc (p, p);
	gsl_matrix_set_zero(covar);
	struct data d = { n, xDat, yDat };
	double x_init[4] = { 1, 1, 1,1 }; /* starting values */
	gsl_vector_view x;
	double chisq, chisq0;
	int status, info;
	size_t i;
	
	const double xtol = 1e-8;
	const double gtol = 1e-8;
	const double ftol = 0.0;
	
	initL4P(xDat,yDat,n,&x_init[0],&x_init[1],&x_init[2],&x_init[3]);
	x =  gsl_vector_view_array (x_init, p);
	/* define the function to be minimized */
	fdf.f = L4P_fit;
	fdf.df = L4P_df;   /* set to NULL for finite-difference Jacobian */
	fdf.fvv = NULL;     /* not using geodesic acceleration */
	fdf.n = n;
	fdf.p = p;
	fdf.params = &d;
	
	/* 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_winit (&x.vector, NULL, &fdf, w);
	
	/* compute initial cost function */
	f = gsl_multifit_nlinear_residual(w);
	gsl_blas_ddot(f, f, &chisq0);
	
	/* solve the system with a maximum of 100 iterations */
	status = gsl_multifit_nlinear_driver(50000, xtol, gtol, ftol,
		callback, NULL, &info, w);
	
	/* compute covariance of best fit parameters */
	J = gsl_multifit_nlinear_jac(w);
	gsl_multifit_nlinear_covar (J, 0.0, covar);
	
	/* compute final cost */
	gsl_blas_ddot(f, f, &chisq);
	
#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
	
	fprintf(stderr, "summary from method '%s/%s'\n",
		gsl_multifit_nlinear_name(w),
		gsl_multifit_nlinear_trs_name(w));
	fprintf(stderr, "number of iterations: %d\n",
		gsl_multifit_nlinear_niter(w));
	fprintf(stderr, "function evaluations: %d\n", fdf.nevalf);
	fprintf(stderr, "Jacobian evaluations: %d\n", fdf.nevaldf);
	fprintf(stderr, "reason for stopping: %s\n",
		(info == 1) ? "small step size" : "small gradient");
	fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
	fprintf(stderr, "final   |f(x)| = %f\n", sqrt(chisq));
	
	{
		double dof = n - p;
		double c = GSL_MAX_DBL(1, sqrt(chisq / dof));		
		fprintf(stderr, "chisq/dof = %f\n", chisq / dof);	
		fprintf (stderr, "A      = %.5f +/- %.5f\n", FIT(0)*FIT(0), c*ERR(0));
		fprintf (stderr, "B      = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
		fprintf (stderr, "C      = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
		fprintf (stderr, "D      = %.5f +/- %.5f\n", FIT(3)*FIT(3), c*ERR(3));
	}
	
	fprintf (stderr, "status = %s\n", gsl_strerror (status));	
	gsl_multifit_nlinear_free (w);
	gsl_matrix_free (covar);
}
int main(int argc, char* argv[])
{
	double yData[] = {3.90,5.30,7.20,9.60,12.9,17.1,23.1,31.4,38.6,50.2,62.9,
		76,92,105.7,122.8,131.7,150.7,179,205,226.5,248.7};
	double xData[]={1790,1800,1810,1820,1830,1840,1850,1860,1870,1880,
		1890,1900,1910,1920,1930,1940,1950,1960,1970,1980,1990};
	int len  = sizeof(xData)/sizeof(double);
	L4P( xData, yData,len);
	return 0;
}

结果

结果如下，和matlab的L4P.m代码结果基本一致

GSL对于对那些关注算法实现和细节的朋友这是一个福利，就是在windows下编译不太容易，网上有很多借助Mingw等工具编译的方法，好像也不太方便。

光是为了跑这个LMF算法就手动从GSL加了这么多的文件，当然这是一个笨方法，也是一个速成的办法。