Python知識(7)--最小二乘求解

這里展示利用python實現的最小二乘的直接求解方法。其求解原理，請參考:最小二乘法擬合非線性函數及其Matlab/Excel 實現

1.一般曲線擬合

代碼如下：

# -*- coding:utf-8 -*-  


'''
Created on Feb 20, 2017
最小二乘擬合 給定的函數 fit_fun(x)
已知數據X（2xN）,Y(1xN),可直接計算直線的參數a和b
直接計算的公式：ab = inv(XXt)*X*Yt
@author: Xiankai Chen
@email: [e-mail]
'''
import numpy as np
import matplotlib.pyplot as plt

def fun2ploy(x,n):
    '''
    數據轉化為[x^0,x^1,x^2,...x^n]
    '''
    lens = len(x)
    X = np.ones([1,lens]);
    for i in range(1,n):
        X = np.vstack((X,np.power(x,i)))
    return X  

def leastseq_byploy(x,y,ploy_dim):
    '''
    最小二乘求解
    '''
    #散點圖
    plt.scatter(x,y,color="r",marker='o',s = 50)

    X = fun2ploy(x,ploy_dim);
    #直接求解
    Xt = X.transpose();
    XXt=X.dot(Xt);
    XXtInv = np.linalg.inv(XXt)
    XXtInvX = XXtInv.dot(X)
    coef = XXtInvX.dot(y.T)
    
    y_est = Xt.dot(coef)
    
    return y_est,coef

    
def fit_fun(x):
    '''
    要擬合的函數
    '''
    #return np.power(x,5)
    return np.sin(x) 
    #return 5*x+3
        
if __name__ == '__main__':
    data_num = 100;
    ploy_dim =10;
    noise_scale = 0.2;
    ## 數據准備
    x = np.array(np.linspace(-2*np.pi,2*np.pi,data_num))    
    y = fit_fun(x)+noise_scale*np.random.rand(1,data_num)
    # 最小二乘擬合
    [y_est,coef] = leastseq_byploy(x,y,10)
    
    #顯示擬合結果
    org_data = plt.scatter(x,y,color="r",marker='o',s = 50)
    est_data = plt.plot(x,y_est,color="b",linewidth= 3)
    
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.title("Fit funtion with leastseq method")
    plt.legend(["Noise data","Fited function"]);
    plt.show()
    

結果圖示

當維度增加求解矩陣的你運算會消耗較大的計算資源，通常采用梯度下降法，牛頓法等數值迭代算法進行求解。

參考資料：

[1]. 最小二乘法擬合非線性函數及其Matlab/Excel 實現