1. 最小二乘拟合原理

参考西安交大数值分析教材

当f(x)是定义在点集X上的列表函数时，内积为

2. 最小二乘拟合问题的求解过程

最小二乘拟合问题的求解过程有三种方法

法一：

（1）根据问题的特点选择一组线性无关的基函数，

（2）解方正规方程组

（3）求得拟合函数

法二：

（1）利用三项递推关系式构造正交多项式

（2）由式5.2.11得到最小二乘拟合多项式

法三：

（1）做变量替换，利用已有的正交基函数构造最优平方逼近函数

（2）转到法二求解

3. python程序以及说明

本程序实现了上述的法一和法二，对于法一，程序选择了xk，(k=1，2，3 ……)作为基函数，然后使用改进平方根法求解正规方程组得到系数，最后利用系数构造拟合函数。对于法二，程序中使用最高项系数为1的正交多项式的三项递推关系构造正交基函数，然后通过式5.2.11求解系数。正交多项式的三项递推关系如下

python代码如下

import numpy as np
from sympy.abc import t
from sympy import simplifyclass LeastSquareFit():'''用于计算列表函数的最小二乘拟合输入参数：X：自变量列表Y：函数值列表n：多项式阶数JiFunction：基函数形式'''def __init__(self,X,Y,n,JiFunction='正交多项式'):self.x=np.array(X)self.y=np.array(Y)self.n=nself.w=self.myset_w(f=1) #设置权函数self.JiFunction=self.generate_JiFuction(n,JiFunction) #生成基函数self.c=[0]*(n+1)  #初始化拟合多项式的系数#计算权函数def myset_w(self,f):if type(f) is int or float:y=[f]*len(self.x)else:y=[f.subs(t,x) for x in self.x]return y#生成正交基函数def generate_JiFuction(self,m,function_name):if function_name == '正交多项式':value_num=mg={}g[0]=1beta0=np.dot([(t*g[0]).subs(t,x) for x in self.x] , np.array(self.w)*([g[0]]*len(self.x)))gama0=np.dot([g[0]]*len(self.x) , np.array(self.w)*([g[0]]*len(self.x)))g[1]=t-beta0/gama0for k in range(1,value_num):beta=np.dot([(t*g[k]).subs(t,x) for x in self.x] , np.array(self.w)*[g[k].subs(t,x) for x in self.x])gama=np.dot([g[k].subs(t,x) for x in self.x] , np.array(self.w)*[g[k].subs(t,x) for x in self.x])if k==1:gama_pre=np.dot([g[k-1]]*len(self.x) , np.array(self.w)*([g[k-1]]*len(self.x)))else:gama_pre=np.dot([g[k-1].subs(t,x) for x in self.x] , np.array(self.w)*[g[k-1].subs(t,x) for x in self.x])b=beta/gamac=gama/gama_preg[k+1]=simplify((t-b)*g[k]-c*g[k-1])self.g=g#生成线性无关的x的k次方函数elif function_name == 'x^k': value_num=mg={}g[0]=1g[1]=tfor k in range(1,value_num):g[k+1]=t**(k+1)self.g=g#正交基函数正规方程组求解def ZhengJiao_solve_ck(self):for k in range(self.n+1):if k==0:c1=np.dot([self.g[k]]*len(self.x) , np.array(self.w)*self.y)c2=np.dot([self.g[k]]*len(self.x) , np.array(self.w)*([self.g[k]]*len(self.x)))else:c1=np.dot([self.g[k].subs(t,x) for x in self.x] , np.array(self.w)*self.y)c2=np.dot([self.g[k].subs(t,x) for x in self.x] , np.array(self.w)*[self.g[k].subs(t,x) for x in self.x])self.c[k]=c1/c2#最小二乘多项式正规方程组求解def ZXEC_solve_ck(self):#计算A和bG_pre=[]for key in self.g.keys():if key==0:passelse:G_pre.append(self.g[key])G1=[[g.subs(t,x) for g in G_pre] for x in self.x]one_temp=np.ones(len(self.x))one_temp=one_temp.reshape(len(self.x),1)G2=np.c_[one_temp,G1]G=np.array(G2)W=np.diag(self.w)A_temp=np.dot(G.T,W)A=np.dot(A_temp,G)b=np.dot(np.dot(G.T,W),self.y)#LDLT分解L,U=self.my_LU(A)D=np.diag(U)D=np.diag(D)#解yy={}y[0]=b[0]for i in range(1,self.n+1):Ly=[L[i][k]*y[k] for k in range(0,i)]sum_Ly=sum(Ly)y[i]=b[i]-sum_Ly#解xx={}x[self.n]=y[self.n]/U[self.n][self.n]for j in range(self.n-1,-1,-1):Ux=[U[j][k]*x[k] for k in range(j+1,self.n+1)]sum_Ux=sum(Ux)x[j]=(y[j]-sum_Ux)/U[j][j]#将x变为cfor i in range(len(x)):self.c[i]=x[i]#LU分解def my_LU(self,A):A=np.asarray(A)U=np.zeros_like(A)L=np.zeros_like(A)U[0,:]=A[0,:]L[1:A.shape[0],0]=A[1:A.shape[0],0]/U[0][0]for i in range(1,A.shape[0]-1):U[i][i]=A[i][i]-sum([L[i][k]*U[k][i] for k in range(0,i)])for j in range(i+1,A.shape[1]):U[i][j]=A[i][j]-sum([L[i][k]*U[k][j] for k in range(0,i)])L[j][i]=(A[j][i]-sum([L[j][k]*U[k][i] for k in range(0,i)]))/U[i][i]U[A.shape[0]-1][A.shape[1]-1]=A[A.shape[0]-1][A.shape[1]-1]-\sum([L[A.shape[0]-1][k]*U[k][A.shape[1]-1] for k in range(0,A.shape[1]-1)])for i in range(L.shape[0]):L[i][i]=1return L,U   #计算逼近多项式def solve_pn(self):pn=[]for k in range(self.n+1):pn.append(simplify(self.c[k]*self.g[k]))result=0for p in pn:result=result+preturn simplify(result)#计算误差def LeastSquareFit_error(self):y_pn=self.solve_pn()e=0for index,value in enumerate(self.x):temp=abs(y_pn.subs(t,value)-self.y[index])e=e+temp**2return pow(e,0.5)######################测试程序####################
if __name__ == "__main__":#------------------使用正交多项式---------------------# X=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]# Y=[5.1234, 5.3057, 5.5687, 5.9375, 6.4370, 7.0978, 7.9493, 9.0253, 10.3627]# test=LeastSquareFit(X,Y,4,JiFunction='正交多项式')# test.myset_w(f=1) #设置权函数# test.ZhengJiao_solve_ck() #求解系数c# print("多项式系数", test.c)# result=test.solve_pn() #求解拟合的多项式# print("拟合多项式", result)# my_error=test.LeastSquareFit_error() #计算误差# print("误差为：", my_error)#------------------使用x^k---------------------------X = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]Y = [5.1234, 5.3057, 5.5687, 5.9375, 6.4370, 7.0978, 7.9493, 9.0253, 10.3627]test=LeastSquareFit(X,Y,4,JiFunction='x^k')test.myset_w(f=1) #设置权函数test.ZXEC_solve_ck() #求解多项式系数cprint("多项式系数", test.c)result=test.solve_pn() #求解拟合的多项式print("拟合多项式", result)my_error=test.LeastSquareFit_error() #计算误差print("误差为：", my_error)

