Python計算&繪圖——曲線擬合問題（轉）

題目來自老師的課后作業，如下所示。很多地方應該可以直接調用函數，但是初學Python，對里面的函數還不是很了解，順便帶着學習的態度，盡量自己動手code。

         

 

測試版代碼，里面帶有很多注釋和測試代碼：

 

# -*- 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][n], m[ik][n] );
        try:
            # 消元
            m[row_pos][n] /= 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][n] -= m[row_pos][n] * 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][n] ) == 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):

 

         初次編寫這么長的代碼，思路不是有一點的混亂。其中有也有。以后會繼續來優化這個程序，作為學習Python的入口。

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