Python计算&绘图——曲线拟合问题(转)

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Python计算&绘图——曲线拟合问题(转)相关的知识,希望对你有一定的参考价值。

题目来自老师的课后作业,如下所示。很多地方应该可以直接调用函数,但是初学Python,对里面的函数还不是很了解,顺便带着学习的态度,尽量自己动手code。

         技术分享

 

测试版代码,里面带有很多注释和测试代码:

 

[python] view plain copy
 
 技术分享技术分享
  1. # -*- coding: cp936 -*-  
  2. import math  
  3. import random  
  4. import matplotlib.pyplot as plt  
  5. import numpy as np  
  6.   
  7.   
  8. ‘‘‘‘‘ 
  9. 在x=[0,1]上均匀采样10个点组成一个数据集D=[a,b] 
  10. ‘‘‘  
  11. a = []  
  12. b = []  
  13. x=0  
  14. def func(x):  
  15.     mu=0  
  16.     sigma=0.1  
  17.     epsilon = random.gauss(mu,sigma) #高斯分布随机数  
  18.     return np.sin(2*np.pi*x)+epsilon  
  19. for i in range(0,10):  
  20.     x=x+1.0/11.0  
  21.     a.append(x)  
  22.     b.append(func(x))  
  23.   
  24.   
  25.   
  26.   
  27. #定义输出矩阵函数  
  28. def print_matrix( info, m ):   
  29.     i = 0; j = 0; l = len(m)  
  30.     print info  
  31.     for i in range( 0, len( m ) ):  
  32.         for j in range( 0, len( m[i] ) ):  
  33.             if( j == l ):  
  34.                 print ‘ |‘,  
  35.             print ‘%6.4f‘ % m[i][j],  
  36.         print  
  37.     print  
  38.   
  39.   
  40. #定义交换变量函数  
  41. def swap( a, b ):  
  42.     t = a; a = b; b = t  
  43.       
  44. #定义线性方程函数,高斯消元法      
  45. def solve( ma, b, n ):  
  46.     global m; m = ma # 这里主要是方便最后矩阵的显示  
  47.     global s;      
  48.     i = 0; j = 0; row_pos = 0; col_pos = 0; ik = 0; jk = 0  
  49.     mik = 0.0; temp = 0.0    
  50.     n = len( m )  
  51. # row_pos 变量标记行循环, col_pos 变量标记列循环  
  52.     while( ( row_pos < n ) and( col_pos < n ) ):  
  53.         # 选主元  
  54.         mik = - 1  
  55.         for i in range( row_pos, n ):  
  56.             if( abs( m[i][col_pos] ) > mik ):  
  57.                 mik = abs( m[i][col_pos] )  
  58.                 ik = i  
  59.         if( mik == 0.0 ):  
  60.             col_pos = col_pos + 1  
  61.             continue  
  62.         # 交换两行  
  63.         if( ik != row_pos ):  
  64.             for j in range( col_pos, n ):  
  65.                 swap( m[row_pos][j], m[ik][j] )  
  66.                 swap( m[row_pos][n], m[ik][n] );    
  67.         try:  
  68.             # 消元  
  69.             m[row_pos][n] /= m[row_pos][col_pos]  
  70.         except ZeroDivisionError:  
  71.             # 除零异常 一般在无解或无穷多解的情况下出现……  
  72.             return 0;       
  73.         j = n - 1  
  74.         while( j >= col_pos ):  
  75.             m[row_pos][j] /= m[row_pos][col_pos]  
  76.             j = j - 1  
  77.         for i in range( 0, n ):  
  78.             if( i == row_pos ):  
  79.                 continue  
  80.             m[i][n] -= m[row_pos][n] * m[i][col_pos]  
  81.             j = n - 1  
  82.             while( j >= col_pos ):  
  83.                 m[i][j] -= m[row_pos][j] * m[i][col_pos]  
  84.                 j = j - 1  
  85.         row_pos = row_pos + 1; col_pos = col_pos + 1  
  86.     for i in range( row_pos, n ):  
  87.         if( abs( m[i][n] ) == 0.0 ):  
  88.             return 0  
  89.     return 1  
  90.   
  91.   
  92.   
  93.   
  94. matrix_A=[]         #将系数矩阵A的所有元素存到a[n-1][n-1]中  
  95. matrix_b=[]  
  96. X=a  
  97. Y=b  
  98. N=len(X)  
  99. M=3    #对于题目中要求的不同M[0,1,3,9]值,需要在这里更改,然后重新编译运行  
  100.   
  101.   
  102. #计算线性方程组矩阵A的第[i][j]个元素A[i][j]      
  103. def matrix_element_A(x,i,j,n):   
  104.     sum_a=0  
  105.     for k in range(0,n):     
  106.         sum_a = sum_a+pow(x[k],i+j-2)   #x[0]到x[n-1],共n个元素求和  
  107.     return sum_a  
  108.   
  109.   
  110. for i in range(0,M+1):    
  111.     matrix_A.append([])    
  112.     for j in range(0,M+1):    
  113.         matrix_A[i].append(0)  
  114.         matrix_A[i][j] = matrix_element_A(X,i+1,j+1,N)  
  115. #计算线性方程组矩阵b的第[i]行元素b[i]  
  116. def matrix_element_b(x,y,i,n):   
  117.     sum_b=0  
  118.     for k in range(0,n):  
  119.         sum_b=sum_b+y[k]*pow(x[k],i-1)  #x[0]到x[n-1],共n个元素求和  
  120.     return sum_b  
  121. for i in range(0,M+1):  
  122.     matrix_b.append(matrix_element_b(X,Y,i+1,N))  
  123.   
  124.   
  125. #函数matrix_element_A_()用来求扩展矩阵A_,array_A表示系数矩阵A,array_b表示方程组右侧常数,A_row表示A的行秩  
  126. def matrix_element_A_(array_A,array_b,A_row):  
  127.     M=A_row  #局部变量M,与全局变量M无关  
  128.     matrix_A_= []  
  129.     for i in range(0,M+1):  
  130.         matrix_A_.append([])  
  131.         for j in range(0,M+2):  
  132.             matrix_A_[i].append(0)  
  133.             if j<M+1:  
  134.                 matrix_A_[i][j] = array_A[i][j]  
  135.             elif j==M+1:     #如果不加这个控制条件,matrix_A_将被array_b刷新  
  136.                 matrix_A_[i][j] = array_b[i]  
  137.     return matrix_A_  
  138. matrix_A_ = matrix_element_A_(matrix_A,matrix_b,M)  
  139.   
  140.   
  141. ‘‘‘‘‘ 
  142. 多项式拟合函数 
  143. ‘‘‘  
  144. #x为自变量,w为多项式系数,m为多项式的阶数  
  145. def poly_fit(x,wp,m):  
  146.     sumf = 0  
  147.     for j in range(0,m+1):  
  148.         sumf=sumf+wp[j]*pow(x,j)  
  149.     return sumf  
  150.   
  151.   
  152. ‘‘‘‘‘ 
  153. sin(2*pi*x)在x=0处的3阶泰勒展开式 
  154. ‘‘‘  
  155. coef_taylor = [] #正弦函数的泰勒展开式系数  
  156. K=3  #展开到K阶  
  157. if K%2==0:  
  158.     print "K必须为正奇数"  
  159. s = 0  
  160. k=(K-1)/2+1  #小k为系数个数  
  161. #求K阶泰勒展开式的系数:  
  162. for i in range(0,k):  
  163.     s = pow(-1,i)*pow(2*np.pi,2*i+1)/math.factorial(2*i+1)  
  164.     coef_taylor.append(s)  
  165. print "%d阶泰勒级数展开式的系数为:" %K  
  166. print coef_taylor  
  167. #tx为泰勒展开式函数的自变量      
  168. def sin_taylor(tx):  
  169.     sum_tay=0  
  170.     for i in range(0,k):  
  171.         sum_tay=sum_tay+coef_taylor[i]*pow(tx,2*k+1)  
  172.     return sum_tay  
  173. poly_taylor_a = []   #泰勒展开式函数的输入值  
  174. poly_taylor_b = []   #泰勒展开式函数的预测值  
  175. for i in range(0,N):  
  176.     poly_taylor_a.append(a[i])  
  177.     poly_taylor_b.append(sin_taylor(poly_taylor_a[i]))  
  178.   
  179.   
  180.   
  181.   
  182. ‘‘‘‘‘ 
  183. 在x=[0,1]上生成100个点,作为测试集 
  184. ‘‘‘  
  185. testa = []  #测试集的横坐标  
  186. testb = []  #测试集的纵坐标  
  187. x=0  
  188. for i in range(0,100):  
  189.     x=x+1.0/101.0  
  190.     testa.append(x)  
  191.     testb.append(np.sin(2*np.pi*x))  
  192.       
  193. ‘‘‘‘‘ 
  194. 计算泰勒展开式模型的训练误差和测试误差 
  195. ‘‘‘  
  196. #定义误差函数:  
  197. #ly为真实值,fx为预测值  
  198. def Lfun(ly,fx):  
  199.     L=0  
  200.     for i in range(0,len(fx)):  
  201.         L=L+pow(ly[i]-fx[i],2)  
  202.     return L  
  203.   
  204.   
  205. ‘‘‘‘‘ 
  206. 主程序 
  207. ‘‘‘  
  208. if __name__ == ‘__main__‘:  
  209. # 求解方程组, 并输出方程组的可解信息  
  210.     ret = solve( matrix_A_, 0, 0 )   
  211.     if( ret== 0 ):  
  212.         print "方 程组无唯一解或无解\n"  
  213.      
  214.     # 输出方程组及其解,解即为w[j]  
  215.     w = []  
  216.     for i in range( 0, len( m ) ):  
  217.         w.append(m[i][len( m )])  
  218.     print "M=%d时的系数w[j]:" %M  
  219.     print w  
  220.       
  221.     #多项式拟合后的预测值:  
  222.     poly_a = []  
  223.     poly_b = []  
  224.     for i in range(0,N):  
  225.         poly_a.append(a[i])  
  226.         poly_b.append(poly_fit(poly_a[i],w,M))  
  227.   
  228.   
  229.     #fxtay为泰勒展开式的预测值,LCtaylor为测试误差:  
  230.     fxtay = []  
  231.     for i in range(0,100):  
  232.          fxtay.append(sin_taylor(testa[i]))  
  233.     LCtaylor = Lfun(testb,fxtay)/100  
  234.     print "三阶泰勒展开式的测试误差为:%f" %LCtaylor  
  235.   
  236.   
  237.     #fxpoly为M阶多项式拟合函数的预测值,LXpoly为训练误差:  
  238.     fxpoly = []  
  239.     for i in range(0,N):   #len(poly_b)=N=10  
  240.         fxpoly.append(poly_fit(a[i],w,M))  
  241.     LXpoly = Lfun(b,fxpoly)/len(poly_b)  
  242.     print "M=%d时多项式拟合函数的训练误差为:%f" % (M,LXpoly)  
  243.   
  244.   
  245.     #fxpolyc为M阶多项式拟合函数的预测值,LCpoly为测试误差:  
  246.     fxpolyc = []  
  247.     for i in range(0,100):  
  248.         fxpolyc.append(poly_fit(testa[i],w,M))  
  249.     LCpoly = Lfun(testb,fxpolyc)/100   
  250.     print "M=%d时多项式拟合函数的测试误差为:%f" % (M,LCpoly)  
  251.       
  252.     #多项式拟合的效果:  
  253.     fig1 = plt.figure(1)  
  254.     plt.plot(poly_a,poly_b,color=‘blue‘,linestyle=‘solid‘,marker=‘o‘)   
  255.     #加入epsilon后的样本:  
  256.     plt.plot(a,b,color=‘red‘,linestyle=‘dashed‘,marker=‘x‘)   
  257.     #泰勒展开式拟合效果:  
  258.     plt.plot(poly_taylor_a,poly_taylor_b,color=‘yellow‘,linestyle=‘dashed‘,marker=‘o‘)  
  259.     #figure(2)对比多项式拟合函数与训练数据:  
  260.     fig2 = plt.figure(2)  
  261.     plt.plot(poly_a,poly_b,color=‘blue‘,linestyle=‘solid‘,marker=‘o‘)  
  262.     plt.plot(a,b,color=‘red‘,linestyle=‘dashed‘,marker=‘x‘)  
  263.     plt.show()  


