GSL数学库解多参数方程

在之前的博文里介绍了gsl库的安装使用，本篇介绍使用gsl库进行多参数方程的线性拟合求解方法，使用手册给出的实例，主要参考和学习调用gsl多参线性拟合的使用步骤：

使用下面的程序来拟合二次方程：，矩阵形式为：y=X.*c主要语句为：gsl_multifit_wlinear()，其中，模型中的X为：

c=[c0,c1,c2]，程序实例：

#include <stdio.h>
#include <gsl/gsl_multifit.h>

int main (int argc, char **argv)
{
  int i, n;
  double xi, yi, ei, chisq;
  gsl_matrix *X, *cov;
  gsl_vector *y, *w, *c;

  if (argc != 2)
    {
      fprintf (stderr,"usage: fit n < data\n");
      exit (-1);
    }

  n = atoi (argv[1]);
//初始化X,y,w,c,cov等参数，申请空间
  X = gsl_matrix_alloc (n, 3);
  y = gsl_vector_alloc (n);
  w = gsl_vector_alloc (n);

  c = gsl_vector_alloc (3);
  cov = gsl_matrix_alloc (3, 3);

  for (i = 0; i < n; i++)//读入n行数据（x,y,ei）,ei为y的方差
    {
      int count = fscanf (stdin, "%lg %lg %lg",
                          &xi, &yi, &ei);

      if (count != 3)
        {
          fprintf (stderr, "error reading file\n");
          exit (-1);
        }

      printf ("%g %g +/- %g\n", xi, yi, ei);

      gsl_matrix_set (X, i, 0, 1.0);//将数值1赋给X（i,0）
      gsl_matrix_set (X, i, 1, xi); //将数值xi赋给X（i,1）
      gsl_matrix_set (X, i, 2, xi*xi); //将数值xi*xi赋给X(i,2)

      gsl_vector_set (y, i, yi);//类似上面
      gsl_vector_set (w, i, 1.0/(ei*ei));
    }

  {
    gsl_multifit_linear_workspace * work
      = gsl_multifit_linear_alloc (n, 3);//使用拟合函数的第一步，申请运行函数的空间
    gsl_multifit_wlinear (X, w, y, c, cov,
                          &chisq, work);
    gsl_multifit_linear_free (work); //释放空间
  }

#define C(i) (gsl_vector_get(c,(i)))
#define COV(i,j) (gsl_matrix_get(cov,(i),(j)))

  {
    printf ("# best fit: Y = %g + %g X + %g X^2\n",
            C(0), C(1), C(2));

    printf ("# covariance matrix:\n");
    printf ("[ %+.5e, %+.5e, %+.5e  \n",
               COV(0,0), COV(0,1), COV(0,2));
    printf ("  %+.5e, %+.5e, %+.5e  \n",
               COV(1,0), COV(1,1), COV(1,2));
    printf ("  %+.5e, %+.5e, %+.5e ]\n",
               COV(2,0), COV(2,1), COV(2,2));
    printf ("# chisq = %g\n", chisq);
  }

  gsl_matrix_free (X);
  gsl_vector_free (y);
  gsl_vector_free (w);
  gsl_vector_free (c);
  gsl_matrix_free (cov);

  return 0;
}

相关语句命令：

gsl_matrix * gsl_matrix_alloc(size_t n1, size_t n2)

This function creates a matrix of size n1 rows by n2 columns, returning a pointer to a newly initialized matrix struct. A new block is allocated for the elements of the matrix, and stored in the block component of the matrix struct. The block is “owned” by the matrix, and will be deallocated when the matrix is deallocated. Requesting zero for n1 or n2 is valid and returns a non-null result.

gsl_vector * gsl_vector_alloc(size_t n)

This function creates a vector of length n, returning a pointer to a newly initialized vector struct. A new block is allocated for the elements of the vector, and stored in the block component of the vector struct. The block is “owned” by the vector, and will be deallocated when the vector is deallocated. Zero-sized requests are valid and return a non-null result.

void gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, double x)

This function sets the value of the (i,j)-th element of a matrix m to x. If i or j lies outside the allowed range of 0 to n1 - 1 and 0 to n2 - 1 then the error handler is invoked. An inline version of this function is used when HAVE_INLINE is defined.

void gsl_vector_set(gsl_vector * v, const size_t i, double x)

This function sets the value of the i-th element of a vector v to x. If i lies outside the allowed range of 0 to size - 1 then the error handler is invoked. An inline version of this function is used when HAVE_INLINE is defined.

gsl_multifit_linear_workspace * gsl_multifit_linear_alloc(const size_t n, const size_t p)

This function allocates a workspace for fitting a model to a maximum of n observations using a maximum of p parameters. The user may later supply a smaller least squares system if desired. The size of the workspace is O(np + p^2)

int gsl_multifit_linear(const gsl_matrix * X, const gsl_vector * y, gsl_vector * c, gsl_matrix * cov, double * chisq, gsl_multifit_linear_workspace * work)

This function computes the best-fit parameters c of the model y = X c for the observations y and the matrix of predictor variables X, using the preallocated workspace provided in work. The p-by-p variance-covariance matrix of the model parameters cov is set to , where σ is the standard deviation of the fit residuals. The sum of squares of the residuals from the best-fit, , is returned in chisq. If the coefficient of determination is desired, it can be computed from the expression, where the total sum of squares (TSS) of the observations y may be computed from gsl_stats_tss().

int gsl_multifit_wlinear(const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, gsl_vector * c, gsl_matrix * cov, double * chisq, gsl_multifit_linear_workspace * work)

This function computes the best-fit parameters c of the weighted model y = X c for the observations y with weights w and the matrix of predictor variables X, using the preallocated workspace provided in work. The p-by-p covariance matrix of the model parameters cov is computed as . The weighted sum of squares of the residuals from the best-fit, , is returned in chisq. If the coefficient of determination is desired, it can be computed from the expression , where the weighted total sum of squares (WTSS) of the observations y may be computed from gsl_stats_wtss().

double gsl_vector_get(const gsl_vector * v, const size_t i)

This function returns the i-th element of a vector v. If i lies outside the allowed range of 0 to size - 1 then the error handler is invoked and 0 is returned. An inline version of this function is used whenHAVE_INLINE is defined.

gsl_vector_view gsl_matrix_row(gsl_matrix * m, size_t i)
gsl_vector_const_view gsl_matrix_const_row(const gsl_matrix * m, size_t i)

These functions return a vector view of the i-th row of the matrix m. The data pointer of the new vector is set to null if i is out of range. The function gsl_matrix_const_row() is equivalent to gsl_matrix_row() but can be used for matrices which are declared const.

int  gsl_multifit_linear_est(const gsl_vector * x, const gsl_vector * c, const gsl_matrix * cov, double * y, double * y_err)

This function uses the best-fit multilinear regression coefficients c and their covariance matrix cov to compute the fitted function value y and its standard deviation y_err for the model y = x.c at the point x.

注释：

程序中最核心的语句是：gsl_multifit_wlinear (X, w, y, c, cov,&chisq, work) ; 其中，输出参数c即为所要求解的参数。其他语句主要负责构建方程：y=Xc。当需要求解实际问题时，即可按照实例中的步骤进行求解。