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

         

Python计算&绘图——曲线拟合问题(转)_线性方程组

 

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

 


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. 0
16. 0.1
17. #高斯分布随机数
18. return np.sin(2*np.pi*x)+epsilon
19. for i in range(0,10):
20. 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. 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. 0; j = 0; row_pos = 0; col_pos = 0; ik = 0; jk = 0
49. 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. 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. 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. 1
74. while( j >= col_pos ):
75. m[row_pos][j] /= m[row_pos][col_pos]
76. 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. 1
82. while( j >= col_pos ):
83. m[i][j] -= m[row_pos][j] * m[i][col_pos]
84. 1
85. 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. 0
105. for k in range(0,n):
106. 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. 0)
114. 1,j+1,N)
115. #计算线性方程组矩阵b的第[i]行元素b[i]
116. def matrix_element_b(x,y,i,n):
117. 0
118. for k in range(0,n):
119. 1) #x[0]到x[n-1],共n个元素求和
120. return sum_b
121. for i in range(0,M+1):
122. 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,与全局变量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. 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. 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. 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. 0
170. for i in range(0,k):
171. 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. 1.0/101.0
190. testa.append(x)
191. 2*np.pi*x))
192.
193. '''''
194. 计算泰勒展开式模型的训练误差和测试误差
195. '''
196. #定义误差函数:
197. #ly为真实值,fx为预测值
198. def Lfun(ly,fx):
199. 0
200. for i in range(0,len(fx)):
201. 2)
202. return L
203.
204.
205. '''''
206. 主程序
207. '''
208. if __name__ == '__main__':
209. # 求解方程组, 并输出方程组的可解信息
210. 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. 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. 100
250. print "M=%d时多项式拟合函数的测试误差为:%f" % (M,LCpoly)
251.
252. #多项式拟合的效果:
253. 1)
254. 'blue',linestyle='solid',marker='o')
255. #加入epsilon后的样本:
256. 'red',linestyle='dashed',marker='x')
257. #泰勒展开式拟合效果:
258. 'yellow',linestyle='dashed',marker='o')
259. #figure(2)对比多项式拟合函数与训练数据:
260. 2)
261. 'blue',linestyle='solid',marker='o')
262. 'red',linestyle='dashed',marker='x')
263. plt.show()

M=3时的运行结果:

 

 



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):

Python计算&绘图——曲线拟合问题(转)_方程组_02

Figure(2):

Python计算&绘图——曲线拟合问题(转)_线性方程组_03

 

         初次编写这么长的代码,思路不是有一点的混乱

。其中有

也有

。以后会继续来优化这个程序,作为学习Python的入口。