GSL使用最小二乘

来源：互联网发布：手机pdf免费编辑软件编辑：程序博客网时间：2024/05/16 18:38

第一部分
1.这一部分介绍如何利用GSL解决线性最小二乘问题。
首先针对一个最标准的线性函数进行计算，形如 y = a * x +b;
利用到的函数是

int gsl_fit_linear (const double * x, const size_t xstride, const double * y, const size_t ystride, size_t n, double * c0, double * c1, double * cov00, double * cov01, double * cov11, double * sumsq)//x是自变量// xstrside是步长// y是因变量//ystride是y的步长//n是数据组数//c0 储存计算出来的a （f(x) = a * x + b）// c1 储存计算出来的b （f(x) = a * x + b）// cov00 cov01 cov11是点的协方差 我的博客中有介绍协方差// 没有cov10是因为它和cov01是相等的//sumsq是残差

测试程序如下

int testGslLinear(){        //sample        size_t norris_n = 11;        size_t xstride = 1;        size_t ystride = 1;        double x[] = { 71,  73,  64,  65,  61,  70,  65,  72,  63,  67,  64};        double y[11];        for(int i = 0; i <= 10; i++){            y[i] = 2 * x[i] + 10.2;        }        //result        double c0 = 0.0;        double c1 = 0.0;        double cov00 = 0.0;        double cov01 = 0.0;        double cov11 = 0.0;        double sumsq = 0.0;        //test infinite.        //for (i = 0;i < norris_n;i++) {        //    x[i] = 10;        //}        //calculate        gsl_fit_linear (x, xstride, y, ystride, norris_n,             &c0, &c1, &cov00, &cov01, &cov11, &sumsq);         {            printf("f(x) = %lf * x + %lf x\n", c0, c1);            printf("cov00 = %lf\n", cov00);            printf("cov01 = %lf\n", cov01);            printf("cov11 = %lf\n", cov11);            return 0;         }         return 0;    }

2.测试形如 f(x) = a * g(x) + b的方程
例如 f(x)=a∗ex+b，直观看他并不是关于x的线性方程，但是如果把ex看做g(x)，那么就是关于g(x)的线性方程，这种观点可以大大扩宽最小二乘的应用范围。
首先介绍要利用的函数

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)//x是自变量 注意这儿是矩阵形式 每一行代表一个自变量的输入向量//y是因变量  每一行代表相应x输入的结果//c是结果 大小等于x的列向量的长度//cov是系数的variance-covariance矩阵 是M*M (M是x的列长）//chisq是残差//work是工作空间 长宽与x的长宽相同

下面是测试代码

int testGSLMultilinear(){        // y = 2 * √x1 + 6 * x0        double xvalue[] = { 33, 66, 145, 68, 23};        double yvalue[5];        for(int i=0; i < 5; i++){            yvalue[i] = sqrt(xvalue[i]) * 2 + 6;        }        //申请矩阵空间        //y = a * √x1 + b * x0        gsl_matrix* A = gsl_matrix_alloc(5, 2);        gsl_vector* b = gsl_vector_alloc(5);        for(int i=0; i<5; i++){            gsl_matrix_set(A, i, 0, 1);            gsl_matrix_set(A, i, 1, pow(xvalue[i], 0.5) );            gsl_vector_set(b, i, yvalue[i]);        }        gsl_vector* c = gsl_vector_alloc(2);        gsl_matrix* cov = gsl_matrix_alloc(2, 2);        double residual;        //申请工作空间        gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(5, 2);        gsl_multifit_linear(A, b, c, cov, &residual, work);        printf("y = %lf + %lf * pow(x, 0.5) \n", c->data[0], c->data[1]);        gsl_matrix_free(A);        gsl_vector_free(b);        gsl_matrix_free(cov);        gsl_vector_free(c);        gsl_multifit_linear_free (work);        return 0;    }

另外，两个都需要的头文件是

#include <gsl/gsl_fit.h>#include <gsl/gsl_linalg.h>#include <gsl/gsl_multifit.h>

1 0