最小二乘法
设(x 1, y 1 ), (x 2, y 2), …, (x n, y n)是直角平面坐标系下给出
的一组数据,若x 1<x 2<…<x n,我们也可以把这组数据看作
是一个离散的函数。根据观察,如果这组数据图象“很象”一条直线
(不是直线),我们的问题是确定一条直线y = bx +a ,使得它能" 最好" 的反映出这组数据的变化。
最小二乘法是处理各种观测数据进行测量平差的一种基本方法。
如果以不同精度多次观测一个或多个未知量,为了求定各未知量的最可靠值,各观测量必须加改正数,使其各改正数的平方乘以观测值的权数的总和为最小。因此称最小二乘法。所谓“权”就是表示观测结果质量相对可靠程度的一种权衡值。
法国数学家勒让德于1806年首次发表最小二乘理论。事实上,德国的高斯于1794年已经应用这一理论推算了谷神星的轨道,但迟至1809年才正式发表。此后他又提出平差三角网的理论,拟定了解法方程式的方法等。为利用最小二乘法测量平差奠定了基础。
最小二乘法也是数理统计中一种常用的方法,在工业技术和其他科学研究中有广泛应用。
在我们研究两个变量(x, y)之间的相互关系时,通常可以得到一系列成对的数据(x1, y1、x2, y2... xm , ym) ;将这些数据描绘在x -y直角坐标系中(如图1), 若发现这些点在一条直线附近,可以令这条直线方程如(式1-1) 。
Y 计= a0 + a1 X (式1-1)
其中:a0、a1 是任意实数
为建立这直线方程就要确定a0和a1,应用《最小二乘法原理》,将实测值Yi 与利用(式1-1) 计算值(Y计= a0 + a1 X) 的离差(Yi - Y 计) 的平方和`〔∑(Yi - Y计)2〕最小为“优化判据”。
令: φ = ∑(Yi - Y计)2 (式1-2)
把(式1-1) 代入(式1-2) 中得:
φ = ∑(Yi - a0 - a1 Xi)2 (式1-3)
当∑(Yi-Y 计) 平方最小时,可用函数 φ 对a0、a1求偏导数,令这两个偏导数等于零。
(式1-4)
(式1-5) (见附图)
亦即:
m a0 + (∑Xi ) a1 = ∑Yi (式1-6)
(∑Xi ) a0 + (∑Xi2 ) a1 = ∑(Xi, Yi) (式1-7)
得到的两个关于a0、 a1为未知数的两个方程组,解这两个方程组得出:
a0 = (∑Yi) / m - a1(∑Xi) / m (式1-8)
a1 = [∑Xi Yi - (∑Xi ∑Yi)/ m] / [∑Xi2 - (∑Xi)2 / m)] (式1-9)
这时把a0、a1代入(式1-1) 中, 此时的(式1-1) 就是我们回归的元线性方程即:数学模型。
在回归过程中,回归的关联式是不可能全部通过每个回归数据点(x1, y1、 x2, y2...xm,ym),为了判断关联式的好坏, 可借助相关系数“R”,统计量“F”,剩余标准偏差“S”进行判断;“R”越趋近于 1 越好;“F”的绝对值越大越好;“S”越趋近于 0 越好。
R = [∑XiYi - m (∑Xi / m)(∑Yi / m)]/ SQR{[∑Xi2 - m (∑Xi / m)2][∑Yi2 - m (∑Yi / m)2]} (式1-10) *
在(式1-1) 中,m 为样本容量,即实验次数;Xi 、Yi 分别任意一组实验X 、Y 的数值。
如果您认为本词条还有待完善,需要补充新内容或修改错误内容,请
编辑词条
开放分类: 数学、物理
参考资料:
1.(2)/shuzhifenxi/ch11ff5.htm
2.
贡献者: 、、、、
本词条在以下词条中被提及: 高斯、回归分析、正交
关于本词条的评论(共3条):
·这个解释的数字式写的不是很好。我再看的时候根本就有一点看不懂。把东西改明朗一下应该会好一点。 飞尔一万年 10-18 00:03
·好复杂 如果我要用最小二乘法来算直流电机的调速时候的那个速度怎么用啊,使他实际的与给定的误差很小呢!是不是之前就要吧那个直线的方程参数选好啊 790732342 08-22 17:29
·原来是这样~ 12-14 12:53
最小二乘法直线拟合
减小字体 增大字体 作者:佚名 来源:本站整理 发布时间:2006-3-27 16:47:03
//最小二乘法直线拟合
BOOL CalculateLineKB(CFoldPointList *m_FoldList,double &k,double &b)
{
//最小二乘法直线拟合
//m_FoldList为关键点(x,y)的链表
//拟合直线方程(Y=kX+b)
if(m_FoldList==NULL)return FALSE;
long lCount=m_FoldList->GetCount();
if(lCount
CFoldPoint *pFold;
double mX,mY,mXX,mXY,n;
mX=mY=mXX=mXY=0;
n=lCount;
POSITION pos=m_FoldList->GetHeadPosition();
while(pos != NULL)
{
pFold=m_FoldList->GetNext(pos);
mX+=pFold->X;
mY+=pFold->Y;
mXX+=pFold->X*pFold->X;
mXY+=pFold->X*pFold->Y;
}
if(mX*mX-mXX*n==0)return FALSE;
k=(mY*mX-mXY*n)/(mX*mX-mXX*n);
b=(mY-mX*k)/n;
return TRUE;
}
用最小二乘法拟合直线的问题 回答清楚追加分!
悬赏分:0 - 解决时间:2006-12-27 13:40
我用最小二乘法拟合直线 y=ax+b ,但是在计算a 的时候,我用偏差为最小,然后求偏导,得出来的公式进行计算。可是运算结果没法得到垂直于x 轴的直线。(当所有的点都在y 轴的时候可以得到a 的分母为0,可以判断)但是如果点分布在y 轴两侧就不行了。会得到a = 0的结果。但是这个结果肯定定不对。 怎么办啊?谢谢大家了!
问题补充:请看清我的问题好吗?我是说用最小二乘法拟合直线时垂直于x 轴的直线怎么处理?请不要随便帖无关的代码。请提供判断方法。谢谢!
提问者: summonerx - 见习魔法师 二级 最佳答案
原理中有一类题目,对测量数据进行处理,然后使用最小二乘法对数据进行处理并且拟合一条曲线,以方便对数据结果进行进一步的处理。这个程序拟合的是直线,用于处理近似线性的数据。下面是源程序,至少可以运行,会不会有问题就不知道了噻。程序是用C 语言写的,但是注释的风格是C++的,在某些编译器下,如TC 可能会有问题,把 // 换成 /* */就可以了。
#include
#include
#define N 20 //定义最多能够处理的数据组数
//变量X ,Y 线性方程系数k 线性方程矩阵m0 m1 m2
double x[N],y[N],k[2][3],m1,m2,m0;
int i=0,j=0;
//求(A1*B1~Ac*Bc)的和
double fsum(double a[],double b[],int c)
{
double sum=0;
for(i=0;i
sum+=a[i]*b[i];
return sum;
}
//求矩阵
double fmatrix(int m,int n)
{
double matrix;
matrix=k[0][m]*k[1][n]-k[0][n]*k[1][m];
return matrix;
}
void main()
{
int limit=0; //数据组数
double mi[N]; //大小为1的数列,矩阵求和时匹配使用
char ch;
//声明
printf("This program will calculate Y=aX+b, with maximum of data group of %d",N);
for(i=0;i
mi[i]=1;
//输入数据
begin:
printf("\n\nPlease input the number of data group:");
scanf("%d",&limit);
if(limit>N)
{
printf("Out of range! Should no more than %d",N);
goto begin;
}
//按X1~Xi Y1~Yi的顺序输入
printf("Please input the data one by one:");
input:
printf("\n");
for(i=0;i
{
printf("X[%d]:",i+1);
scanf("%lf",&x[i]);
}
for(i=0;i
{
printf("Y[%d]:",i+1);
scanf("%lf",&y[i]);
}
//显示输入数据以供检查
printf("Please check the data:\n");
for(i=0;i
printf("X%2d: %-11f ",i+1,x[i]);
printf("\n");
for(i=0;i
printf("Y%2d: %-11f ",i+1,y[i]);
printf("\nAre the data ok?(Y/N)\n");
ch=getch();
if(ch=='n')
goto input;
//求线性方程系数
k[0][0]=fsum(x,x,limit);
k[0][1]=fsum(x,mi,limit);
k[0][2]=-fsum(x,y,limit);
k[1][0]=fsum(x,mi,limit);
k[1][1]=limit;
k[1][2]=-fsum(mi,y,limit);
//输出线性方程系数
printf("\nThe modulus is:\n");
for(i=0; i
{
for(j=0;j
printf("%15lf",k[i][j]);
printf("\n");
}
m0=fmatrix(0,1);
m1=fmatrix(1,2);
m2=fmatrix(2,0);
printf("\n%lf %lf %lf\n",m0,m1,m2);
if(m0==0)
printf("An error has occured! Matrix0 is zero!");//分母上的线性方程矩阵为零
else
printf("The function should be:\nY= %lf X %+lf\n",m1/m0,m2/m0);
}
回答者:hwsnet - 举人 四级 12-22 15:00
提问者对于答案的评价:
程序没什么用,但还是谢谢你。
评价已经被关闭 目前有 3 个人评价
好
33% (1)
不好 66% (2) 最小二乘法线性回归 最小二乘法线性回归
点击数:11604 发布日期:2005-9-23 21:28:00
【收藏】 【评论】 【打印】 【编程爱好者论坛】 【关闭】
最小二乘法线性回归
点击数:11604 发布日期:2005-9-23 21:28:00
【收藏】 【评论】 【打印】 【编程爱好者论坛】 【关闭】
最小二乘法线性回归
点击数:11604 发布日期:2005-9-23 21:28:00
【收藏】 【评论】 【打印】 【编程爱好者论坛】 【关闭】
最小二乘法线性回归
点击数:11604 发布日期:2005-9-23 21:28:00
【收藏】 【评论】 【打印】 【编程爱好者论坛】 【关闭】
最小二乘法线性回归
点击数:11604 发布日期:2005-9-23 21:28:00
【收藏】 【评论】 【打印】 【编程爱好者论坛】 【关闭】
Tag:数学 算法
在现实中经常遇到这样的问题,一个函数并不是以某个数学表达式的形式给出,而是以一些自变量与因变量的对应表给出,老师讲课的时候举的个例子是犯罪人的身高和留下的脚印长,可以测出一些人的数据然后得到一张表,它反应的是一个函数,回归的意思就是将它还原成数学表达式,这个式子也称为经验表达式,之所以叫经验就是说它不完全是实际中的那样准确,是有一定偏差的,只是偏差很小罢了。
最小二乘法
设经验方程是y=F(x),方程中含有一些待定系数an ,给出真实值{(xi,yi)|i=1,2,...n},将这些x, y 值代入方程然后作差,可以描述误差:yi-F(xi),为了考虑整体的误差,可以取平方和,之所以要平方是考虑到误差可正可负直接相加可以相互抵消,所以记误差为:
e=∑(yi-F(xi))^2
它是一个多元函数,有an 共n 个未知量,现在要求的是最小值。所以必然满足对各变量的偏导等于0,于是得到n 个方程:
de/da1=0
2007-12-10 14:24:27 分类:VC 学习评论阅读(313) 最小二乘法曲线拟合 //最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
最小二乘法曲线拟合 //最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
最小二乘法曲线拟合 //最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
//最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
//最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C
DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
//最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
//最小二乘法曲线拟合
typedef CArrayCDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,C DoubleArray *A)
{
//X,Y -- X,Y 两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//A -- 结果参数 改进的移动最小二乘法
改进的移动最小二乘法 陈美娟 程玉民
下载全文
[全文大小:163 K]
移动用户 编辑短信 020 到 [1**********]61 | 联通用户 点击链接 联通手机支付通道 进行支付 获得全文提取码!(2.00 元/篇)
下载安装 PDF 阅读器
支付说明 收不到提取码不会产生扣费。 其他的支付方式 包括神州行充值卡、支付宝、银行卡、电汇等。
国际标准刊号:ISSN 0254-0053
国内统一刊号:CN 31-1829
摘 要:
近年来发展的无网格方法大多采用移动员小二乘法来构造试函数,而应用移动最小二乘法形成的方程组有时会是病态的甚至奇异的,从而限制了它的发展和应用。本文采用带权正交函数作为基函数对移动最小二乘法做了改进,避免出现病态方程组,且在计算过程中不需要进
行短阵求逆运算,提高了计算速度。之后,借鉴牛顿法、平衡法和摄动法对由移动最小二乘
法得到的非线性代数方程组提出了新的求解方法。
将文章收藏到好诶
关 键 词:
移动最小二乘法 带权的正交基函数 改进的移动最小二乘法 非线性代数方程组 已解决问题 - 浏览941次
下一个已解决问题
找寻他的味道
穷酸秀才
最小二乘法怎么回事?
最小二乘法怎么回事?
有没有计算软件啊?
1年前 - 2个回答 - 15 收藏 寄给朋友 检举
添加/查看意见(0)
还可输入300个字
请输入上图中的验证码,字母不区分大小写。
点击查看更多 乘法 回事 最小 相关信息
雅虎知识堂连环画PK 大赛,手机、音箱、宝伞等你拿!
进士
最佳答案 - 由提问者1年前选出
最小二乘法
在我们研究两个变量(x, y) 之间的相互关系时,通常可以得到一系列成对的数据(x1, y1、x2, y2... xm , ym);将这些数据描绘在x -y直角座标系中(如图1), 若发现这些点在一条直线附近,可以令这条直线方程如(式1-1) 。
Y 计= a0 + a1 X (式1-1)
其中:a0、a1 是任意实数
为建立这直线方程就要确定a0和a1,应用《最小二乘法原理》,将实测值Yi 与利用(式1-1) 计算值(Y计=a0+a1X)的离差(Yi-Y计) 的平方和〔∑(Yi - Y 计)2〕最小为“优化判据”。
令: φ = ∑(Yi - Y计)2 (式1-2)
把(式1-1) 代入(式1-2) 中得:
φ = ∑(Yi - a0 - a1 Xi)2 (式1-3)
当∑(Yi-Y 计) 平方最小时,可用函数 φ 对a0、a1求偏导数,令这两个偏导数等于零。
(式1-4)
(式1-5)
亦即:
m a0 + (∑Xi ) a1 = ∑Yi (式1-6)
(∑Xi ) a0 + (∑Xi2 ) a1 = ∑(Xi, Yi) (式1-7)
得到的两个关于a0、 a1为未知数的两个方程组,解这两个方程组得出:
a0 = (∑Yi) / m - a1(∑Xi) / m (式1-8)
a1 = [∑Xi Yi - (∑Xi ∑Yi)/ m] / [∑Xi2 - (∑Xi)2 / m)] (式1-9)
这时把a0、a1代入(式1-1) 中, 此时的(式1-1) 就是我们回归的元线性方程即:数学模型。
在回归过程中,回归的关联式是不可能全部通过每个回归数据点(x1, y1、 x2, y2...xm,ym), 为了判断关联式的好坏, 可借助相关系数“R”,统计量“F”,剩余标准偏差“S”进行判断;“R”越趋近于 1 越好;“F”的绝对值越大越好;“S”越趋近于 0 越好。
R = [∑XiYi - m (∑Xi / m)(∑Yi / m)]/ SQR{[∑Xi2 - m (∑Xi / m)2][∑Yi2 - m (∑Yi / m)2]} (式1-10) *
在(式1-1) 中,m 为样本容量,即实验次数;Xi 、Yi 分别任意一组实验X 、Y 的数值。微积分应用课题一 最小二乘法
从前面的学习中, 我们知道最小二乘法可以用来处理一组数据, 可以从一组测定的数据中寻求变量之间的依赖关系, 这种函数关系称为经验公式. 本课题将介绍最小二乘法的精确定义及如何寻求 与 之间近似成线性关系时的经验公式. 假定实验测得变量之间的 个数据 , , „, , 则在 平面上, 可以得到 个点 , 这种图形称为“散点图”, 从图中可以粗略看出这些点大致散落在某直线近旁, 我们认为 与 之间近似为一线性函数, 下面介绍求解步骤.
考虑函数 , 其中 和 是待定常数. 如果 在一直线上, 可以认为变量之间的关系为 . 但一般说来, 这些点不可能在同一直线上. 记 , 它反映了用直线 来描述 , 时, 计算值 与实际值 产生的偏差. 当然要求偏差越小越好, 但由于 可正可负, 因此不能认为总偏差 时, 函数 就很好地反映了变量之间的关系, 因为此时每个偏差的绝对值可能很大. 为了改进这一缺陷, 就考虑用 来代替 . 但是由于绝对值不易作解析运算, 因此, 进一步用 来度量总偏差. 因偏差的平方和最小可以保证每个偏差都不会很大. 于是问题归结为确定 中的常数 和 , 使 为最小. 用这种方法确定系数 , 的方法称为最小二乘法.
由极值原理得 , 即
解此联立方程得
(*)
问题 I 为研究某一化学反应过程中, 温度 ℃)对产品得率 (%) 的影响, 测得数据如下:
温度 ℃)
100 110 120 130 140 150 160 170 180 190
得率 (%)
45 51 54 61 66 70 74 78 85 89
(1) 利用“ListPlot”函数, 绘出数据 的散点图(采用格式:
ListPlot[{ , , „, }, Prolog->AbsolutePointSize[3]] );
(2) 利用“Line”函数, 将散点连接起来, 注意观察有何特征? (采用格式: Show[Graphics[Line[{ , , „, }]] , Axes->True ]) ;
(3) 根据公式(*), 利用“Apply”函数及集合的有关运算编写一个小的程序, 求经验公式 ;
(程序编写思路为: 任意给定两个集合A (此处表示温度) 、B(此处表示得率), 由公式(*)可定义两个二元函数(集合A 和B 为其变量) 分别表示 和 . 集合A 元素求和: Apply[Plus,A] 表示将加法施加到集合A 上, 即各元素相加, 例如
Apply[Plus,{1,2,3}]=6;Length[A]表示集合A 元素的个数, 即为n; A.B表示两集合元素相乘相加;A*B表示集合A 与B 元素对应相乘得到的新的集合.)
(4) 在同一张图中显示直线 及散点图;
(5) 估计温度为200时产品得率.
然而, 不少实际问题的观测数据 , , „, 的散点图明显地不能用线性关系来描叙, 但确实散落在某一曲线近旁, 这时可以根据散点图的轮廓和实际经验, 选一条曲线来近似表达 与 的相互关系.
问题 II 下表是美国旧轿车价格的调查资料, 今以 表示轿车的使用年数, (美元) 表示相应的平均价格, 求 与 之间的关系.
使用年数
1 2 3 4 5 6 7 8 9 10
平均价格
2651 1943 1494 1087 765 538 484 290 226 204
(1) 利用“ListPlot”函数绘出数据 的散点图, 注意观察有何特征?
(2) 令 , 绘出数据 的散点图, 注意观察有何特征?
(3) 利用“Line”函数, 将散点 连接起来, 说明有何特征?
(4) 利用最小二乘法, 求 与 之间的关系;
(5) 求 与 之间的关系;
(6) 在同一张图中显示散点图 及 关于 的图形.
思考与练习
1. 假设一组数据 : , , „, 变量之间近似成线性关系, 试利用集合的有关运算, 编写一简单程序: 对于任意给定的数据集合 , 通过求解极值原理所包含的方程组, 不需要给出 、 计算的表达式, 立即得到 、 的值, 并就本课题 I /(3)进行实验.
注: 利用Transpose 函数可以得到数据A 的第一个分量的集合, 命令格式为:
先求A 的转置, 然后取第一行元素, 即为数据A 的第一个分量集合, 例如
(A即为矩阵 )
= (数据A 的第一个分量集合)
= (数据A 的第二个分量集合)
B-C 表示集合B 与C 对应元素相减所得的集合, 如 = .
2. 最小二乘法在数学上称为曲线拟合, 请使用拟合函数“Fit”重新计算 与 的值, 并与先前的结果作一比较.
注: Fit函数使用格式:
设变量为x, 对数据A 进行线性拟合, 如对题1中的A 拟合函数为:
=
参考资料:http://www.eaihua.com/ahrjyl.htm
www.math.sjtu.edu.cn/gskc/syal/最小二乘法 ... 94K
1年前
添加/查看意见(0)
还可输入300个字
请输入上图中的验证码,字母不区分大小写。
提问者对最佳答案的评价
Thank you!
其他回答(1)
高级书童
不知道
最小二乘法建模
最小二乘法建模的一般步骤是什么? 怎样才能得到最好的模型?
谢谢!
7个月前 - 2个回答 - 15 收藏 寄给朋友 检举
添加/查看意见(0)
还可输入300个字
请输入上图中的验证码,字母不区分大小写。
点击查看更多 乘法 建模 最小 相关信息
雅虎知识堂连环画PK 大赛,手机、音箱、宝伞等你拿!
大傻瓜
传胪
最佳答案 - 由投票者6个月前选出
呵呵,实际上就是说,通过合理设置参数,使得拟合点和实际点的距离的平方和的值最小。关键是要有好的建模假设。
7个月前
添加/查看意见(0)
还可输入300个字
请输入上图中的验证码,字母不区分大小写。
其他回答(1)
榜眼
百度看看
最小二乘法多次曲线拟合
来源:CSDN 发布会员:新书城收集整理 发布时间:2006-9-17 人气:709
因本人项目用到最小二乘法多次曲线拟合,便抽时间写了一个实现的程序,不敢独享,如下:
Option Explicit
'****************************************************************************************************
'** X()------Double 实型一维数组,长度为 n 。存放给定 n 个数据点的 X 坐标。 **
'** Y()------Double 实型一维数组,长度为 n 。存放给定 n 个数据点的 Y 坐标。 **
'** n -------Integer 变量。给定数据点的个数。 **
'** a()------Double 实型一维数组,长度为 m 。返回 m-1 次拟合多项式的 m 个系数。 **
'** m -------Integer 变量。拟合多项式的项数,即拟合多项式的最高次数为 m-
1。 **
'** 要求 mn 或 m>20 ,则本函数自动按 m=min{n,20} 处理。 **
'**rdblAverageX--Double 变量,返回给定n 个数据点的 X 坐标的平均值 ** '** dt()------Double 实型一维数组,长度为 3。其中: **
'** dt(0) 返回拟合多项式与数据点误差的平方和; **
'** dt(1) 返回拟合多项式与数据点误差的绝对值之和; **
'** dt(2) 返回拟合多项式与数据点误差绝对值的最大值。 **
'*****************************************************************************************************
Public Sub Iapcir(X() As Double, _
Y() As Double, _
ByVal n As Integer, _
ByRef a() As Double, _
ByVal m As Integer, _
ByRef rdblAverageX As Double, _
ByRef dt() As Double)
Dim I As Integer, J As Integer, K As Integer
Dim Z As Double, P As Double, C As Double, G As Double, Q As Double, D1 A s Double, D2 As Double
Dim S(19) As Double, T(19) As Double, B(19) As Double
For I = 0 To m - 1
a(I) = 0
Next I
If m > n Then m = n
If m > 20 Then m = 20
Z = 0#
For I = 0 To n - 1
rdblAverageX = rdblAverageX + X(I)
Z = Z + X(I) / (1# * n)
Next I
rdblAverageX = rdblAverageX / n
B(0) = 1#
D1 = 1# * n
P = 0#
C = 0#
For I = 0 To n - 1
P = P + (X(I) - Z)
C = C + Y(I)
Next I
C = C / D1
P = P / D1
a(0) = C * B(0)
If m > 1 Then
T(1) = 1#
T(0) = (-1) * P
D2 = 0#
C = 0#
G = 0#
For I = 0 To n - 1
Q = X(I) - Z - P
D2 = D2 + Q * Q
C = C + Y(I) * Q
G = G + (X(I) - Z) * Q * Q
Next I
C = C / D2
P = G / D2
Q = D2 / D1
D1 = D2
a(1) = C * T(1)
a(0) = C * T(0) + a(0)
End If
For J = 2 To m - 1
S(J) = T(J - 1)
S(J - 1) = (-1) * P * T(J - 1) + T(J - 2)
If J >= 3 Then
For K = J - 2 To 1 Step -1
S(K) = (-1) * P * T(K) + T(K - 1) - Q * B(K) Next K
End If
S(0) = (-1) * P * T(0) - Q * B(0)
D2 = 0#
C = 0#
G = 0#
For I = 0 To n - 1
Q = S(J)
For K = J - 1 To 0 Step -1
Q = Q * (X(I) - Z) + S(K)
Next K
D2 = D2 + Q * Q
C = C + Y(I) * Q