python回歸分析

假設原函數由一個三角函數和一個線性項組成

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def f(x):
    return np.sin(x) + 0.5 * x

x = np.linspace(-2 * np.pi, 2 * np.pi, 50)

plt.plot(x, f(x), 'b')
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')

一、用回歸方式逼近

1. 作為基函數的單項式

最簡單的情況是以單項式為基函數——也就是說，b₁=1,b₂=x,b₃=x²,b₄=x³,... 在這種情況下，Numpy有確定最優參數（polyfit）和以一組輸入值求取近似值（ployval）的內建函數。ployfit函數參數如下：

參數	描述
x	x坐標（自變量值）
y	y坐標（因變量值）
deg	多項式擬合度
full	如果有真，返回額外的診斷信息
w	應用到y坐標的權重
cov	如果為真，返回協方差矩陣

典型向量化風格的polyfit和polyval線性回歸（deg=7）應用方式如下：

reg = np.polyfit(x, f(x), deg=7)
ry = np.polyval(reg, x)
plt.plot(x,f(x),'b',label='f(x)')
plt.plot(x,ry,'r.',label='regression')
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')

 

2. 單獨的基函數

當選擇更好的基函數組時，可以得到更好的回歸結果。單獨的基函數必須能通過一個矩陣方法定義（使用Numpy ndarray對象）。

numpy.linalg子庫提供lstsq函數，以解決最小二乘化做問題

matrix = np.zeros((3+1,len(x)))
matrix[3,:] = np.sin(x)
matrix[2,:] = x **2
matrix[1,:] = x
matrix[0,:] = 1

reg = np.linalg.lstsq(matrix.T,f(x))[0]
ry = np.dot(reg,matrix)

plt.plot(x,f(x),'b',label='f(x)')
plt.plot(x,ry,'r.',label='regression')
plt.legend(loc = 0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')

 

 對於有噪聲的數據，未排序的數據，回歸法都可處理。

3. 多維

以fm函數為例

def fm(x,y):
    return np.sin(x) + 0.25 *x +np.sqrt(y) + 0.05 * y **2

x = np.linspace(0,10,20)
y = np.linspace(0,10,20)
X, Y = np.meshgrid(x,y)
Z = fm(X,Y)
x = X.flatten()
y = Y.flatten()

from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl

fig = plt.figure(figsize=(9, 6))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X,Y,Z,rstride=2,cstride=2,cmap=mpl.cm.coolwarm,linewidth=0.5,antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x,y)')
fig.colorbar(surf,shrink=0.5,aspect=5)

為了獲得好的回歸結果，編譯一組基函數，包括一個sin函數和sqrt函數。statsmodels庫提供相當通過和有益的函數OLS，可以用於一維和多維最小二乘回歸。

matrix = np.zeros((len(x),6+1))
matrix[:,6] = np.sqrt(y)
matrix[:,5] = np.sin(x)
matrix[:,4] = y ** 2
matrix[:,3] = x **2 
matrix[:,2] = y
matrix[:,1] = x
matrix[:,0] = 1

import statsmodels.api as sm 
model = sm.OLS(fm(x,y),matrix).fit()

OSL函數的好處之一是提供關於回歸及其而非結果的大量幾百萬來看信息。調用model.summary可以訪問結果的一個摘要。單獨統計數字（如確定系數）通常可以直接訪問model.rsquared。最優回歸系數，保存在model對象的params屬性中。

model.rsquared
a = model.params

def reg_func(a,x,y):
    f6 = a[6] * np.sqrt(y)
    f5 = a[5] * np.sin(x)
    f4 = a[4] * y ** 2
    f3 = a[3] * x ** 2
    f2 = a[2] * y
    f1 = a[1] * x
    f0 = a[0] * 1
    return (f6+f5+f4+f3+f2+f1+f0)
# reg_func返回給定最優回歸參數和自變量數據點的回歸函數值
RZ = reg_func(a,X,Y)
fig = plt.figure(figsize=(9, 6))
ax = fig.gca(projection='3d')
surf1 = ax.plot_surface(X,Y,Z,rstride=2,cstride=2,cmap=mpl.cm.coolwarm,linewidth=0.5,antialiased=True)
surf2 = ax.plot_wireframe(X,Y,RZ,rstride=2,cstride=2,label = 'regression')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x,y)')
fig.colorbar(surf,shrink=0.5,aspect=5)