4. 程序使用说明

首先给定输入数据X和Y，然后确定多项式阶数比如4，最后选择基函数，可选“正交基函数”或“x^k”。确定完参数后，具体使用步骤为：（1）类实例化；（2）调用myset_w()设置权函数；（3）调用ZXEC_solve_ck()（或ZhengJiao_solve_ck()）求解多项式系数；（4）调用solve_pn()生成拟合多项式；（5）调用LeastSquareFit_error()计算误差。

最小二乘拟合问题求解算法（含python代码）相关推荐

12种降维方法终极指南（含Python代码）
12种降维方法终极指南(含Python代码) 你遇到过特征超过1000个的数据集吗?超过5万个的呢?我遇到过.降维是一个非常具有挑战性的任务,尤其是当你不知道该从哪里开始的时候.拥有这么多变量既是一个 ...
一文看懂Stacking！(含Python代码)
一文看懂Stacking!(含Python代码) https://mp.weixin.qq.com/s/faQNTGgBZdZyyZscdhjwUQ 转载于:https://www.cnblogs.c ...
感知机算法之Python代码实现
感知机算法之Python代码实现 1.算法简介感知机学习算法原始形式: 输入:训练集T 输出:w,b 感知机模型:f(x)=sign(w·x+b) 算法步骤: 1.初始化参数w0,b0 2.在训练集 ...
BPR贝叶斯个性化推荐算法—推荐系统基础算法（含python代码实现以及详细例子讲解）
BPR贝叶斯个性化排序算法一.问题导入二.显示反馈与隐式反馈 2.1 显式反馈与隐式反馈基本概念 2.2 显式反馈与隐式反馈的比较 2.3 显式反馈与隐式反馈的评价方法 2.3.1 显式反馈数据模 ...
密码学——保序加密算法（OPE算法-2009年提出）通俗易懂解析（小学生都能懂！）含python代码
一. 预备知识保序加密算法:最初是由2009年,Boldyreva等四个人提出来的,可简称BCLO-09算法,论文题目为<Order-Preserving Symmetric Encrypti ...
如何用机器学习算法来进行电影分类?(含Python代码)
电影分析--K近邻算法周末,小迪与女朋友小西走出电影院,回味着刚刚看过的电影. 小迪:刚刚的电影很精彩,打斗场景非常真实,又是一部优秀的动作片! 小西:是吗?我怎么感觉这是一部爱情片呢?真心被男主女 ...
图解匈牙利算法（含python代码）
文章目录分配问题匈牙利算法算法步骤算法实现 python版本 C++版本分配问题分配问题/指派问题(Assignment Problem)作为线性规划问题的一个特例,在运筹学研究中占有重要 ...
算法笔记（11）逻辑回归算法及Python代码实现
逻辑回归算法是一种被广泛使用的分类算法,通过训练数据中的正负样本,学习样本特征到样本标签之间的假设函数.逻辑回归假设因变量 y 服从伯努利分布,而线性回归假设因变量 y 服从高斯分布. 因此与线性回归 ...
决策树模型 - （ID3算法、C4.5算法） - Python代码实现
目录算法简介信息熵(Entropy) 信息增益(Information gain) - ID3算法信息增益率(gain ratio) - C4.5算法源数据代码实现 - ID3算法代码实现 ...
维吉尼亚密码的破解算法及python代码实现
目录 1. 密文描述 1.1 密文1 1.2 密文2 2. 破解原理 2.1 重合指数法确定密钥长度 2.2 互重合指数确定子串间相对偏移 2.3 密钥字的确定 2.4 密文破解 3. 破解代码参考 ...

最小二乘拟合问题求解算法（含python代码）

1. 最小二乘拟合原理

2. 最小二乘拟合问题的求解过程

3. python程序以及说明

4. 程序使用说明

最小二乘拟合问题求解算法（含python代码）相关推荐

最新文章

热门文章