整個過程分七步，為了方便喜歡直接copy代碼看結果的同學，每步都放上了完整的代碼。
實驗數據：
第一步：准備樣本數據並繪制散點圖
1）代碼及其說明
import numpy as np import scipy as sp import matplotlib.pyplot as plt from scipy.optimize import leastsq ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) #身高 Yi=np.array([58,63,57,65,62,66,58,59,62])#體重 #畫樣本點 plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 plt.scatter(Xi,Yi,color="green",label="樣本數據",linewidth=1) plt.show()
2）結果圖
3）分析
從散點圖可以看出，樣本點基本是圍繞箭頭所示的直線分布的。所以先以直線模型對數據進行擬合
第二步： 使用最小二乘法算法求擬合直線
1）代碼及其說明
import numpy as np import scipy as sp import matplotlib.pyplot as plt from scipy.optimize import leastsq ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) ##需要擬合的函數func :指定函數的形狀 k= 0.42116973935 b= -8.28830260655 def func(p,x): k,b=p return k*x+b ##偏差函數：x,y都是列表:這里的x,y更上面的Xi,Yi中是一一對應的 def error(p,x,y): return func(p,x)-y #k,b的初始值，可以任意設定,經過幾次試驗，發現p0的值會影響cost的值：Para[1] p0=[1,20] #把error函數中除了p0以外的參數打包到args中(使用要求) Para=leastsq(error,p0,args=(Xi,Yi)) #讀取結果 k,b=Para[0] print("k=",k,"b=",b) #畫樣本點 plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 plt.scatter(Xi,Yi,color="green",label="樣本數據",linewidth=2) #畫擬合直線 x=np.linspace(150,190,100) ##在150-190直接畫100個連續點 y=k*x+b ##函數式 plt.plot(x,y,color="red",label="擬合直線",linewidth=2) plt.legend() #繪制圖例 plt.show()
2）結果圖
3）分析
從圖上看，擬合效果還是不錯的。樣本點基本均勻的分布在回歸線兩邊，沒有出現數據點嚴重偏離回歸線的情況。
第三步： 驗證回歸線的擬合程度—殘差分布圖
1）代碼及其說明
import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.graphics.api import qqplot ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) xy_res=[] ##計算殘差 def residual(x,y): res=y-(0.42116973935*x-8.28830260655) return res ##讀取殘差 for d in range(0,len(Xi)): res=residual(Xi[d],Yi[d]) xy_res.append(res) ##print(xy_res) ##計算殘差平方和:22.8833439288 -->越小擬合情況越好 xy_res_sum=np.dot(xy_res,xy_res) #print(xy_res_sum) ##如果數據擬合模型效果好，殘差應該遵從正態分布(0,d*d:這里d表示殘差) #畫樣本點 fig=plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 ax=fig.add_subplot(111) fig=qqplot(np.array(xy_res),line='q',ax=ax) plt.show()
2）結果圖
3）分析
上圖為Q-Q圖，原理：如果兩個分布相似，則該Q-Q圖趨近於落在y=x線上。如果兩分布線性相關，則點在Q-Q圖上趨近於落在一條直線上，但不一定在y=x線上。Q-Q圖可以用來可在分布的位置-尺度范疇上可視化的評估參數。
從圖上可以看出，回歸效果比較理想，但不是最理想的
4）以下代碼可以同樣實現上述圖示：
import numpy as np import scipy.stats as stats import pylab ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) xy_res=[] ##計算殘差 def residual(x,y): res=y-(0.42116973935*x-8.28830260655) return res ##讀取殘差 for d in range(0,len(Xi)): res=residual(Xi[d],Yi[d]) xy_res.append(res) ##print(xy_res) ##計算殘差平方和:22.8833439288 -->越小擬合情況越好 xy_res_sum=np.dot(xy_res,xy_res) #print(xy_res_sum) ##如果數據擬合模型效果好，殘差應該遵從正態分布(0,d*d:這里d表示殘差) #畫樣本點 stats.probplot(np.array(xy_res),dist="norm",plot=pylab) pylab.show()
第四步： 驗證回歸線的擬合程度—標准化殘差
1）代碼及其說明
import numpy as np import matplotlib.pyplot as plt ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) xy_res=[] ##計算殘差 def residual(x,y): res=y-(0.42116973935*x-8.28830260655) return res ##讀取殘差 for d in range(0,len(Xi)): res=residual(Xi[d],Yi[d]) xy_res.append(res) ##print(xy_res) ##計算殘差平方和:22.8833439288 -->越小擬合情況越好 xy_res_sum=np.dot(xy_res,xy_res) ''' 標准殘差: （殘差-殘差平均值）/殘差的標准差 ''' ''' 標准殘差圖： 1.標准殘差是以擬合模型的自變量為橫坐標，以標准殘差為縱坐標形成的平面坐標圖像 2.試驗點的標准殘差落在殘差圖的(-2,2)區間以外的概率<=0.05.若某一點落在區間外，可判為異常點 3.有效標准殘差點圍繞y=0的直線上下完全隨機分布，說明擬合情況良好 4.如果擬合方程原本是非線性模型，但擬合時卻采用了線性模型，標准化殘差圖就會表現出曲線形狀，產生 系統性偏差 ''' ##計算殘差平均值 xy_res_avg=0 for d in range(0,len(xy_res)): xy_res_avg+=xy_res[d] xy_res_avg/=len(xy_res) #殘差的標准差 xy_res_sd=np.sqrt(xy_res_sum/len(Xi)) ##標准化殘差 xy_res_sds=[] for d in range(0,len(Xi)): res=(xy_res[d]-xy_res_avg)/xy_res_sd xy_res_sds.append(res) #print(xy_res_sds) #標准化殘差分布 plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 plt.scatter(Xi,xy_res_sds) plt.show() ''' 1.絕大部分數據分布在(-2,+2)的水平帶狀區間內，因此模型擬合較充分 2.數據點分布稍均勻，但沒有達到隨機均勻分布的狀態。此外，部分數據點還呈現某種曲線波動形狀， 有少許系統性偏差。因此可能采用非線性擬合效果會更好 '''
2）結果圖
3）分析
數據點分布還是存在一定的變化趨勢的。
第五步：使用曲線模型擬合數據
1）代碼及其說明
import numpy as np import scipy as sp import matplotlib.pyplot as plt from scipy.optimize import leastsq ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) ##需要擬合的函數func :指定函數的形狀 k= 0.860357336936 b= -19.6659389666 def func(p,x): k,b=p return x**k+b ##偏差函數：x,y都是列表:這里的x,y更上面的Xi,Yi中是一一對應的 def error(p,x,y): return func(p,x)-y #k,b的初始值，可以任意設定,經過幾次試驗，發現p0的值會影響cost的值：Para[1] p0=[1,20] #把error函數中除了p0以外的參數打包到args中(使用要求) Para=leastsq(error,p0,args=(Xi,Yi)) #讀取結果 k,b=Para[0] print("k=",k,"b=",b) #畫樣本點 plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 plt.scatter(Xi,Yi,color="green",label="樣本數據",linewidth=2) #畫擬合直線 x=np.linspace(150,190,100) ##在150-190直接畫100個連續點 y=x**k+b ##函數式 plt.plot(x,y,color="red",label="擬合直線",linewidth=2) plt.legend() #繪制圖例 plt.show()
2）結果圖
3）分析
由於標准化殘差的分布圖，部分數據的趨勢與冪函數在第一象限的圖像類似， 所以采用了y=xa +b的函數形式，截距b是為了圖像可以在Y軸上下移動
第六步：驗證回歸線的擬合程度—殘差分布圖
1）代碼及其說明
import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.graphics.api import qqplot ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) xy_res=[] ##計算殘差 def residual(x,y): res=y-(x**0.860357336936-19.6659389666) return res ##讀取殘差 for d in range(0,len(Xi)): res=residual(Xi[d],Yi[d]) xy_res.append(res) ##print(xy_res) ##計算殘差平方和:22.8833439288 -->越小擬合情況越好 xy_res_sum=np.dot(xy_res,xy_res) #print(xy_res_sum) ##如果數據擬合模型效果好，殘差應該遵從正態分布(0,d*d:這里d表示殘差) #畫樣本點 fig=plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 ax=fig.add_subplot(111) fig=qqplot(np.array(xy_res),line='q',ax=ax) plt.show()
2）結果圖
3）分析
從圖上可以看出，回歸效果也比較理想
第七步：驗證回歸線的擬合程度—標准化殘差
1）代碼及其說明
import numpy as np import matplotlib.pyplot as plt ##樣本數據(Xi,Yi)，需要轉換成數組(列表)形式 Xi=np.array([160,165,158,172,159,176,160,162,171]) Yi=np.array([58,63,57,65,62,66,58,59,62]) xy_res=[] ##計算殘差 def residual(x,y): res=y-(x**0.860357336936-19.6659389666) return res ##讀取殘差 for d in range(0,len(Xi)): res=residual(Xi[d],Yi[d]) xy_res.append(res) ##print(xy_res) ##計算殘差平方和:22.881076636 -->越小擬合情況越好 xy_res_sum=np.dot(xy_res,xy_res) #print(xy_res_sum) ''' 標准殘差: （殘差-殘差平均值）/殘差的標准差 ''' ##計算殘差平均值 xy_res_avg=0 for d in range(0,len(xy_res)): xy_res_avg+=xy_res[d] xy_res_avg/=len(xy_res) #殘差的標准差 xy_res_sd=np.sqrt(xy_res_sum/len(Xi)) ##標准化殘差 xy_res_sds=[] for d in range(0,len(Xi)): res=(xy_res[d]-xy_res_avg)/xy_res_sd xy_res_sds.append(res) print(xy_res_sds) #標准化殘差分布 plt.figure(figsize=(8,6)) ##指定圖像比例： 8：6 plt.scatter(Xi,xy_res_sds) plt.show() ''' 1.絕大部分數據分布在(-2,+2)的水平帶狀區間內，因此模型擬合較充分 2.數據點分布稍均勻，但沒有達到隨機均勻分布的狀態。此外，部分數據點還呈現某種曲線波動形狀， 有少許系統性偏差。因此可能采用非線性擬合效果會更好 '''
2）結果圖
3）分析
數據點分布趨和直線回歸方程基本一樣
補充說明：
其實整個實驗過程並沒有達到預期效果。
1）如果對實驗過程的5-7步使用R語言重新實驗(R語言提供了所有相關函數)，第7步的效果如下：
也就說所有的標准化殘差都是落在(-2,+2)區間內的，即曲線方程才是最佳擬合方程。
2）標准化殘差沒有找到具體的定義，網上對這個定義有多種解釋
3）標准化殘差的計算方式沒有找到相應的python包，只能按照其中某一個定義自己寫代碼計算(估計是浮點數計算產生的誤差)