Python计算&绘图——曲线拟合问题
2014-03-16 11:19
471 查看
题目来自老师的课后作业,如下所示。很多地方应该可以直接调用函数,但是初学Python,对里面的函数还不是很了解,顺便带着学习的态度,尽量自己动手code。
![](http://img.blog.csdn.net/20140316112532062?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvenV5dWFuemh1/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
测试版代码,里面带有很多注释和测试代码:
M=3时的运行结果:
Figure(1):
![](http://img.blog.csdn.net/20140316113410734?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvenV5dWFuemh1/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
Figure(2):
![](http://img.blog.csdn.net/20140316113422015?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvenV5dWFuemh1/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
初次编写这么长的代码,思路不是有一点的混乱
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/sad.gif)
。其中有
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/fastcry.gif)
也有
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/proud.gif)
。以后会继续来优化这个程序,作为学习Python的入口。
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/struggle.gif)
测试版代码,里面带有很多注释和测试代码:
# -*- coding: cp936 -*- import math import random import matplotlib.pyplot as plt import numpy as np ''' 在x=[0,1]上均匀采样10个点组成一个数据集D=[a,b] ''' a = [] b = [] x=0 def func(x): mu=0 sigma=0.1 epsilon = random.gauss(mu,sigma) #高斯分布随机数 return np.sin(2*np.pi*x)+epsilon for i in range(0,10): x=x+1.0/11.0 a.append(x) b.append(func(x)) #定义输出矩阵函数 def print_matrix( info, m ): i = 0; j = 0; l = len(m) print info for i in range( 0, len( m ) ): for j in range( 0, len( m[i] ) ): if( j == l ): print ' |', print '%6.4f' % m[i][j], print print #定义交换变量函数 def swap( a, b ): t = a; a = b; b = t #定义线性方程函数,高斯消元法 def solve( ma, b, n ): global m; m = ma # 这里主要是方便最后矩阵的显示 global s; i = 0; j = 0; row_pos = 0; col_pos = 0; ik = 0; jk = 0 mik = 0.0; temp = 0.0 n = len( m ) # row_pos 变量标记行循环, col_pos 变量标记列循环 while( ( row_pos < n ) and( col_pos < n ) ): # 选主元 mik = - 1 for i in range( row_pos, n ): if( abs( m[i][col_pos] ) > mik ): mik = abs( m[i][col_pos] ) ik = i if( mik == 0.0 ): col_pos = col_pos + 1 continue # 交换两行 if( ik != row_pos ): for j in range( col_pos, n ): swap( m[row_pos][j], m[ik][j] ) swap( m[row_pos] , m[ik] ); try: # 消元 m[row_pos] /= m[row_pos][col_pos] except ZeroDivisionError: # 除零异常 一般在无解或无穷多解的情况下出现…… return 0; j = n - 1 while( j >= col_pos ): m[row_pos][j] /= m[row_pos][col_pos] j = j - 1 for i in range( 0, n ): if( i == row_pos ): continue m[i] -= m[row_pos] * m[i][col_pos] j = n - 1 while( j >= col_pos ): m[i][j] -= m[row_pos][j] * m[i][col_pos] j = j - 1 row_pos = row_pos + 1; col_pos = col_pos + 1 for i in range( row_pos, n ): if( abs( m[i] ) == 0.0 ): return 0 return 1 matrix_A=[] #将系数矩阵A的所有元素存到a[n-1][n-1]中 matrix_b=[] X=a Y=b N=len(X) M=3 #对于题目中要求的不同M[0,1,3,9]值,需要在这里更改,然后重新编译运行 #计算线性方程组矩阵A的第[i][j]个元素A[i][j] def matrix_element_A(x,i,j,n): sum_a=0 for k in range(0,n): sum_a = sum_a+pow(x[k],i+j-2) #x[0]到x[n-1],共n个元素求和 return sum_a for i in range(0,M+1): matrix_A.append([]) for j in range(0,M+1): matrix_A[i].append(0) matrix_A[i][j] = matrix_element_A(X,i+1,j+1,N) #计算线性方程组矩阵b的第[i]行元素b[i] def matrix_element_b(x,y,i,n): sum_b=0 for k in range(0,n): sum_b=sum_b+y[k]*pow(x[k],i-1) #x[0]到x[n-1],共n个元素求和 return sum_b for i in range(0,M+1): matrix_b.append(matrix_element_b(X,Y,i+1,N)) #函数matrix_element_A_()用来求扩展矩阵A_,array_A表示系数矩阵A,array_b表示方程组右侧常数,A_row表示A的行秩 def matrix_element_A_(array_A,array_b,A_row): M=A_row #局部变量M,与全局变量M无关 matrix_A_= [] for i in range(0,M+1): matrix_A_.append([]) for j in range(0,M+2): matrix_A_[i].append(0) if j<M+1: matrix_A_[i][j] = array_A[i][j] elif j==M+1: #如果不加这个控制条件,matrix_A_将被array_b刷新 matrix_A_[i][j] = array_b[i] return matrix_A_ matrix_A_ = matrix_element_A_(matrix_A,matrix_b,M) ''' 多项式拟合函数 ''' #x为自变量,w为多项式系数,m为多项式的阶数 def poly_fit(x,wp,m): sumf = 0 for j in range(0,m+1): sumf=sumf+wp[j]*pow(x,j) return sumf ''' sin(2*pi*x)在x=0处的3阶泰勒展开式 ''' coef_taylor = [] #正弦函数的泰勒展开式系数 K=3 #展开到K阶 if K%2==0: print "K必须为正奇数" s = 0 k=(K-1)/2+1 #小k为系数个数 #求K阶泰勒展开式的系数: for i in range(0,k): s = pow(-1,i)*pow(2*np.pi,2*i+1)/math.factorial(2*i+1) coef_taylor.append(s) print "%d阶泰勒级数展开式的系数为:" %K print coef_taylor #tx为泰勒展开式函数的自变量 def sin_taylor(tx): sum_tay=0 for i in range(0,k): sum_tay=sum_tay+coef_taylor[i]*pow(tx,2*k+1) return sum_tay poly_taylor_a = [] #泰勒展开式函数的输入值 poly_taylor_b = [] #泰勒展开式函数的预测值 for i in range(0,N): poly_taylor_a.append(a[i]) poly_taylor_b.append(sin_taylor(poly_taylor_a[i])) ''' 在x=[0,1]上生成100个点,作为测试集 ''' testa = [] #测试集的横坐标 testb = [] #测试集的纵坐标 x=0 for i in range(0,100): x=x+1.0/101.0 testa.append(x) testb.append(np.sin(2*np.pi*x)) ''' 计算泰勒展开式模型的训练误差和测试误差 ''' #定义误差函数: #ly为真实值,fx为预测值 def Lfun(ly,fx): L=0 for i in range(0,len(fx)): L=L+pow(ly[i]-fx[i],2) return L ''' 主程序 ''' if __name__ == '__main__': # 求解方程组, 并输出方程组的可解信息 ret = solve( matrix_A_, 0, 0 ) if( ret== 0 ): print "方 程组无唯一解或无解\n" # 输出方程组及其解,解即为w[j] w = [] for i in range( 0, len( m ) ): w.append(m[i][len( m )]) print "M=%d时的系数w[j]:" %M print w #多项式拟合后的预测值: poly_a = [] poly_b = [] for i in range(0,N): poly_a.append(a[i]) poly_b.append(poly_fit(poly_a[i],w,M)) #fxtay为泰勒展开式的预测值,LCtaylor为测试误差: fxtay = [] for i in range(0,100): fxtay.append(sin_taylor(testa[i])) LCtaylor = Lfun(testb,fxtay)/100 print "三阶泰勒展开式的测试误差为:%f" %LCtaylor #fxpoly为M阶多项式拟合函数的预测值,LXpoly为训练误差: fxpoly = [] for i in range(0,N): #len(poly_b)=N=10 fxpoly.append(poly_fit(a[i],w,M)) LXpoly = Lfun(b,fxpoly)/len(poly_b) print "M=%d时多项式拟合函数的训练误差为:%f" % (M,LXpoly) #fxpolyc为M阶多项式拟合函数的预测值,LCpoly为测试误差: fxpolyc = [] for i in range(0,100): fxpolyc.append(poly_fit(testa[i],w,M)) LCpoly = Lfun(testb,fxpolyc)/100 print "M=%d时多项式拟合函数的测试误差为:%f" % (M,LCpoly) #多项式拟合的效果: fig1 = plt.figure(1) plt.plot(poly_a,poly_b,color='blue',linestyle='solid',marker='o') #加入epsilon后的样本: plt.plot(a,b,color='red',linestyle='dashed',marker='x') #泰勒展开式拟合效果: plt.plot(poly_taylor_a,poly_taylor_b,color='yellow',linestyle='dashed',marker='o') #figure(2)对比多项式拟合函数与训练数据: fig2 = plt.figure(2) plt.plot(poly_a,poly_b,color='blue',linestyle='solid',marker='o') plt.plot(a,b,color='red',linestyle='dashed',marker='x') plt.show()
M=3时的运行结果:
3阶泰勒级数展开式的系数为: [6.283185307179586, -41.341702240399755] M=3时的系数w[j]: [-0.28492708632293295, 13.031310645420685, -37.730992850050448, 25.464782221275197] 三阶泰勒展开式的测试误差为:100.889335 M=3时多项式拟合函数的训练误差为:0.008933 M=3时多项式拟合函数的测试误差为:0.007886
Figure(1):
Figure(2):
初次编写这么长的代码,思路不是有一点的混乱
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/sad.gif)
。其中有
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/fastcry.gif)
也有
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/proud.gif)
。以后会继续来优化这个程序,作为学习Python的入口。
![](http://static.blog.csdn.net/xheditor/xheditor_emot/default/struggle.gif)
相关文章推荐
- Python计算&绘图——曲线拟合问题(转)
- python科学计算学习二:matplotlib绘图,图标注释(2)
- Python 网页爬虫 & 文本处理 & 科学计算 & 机器学习 & 数据挖掘兵器谱
- Stanford机器学习---第三讲. 逻辑回归和过拟合问题的解决 logistic Regression & Regularization
- python绘图工具包 matplotlib 中文乱码问题
- [Python Tip]如何计算时间差
- Python基于最小二乘法实现曲线拟合示例
- 解决多个版本的python共存时的问题 => 持续更新
- python matplotlib绘图时图例显示问题
- 【Python开发】matplotlib绘图不显示问题解决plt.show()
- python第三方库安装问题-'ascii' codec can't decode byte 0xb0
- 解决Python代码编码问题 SyntaxError: Non-UTF-8 code starting with '\xc1'
- Stanford机器学习---第三讲. 逻辑回归和过拟合问题的解决 logistic Regression & Regularization
- Python r”\" SyntaxError 问题分析
- matlap 指数曲线拟合及其R2的计算
- Python数值计算:一 使用Pylab绘图(1)
- Python实现二维曲线拟合
- python编码问题 decode('unicode-escape')
- Python 网页爬虫 & 文本处理 & 科学计算 & 机器学习 & 数据挖掘兵器谱(转)
- python matplotlib绘图时图例显示问题