一种极小卡方估计方法与流程

本发明涉及一种极小卡方估计方法，属于数值分析和可靠性技术领域。

背景技术：

在科学试验和工程应用中有很多函数，它的解析表达式是不知道的。比如，惯性器件随着贮存时间的延长，其性能和参数会有所变化，但该变化过程仅能通过实验观测的方法测得一系列节点xi上的值fi。现在的问题是寻求测试点序列的逼近函数y(x)，或用几何语言来说就是寻求一条曲线y(x)来拟合(平滑)这n个点。

在拟合度检验的相关理论中，χ²统计量经常用于检验一组独立样本的共同分布是否属于某一具有特定性质的分布族。极小卡方估计是指将极小化χ²统计量所得到的参数作为真值的最佳估计，写成表达式为

其中，

则称y(x)为f(x)关于权系数{ωi}的卡方逼近，并称上述准则求逼近函数y(x)的方法为卡方估计法。在上式中，{ωj(x)}为线性无关的函数集，{aj}为对应的系数集。

定义m+1元函数h(a0,a1,…,am)为

由多元函数极值的必要条件可知，最小值点应满足

上式等效为

但由于上式为一个多元非线性方程组，直接求解系数a0,a1,…,am很困难。

在参考文献“温度应力下的加速度计贮存寿命评定”(《装甲兵工程学院学报》，vol.28no.3,2014，p28)中提出一种方法，通过matlab软件中的fminsearch函数求解最小值。在文中介绍到，“在matlab中，函数fminsearch用来求解无约束多维极值问题，其采用单纯形搜索法，不需要计算梯度”。

通过在实际应用中对该函数的应用，发现其存在如下缺点：

(1)解不唯一，这就意味着该方法并不能确定使卡方最小的系数值；

(2)方法机理不清楚，算法不是最优。

(3)容易造成信息丢失。

因此，需要针对上述情况，提出一种改进的明确的极小卡方估计方法。

技术实现要素：

本发明的目的在于克服现有技术的上述缺陷，提供一种极小卡方估计方法，该方法能够通过精确迭代解算的方式，给出拟合函数的各项系数，满足卡方极小化的要求。

本发明的上述目的主要是通过如下技术方案予以实现的：

一种极小卡方估计方法，包括：

给定测量序列{fi}和权系数序列{ωi}，i＝1,2,…,n；

确定线性无关的基函数集

构造测量序列{fi}的曲线拟合函数ψ(x)，令其中{aj}为对应的未知系数集；

根据最优准则确定所述系数集{aj}的最优值{a*j}，即根据如下公式求解最优值{a*j}：

将所述最优值{a*j}代入曲线拟合函数ψ(x)，得到极小卡方估计的最优拟合曲线的函数表达式y(x)：

在上述估计方法中，求解最优值{a*j}的具体方法如下：

(1)、根据如下公式求解系数a0,a1,…,am的初值；

(2)、把系数a0,a1,…,am的初值代入如下公式求解gi：

(3)、将gi代入如下公式(5)，并根据如下公式(5)求解得到系数a0,a1,…,am的一组新值：

(4)、返回步骤(2)，将所述系数a0,a1,…,am的一组新值作为初值代入公式(4)，循环计算，直至计算得到的系数a0,a1,…,am的第l组新值与上一次计算得到的系数a0,a1,…,am的第l-1组新值之间的偏差满足要求，则取所述第l组新值作为最优值{a*j}；其中：l为正整数。

在上述估计方法中，所述步骤(4)中系数a0,a1,…,am的第l组新值与上一次计算得到的系数a0,a1,…,am的第l-1组新值之间的偏差满足要求是指：系数a0,a1,…,am的第l组新值与系数a0,a1,…,am的第l-1组新值中相对应的每一个系数在小数点后前p位的值均相同，即第l组新值中的a0与第l-1组新值中的a0在小数点后前p位的值相同，第l组新值中的a1与第l-1组新值中的a1在小数点后前p位的值相同，……，第l组新值中的am与第l-1组新值中的am在小数点后前p位的值相同，其中：p为正整数。

在上述估计方法中，述步骤(1)中系数a0,a1,…,am的初值通过如下公式求解得到：

在上述估计方法中，所述步骤(3)中系数a0,a1,…,am的一组新值通过如下公式求解得到：

本发明与现有技术相比的有益效果如下：

(1)、本发明提出的极小卡方估计方法，在线性无关的基函数集的基础上构造拟合曲线的函数形式，采用迭代计算的方法给出拟合曲线的函数中各项系数，把系数代入拟合函数后得到的理论值具有卡方最小的特点。

(2)、本发明提出的极小卡方估计方法，给出了一种不同于运用matlab软件中的fminsearch函数的方法，避免了fminsearch函数存在多解问题，具有唯一解且能保证卡方最小的特点。

(3)、本发明采用迭代计算求解拟合函数的系数，解决了非线性方程组无法求解的问题；并且保证了求解的准确性，避免了求解过程中存在的信息缺失问题。

(4)、本发明给出的极小卡方估计曲线拟合方法，逻辑性强，易于编程，且方法简单，易于实现。

附图说明

