GSL中的线性最小二乘拟合

线性最小二乘拟合

本章描述使用线性组合函数对实验数据进行最小二乘拟合的例程。数据可以是加权的，也可以是未加权的，即带有已知的或未知的错误。对于加权数据，函数计算最佳拟合参数及其相关的协方差矩阵。对于未加权数据，由点的离散度估计协方差矩阵，给出一个方差-协方差矩阵。

这些函数分为简单单参数或双参数回归和多参数拟合的版本。

40.1 概述

最小二乘拟合是，通过χ2 (卡方分布)的最小化得出的，卡方分布是模型Y (c, x)的n个实验数据(xi, yi)的残差平方的加权和。

 模型的p参数为c = {c0, c1，…}。加权因子wi由wi = 1/σi2给出，其中σi为数据点yi上的实验误差。假设误差为高斯分布且不相关。。对于未加权数据，卡方和的计算没有任何权重因子。

拟合程序返回最佳拟合参数c及其p×p协方差矩阵。协方差矩阵度量由数据上的误差引起的最佳拟合参数上的统计误差σi，定义为

 式中 表示基础数据点的高斯误差分布的平均值。

     协方差矩阵由数据误差σi通过误差传播计算得到。给出了由数据δyi 的微小变化引起的拟合参数δca 的变化，

 允许根据数据上的误差来写协方差矩阵，

 对于不相关的数据，基础数据点的波动满足

 给出相应的参数协方差矩阵

 当计算协方差矩阵为无权重数据时，即数据与未知错误，和里面的权重因子wi 使用的是单一估计w = 1 /σ2，其中σ2是最佳模型计算残差的方差， 。这被称为方差-协方差矩阵。

最佳拟合参数的标准差由协方差矩阵对应对角元素σca =Caa 的平方根给出。拟合参数ca和cb的相关系数为

40.2 线性回归

本节中的函数用于拟合简单的一或两个参数线性回归模型。这些函数声明在头文件gsl_fit.h中。

40.2.1 带常数项的线性回归

本节描述的函数可用于对直线模型Y (c, x) = c0 + c1x执行最小二乘拟合。

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, y)计算模型Y = c0 + c1X的最佳拟合线性回归系数(c0, c1)，(x, y)是步长xstride和ystride的两个长度为n的向量。假设y上的误差是未知的，因此参数(c0, c1)的方差-协方差矩阵是从最佳拟合线周围的点的散点估计出来的，并通过参数(cov00, cov01, cov11)返回。最佳拟合直线的残差平方和在sumsq中返回。注意：数据的相关系数可以使用gsl_stats_correlation()计算，它不依赖于拟合。

int gsl_fit_wlinear(const double * x, const size_t xstride, const double * w,

const size_t wstride, const double * y, const size_t ystride,

size_t n, double * c0, double * c1, double * cov00,

double * cov01, double * cov11, double * chisq)

本函数用加权数据集(x, y)计算模型Y = c0 + c1X的最佳拟合线性回归系数(c0, c1)，(x, y)是步长xstride和ystride的两个长度为n的向量。向量w长度是n，步长wstride，指定了每个数据点的权重。权重是y中每个数据点方差的倒数。

参数(c0, c1)的协方差矩阵使用权重计算，并通过参数(cov00, cov01, cov11)返回。最佳拟合直线的残差的加权平方和χ2 以chisq返回。

int gsl_fit_linear_est(double x, double c0, double c1, double cov00, double cov01,

double cov11, double * y, double * y_err)

    本函数使用最佳拟合线性回归系数c0、c1及其协方差cov00、cov01、cov11计算模型Y = c0 + c1X在x点处的拟合函数y及其标准差y_err。

40.2.2 不带常数项的线性回归

    本节描述的函数可用于对没有常数项Y = c1X的直线模型执行最小二乘拟合。

int gsl_fit_mul(const double * x, const size_t xstride, const double * y,

const size_t ystride, size_t n, double * c1, double * cov11,

double * sumsq)

    本函数计算数据集(x, y)的模型Y = c1X的最佳拟合线性回归系数c1，(x, y)是步长xstride和ystride的两个长度为n的向量。假设y上的误差是未知的，因此参数c1的方差是由最佳拟合线周围的点的散点估计的，并通过参数cov11返回。最佳拟合直线的残差平方和在sumsq中返回。

int gsl_fit_wmul(const double * x, const size_t xstride, const double * w,

const size_t wstride, const double * y, const size_t ystride,

size_t n, double * c1, double * cov11, double * sumsq)

本函数计算加权数据集(x, y)模型Y = c1X的最佳拟合线性回归系数c1，(x, y)是步长xstride和ystride的两个长度为n的向量。向量w长度是n，步幅wstride，指定每个数据点的权重。权重是y中每个数据点方差的倒数。

使用权重计算参数c1的方差，并通过参数cov11返回。最佳拟合直线的残差的加权平方和χ2 以chisq返回。

int gsl_fit_mul_est(double x, double c1, double cov11, double * y, double * y_err)

    本函数使用最佳拟合线性回归系数c1及其协方差cov11计算模型Y = c1X在x点处的拟合函数y及其标准差y_err。

40.3 多参数回归

本节描述通过最小化成本函数来执行最小二乘拟合线性模型的例程，

     其中，y是n个观测值的向量，X是n乘p的预测变量矩阵，c是p个待估计的未知最佳拟合参数的向量，。矩阵W=diag(w1,w2,...,wn)定义了观测向量的权重或不确定性。

     通过适当地准备n乘p矩阵X，这个公式可以用于适合任何数量的函数和(或)变量。例如，要拟合x的一个p阶多项式，使用下面的矩阵，

     索引i遍历观察值，索引j从0到p−1。