M=3时的运行结果:

 

 

[python] view plain copy
 
 技术分享技术分享
  1. 3阶泰勒级数展开式的系数为:  
  2. [6.283185307179586, -41.341702240399755]  
  3. M=3时的系数w[j]:  
  4. [-0.28492708632293295, 13.031310645420685, -37.730992850050448, 25.464782221275197]  
  5. 三阶泰勒展开式的测试误差为:100.889335  
  6. M=3时多项式拟合函数的训练误差为:0.008933  
  7. M=3时多项式拟合函数的测试误差为:0.007886  

 

 

Figure(1):

技术分享

Figure(2):

技术分享

 

         初次编写这么长的代码,思路不是有一点的混乱技术分享。其中有技术分享也有技术分享。以后会继续来优化这个程序,作为学习Python的入口。

http://blog.csdn.net/zuyuanzhu/article/details/21321007

以上是关于Python计算&绘图——曲线拟合问题(转)的主要内容,如果未能解决你的问题,请参考以下文章

Python 怎么用曲线拟合数据

最小二乘法多项式曲线拟合原理与实现(转)

添加拟合二次曲线

origin8.0中怎么在一条拟合好的曲线图上再添加一条拟合曲线?

matlab中如何拟合光滑曲线

Apache Commons Math3学习笔记 - 多项式曲线拟合(转)