图1为本发明基于相对误差极小的曲线拟合流程图；

图2为本发明实施例中给定的测量序列；

图3为本发明实施例的迭代过程；

图4为本发明实施例得到的拟合曲线；

图5为本发明实施例得到的拟合曲线与最小二乘法比较图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细的描述：

本发明提供了一种极小卡方估计方法，在线性无关的基函数集的基础上构造拟合曲线的函数形式，采用迭代计算的方法给出拟合曲线的函数中各项系数，把系数代入拟合函数后得到的理论值具有卡方最小的特点。

如图1所示为本发明基于相对误差极小的曲线拟合流程图，本发明极小卡方估计方法包括如下步骤，测量序列与拟合序列的相对误差平方和最小：

一、给定测量序列{fi}和权系数序列{ωi}，i＝1,2,…,n；

二、确定线性无关的基函数集

三、构造测量序列{fi}的曲线拟合函数ψ(x)，令其中{aj}为对应的未知系数集；

四、根据最优准则确定所述系数集{aj}的最优值{a*j}。为使卡方最小，有

定义m+1元函数h(a0,a1,…,am)为

由多元函数极值的必要条件可知，最小值点应满足

上式等效为

但由于在上式中，各项系数a0,a1,…,am出现在分母中，难以求解。拟采用迭代求解。可以通过如下方法求解最优值{a*j}，具体包括如下步骤：

(1)、对下式求解系数a0,a1,…,am的初值；

定义m+1元函数h1(a0,a1,…,am)为

由多元函数极值的必要条件可知，最小值点应满足

上式等效为

可得对系数初值a0,a1,…,am的一种计算方法为

(2)、把系数a0,a1,…,am的初值代入如下公式求解gi：

(3)、将gi代入如下公式(5)，并根据如下公式(5)求解得到系数a0,a1,…,am的一组新值：

公式(5)推导过程如下：

上式等效为：

即：

可得对系数a0,a1,…,am的一种计算方法为：

(4)、返回步骤(2)，将该系数a0,a1,…,am的一组新值作为初值代入公式(4)，再次计算gi，将gi代入公式(5)中求解得到系数a0,a1,…,am的第二组新值，再返回步骤(2)，将该第二组新值作为初值代入公式(4)，依次类推，循环计算，直至计算得到的系数a0,a1,…,am的第l组新值与上一次计算得到的系数a0,a1,…,am的第l-1组新值之间的偏差满足要求，即各项误差系数趋于一个固定值，则取所述第l组新值作为最优值{a*j}；其中：l为正整数。

上述步骤(4)中系数a0,a1,…,am的第l组新值与上一次计算得到的系数a0,a1,…,am的第l-1组新值之间的偏差满足要求是指：系数a0,a1,…,am的第l组新值与系数a0,a1,…,am的第l-1组新值中相对应的每一个系数在小数点后前p位的值均相同，即第l组新值中的a0与第l-1组新值中的a0在小数点后前p位的值相同，第l组新值中的a1与第l-1组新值中的a1在小数点后前p位的值相同，……，第l组新值中的am与第l-1组新值中的am在小数点后前p位的值相同，其中：p为正整数，例如本实施例中p取10，即m+1个系数小数点后前10均对应相同，则认为满足偏差要求。

五、将所述最优值{a*j}代入曲线拟合函数ψ(x)，得到极小卡方估计的最优拟合曲线的函数表达式y(x)：

采用本发明方法在实际应用中可以对测试的数据进行曲线拟合，用于对设备进行寿命评估，判断和预测设备的有效期和贮存寿命。

实施例1

设给定自变量为{xi}对应的测量序列{fi}，i＝1,2,…,n；其中n＝10，如图2所示为本发明实施例中给定的测量序列；数据具体值见表1。采用线性曲线进行拟合，选取两个基函数和由此构造的拟合函数为y(x)＝a1x-a2。

表1

利用本发明的求解方法，初值为a1＝1.857130151514106、a2＝6.939083333409633，再经过8次迭代计算后最终求得a1＝1.854996089750836、a2＝6.975801527797882。即与第7次迭代计算得到的a1、a2小数点后前10位对应相等，为最优值。迭代过程如图3所示，其中图3a与图3b分别表示a1、a2的迭代过程。把a1和a2代入拟合函数为y(x)＝a1x-a2中求得的理想数据如图4中“*”点所示，图4所示为本发明实施例得到的拟合曲线，图中“。”表示为测量值。

对表1的数据也可采用最小二乘法，求解的系数为a1＝1.696152760852661、a2＝6.574443934342034。如图5所示为本发明实施例得到的拟合曲线与最小二乘法比较图，图5中“虚线”为最小二乘法得到的拟合曲线，“实线”为本发明方法得到的拟合曲线，“点划线”为测量序列{fi}的拟合曲线。为了比较两种方法的区别，分别求取

(1)最小二乘法：-0.4517058218079164；

(2)本发明方法：-0.5265079715177377。

可以看出，本发明给出的卡方最小，具有唯一解，且无信息缺失问题。

以上所述，仅为本发明一个具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明未详细说明部分属于本领域技术人员公知常识。