拟合一组固定频率的p正弦函数ω1, ω2, . . ., ωp，使用，

     为了拟合p个自变量x1,x2,…,xp，使用

Xij=xj(i)

式中，xj(i) 为预测变量xj 的第i个值。

一般线性最小二乘系统的解需要为中间结果提供额外的工作空间，例如矩阵X的奇异值分解。

    这些函数声明在头文件gsl_multifit.h中。

gsl_multifit_linear_workspace

此工作空间包含用于拟合多参数模型的内部变量。

gsl_multifit_linear_workspace * gsl_multifit_linear_alloc(const size_t n,

const size_t p)

    本函数分配了一个工作空间，用于拟合一个模型，最大n个观测值，使用最大p个参数。如果需要，用户可以稍后提供一个更小的最小二乘系统。工作空间的大小为O(np + p2)。

void gsl_multifit_linear_free(gsl_multifit_linear_workspace * work)

    本函数释放与工作空间w相关的内存。

int gsl_multifit_linear_svd(const gsl_matrix * X,

gsl_multifit_linear_workspace * work)

    本函数对矩阵X进行奇异值分解，并在work内部存储SVD因子。

int gsl_multifit_linear_bsvd(const gsl_matrix * X,

gsl_multifit_linear_workspace * work)

本函数对矩阵X进行奇异值分解，并在work内部存储SVD因子。首先利用列比例因子对矩阵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)

    本函数利用work中提供的预先分配的工作空间，根据观测值y和预测变量X的矩阵，计算模型y = Xc的最佳拟合参数c。将模型参数cov的p乘p方差-协方差矩阵设为σ2(XTX)−1，其中σ为拟合残差的标准差。最佳拟合的残差平方和(χ2 )以chisq的形式返回。如果需要确定系数，可以用表达式R2 = 1−χ2 /TSS计算，其中观测值y的总平方和(TSS)可以用gsl_stats_tss()计算。

采用改进的Golub-Reinsch SVD算法对矩阵X进行奇异值分解，获得最佳拟合，并进行列缩放，提高奇异值的精度。任何奇异值为零(机器精度)的部件在拟合中被丢弃。

int gsl_multifit_linear_tsvd(const gsl_matrix * X, const gsl_vector * y,

const double tol, gsl_vector * c, gsl_matrix * cov,

double * chisq, size_t * rank,

gsl_multifit_linear_workspace * work)

    本函数对观测值y和预测变量X的矩阵使用截断SVD展开，计算模型y = Xc的最佳拟合参数c。在拟合中丢弃满足si≤tol×s0的奇异值，其中s0为最大奇异值。将模型参数cov的p乘p方差-协方差矩阵设为σ2(XTX)−1，其中σ为拟合残差的标准差。最佳拟合的残差平方和(χ2 )以chisq的形式返回。有效秩(解中使用的奇异值的数量)按秩返回。如果需要确定系数，可以用表达式R2 = 1−χ2 /TSS计算，其中观测值y的总平方和(TSS)可以用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)

本函数利用work中提供的预先分配的工作空间，计算加权模型y = Xc对观测值y的最佳拟合参数c，权重为w，预测值变量矩阵为X。模型参数cov的p乘p协方差矩阵计算为(XTWX)−1。最佳拟合残差的加权平方和(χ2 )，以chisq返回。如果想要确定系数，可以用表达式R2 = 1−χ2 /WTSS计算，其中观测值y的加权总平方和(WTSS)可以用gsl_stats_wtss()计算。

int gsl_multifit_wlinear_tsvd(const gsl_matrix * X, const gsl_vector * w,

const gsl_vector * y, const double tol, gsl_vector * c,

gsl_matrix * cov, double * chisq, size_t * rank,

gsl_multifit_linear_workspace * work)

    本函数对观测值y和预测变量X的矩阵使用截断SVD展开，计算模型y = Xc的最佳拟合参数c。在拟合中丢弃满足si≤tol×s0的奇异值，其中s0为最大奇异值。模型参数cov的p乘p协方差矩阵计算为(XTWX)−1。最佳拟合的残差的加权平方和(χ2 )以chisq的形式返回。系统的有效秩(解中使用的奇异值的数量)按秩返回。如果想要确定系数，可以用表达式R2 = 1−χ2 /WTSS计算，其中观测值y的加权总平方和(WTSS)可以用gsl_stats_wtss()计算。

int gsl_multifit_linear_est(const gsl_vector * x, const gsl_vector * c,

const gsl_matrix * cov, double * y, double * y_err)

    本函数使用最佳拟合多元线性回归系数c及其协方差矩阵cov计算模型y = x.c在x点处的拟合函数值y及其标准差y_err。

int gsl_multifit_linear_residuals(const gsl_matrix * X, const gsl_vector * y,

const gsl_vector * c, gsl_vector * r)

    本函数用观测值y、系数c和预测变量X的矩阵，计算残差向量r = y−Xc。

size_t gsl_multifit_linear_rank(const double tol,

const gsl_multifit_linear_workspace * work)

    本函数返回矩阵X的秩，它必须首先计算其奇异值分解。秩是通过计算满足σj> tol×σ0的奇异值σj的个数来计算的，其中σ0为最大奇异值。