標籤:
Dataset
比薩斜塔是意大利最大的旅遊景點之一。幾百年來這座塔慢慢靠向一邊,最終達到5.5度的傾斜角度,在頂端水平偏離了近3米。年度資料pisa.csv檔案記錄了從1975年到1987年測量塔的傾斜,其中lean代表了偏離的角度。在這個任務,我們將嘗試使用線性迴歸來估計傾斜率以及解釋其係數和統計資料。
# 讀取資料import pandasimport matplotlib.pyplot as pltpisa = pandas.DataFrame({"year": range(1975, 1988), "lean": [2.9642, 2.9644, 2.9656, 2.9667, 2.9673, 2.9688, 2.9696, 2.9698, 2.9713, 2.9717, 2.9725, 2.9742, 2.9757]})print(pisa)‘‘‘ lean year0 2.9642 19751 2.9644 19762 2.9656 19773 2.9667 19784 2.9673 19795 2.9688 19806 2.9696 19817 2.9698 19828 2.9713 19839 2.9717 198410 2.9725 198511 2.9742 198612 2.9757 1987‘‘‘plt.scatter(pisa["year"], pisa["lean"])
Fit The Linear Model
從散佈圖中我們可以看到用曲線可以很好的擬合該資料。在之前我們利用線性迴歸來分析葡萄酒的品質以及股票市場,但在這個任務中,我們將學習如何理解關鍵的統計學概念。Statsmodels是Python中進行嚴格統計分析的一個庫,對於線性模型,Statsmodels提供了足夠多的統計方法以及適當的評估方法。sm.OLS這個類用於擬合線性模型,採取的最佳化方法是最小二乘法。OLS()不會自動添加一個截距到模型中,但是可以自己添加一列屬性,使其值都是1即可產生截距。
import statsmodels.api as smy = pisa.lean # targetX = pisa.year # featuresX = sm.add_constant(X) # add a column of 1‘s as the constant term# OLS -- Ordinary Least Squares Fitlinear = sm.OLS(y, X)# fit modellinearfit = linear.fit()print(linearfit.summary())‘‘‘ OLS Regression Results ==============================================================================Dep. Variable: lean R-squared: 0.988Model: OLS Adj. R-squared: 0.987Method: Least Squares F-statistic: 904.1Date: Mon, 25 Apr 2016 Prob (F-statistic): 6.50e-12Time: 13:30:20 Log-Likelihood: 83.777No. Observations: 13 AIC: -163.6Df Residuals: 11 BIC: -162.4Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [95.0% Conf. Int.]------------------------------------------------------------------------------const 1.1233 0.061 18.297 0.000 0.988 1.258year 0.0009 3.1e-05 30.069 0.000 0.001 0.001==============================================================================Omnibus: 0.310 Durbin-Watson: 1.642Prob(Omnibus): 0.856 Jarque-Bera (JB): 0.450Skew: 0.094 Prob(JB): 0.799Kurtosis: 2.108 Cond. No. 1.05e+06==============================================================================Warnings:[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.[2] The condition number is large, 1.05e+06. This might indicate that there arestrong multicollinearity or other numerical problems.‘‘‘
Define A Basic Linear Model
# Our predicted values of yyhat = linearfit.predict(X)print(yhat)‘‘‘[ 2.96377802 2.96470989 2.96564176 2.96657363 2.96750549 2.96843736 2.96936923 2.9703011 2.97123297 2.97216484 2.9730967 2.97402857 2.97496044]‘‘‘residuals = yhat - y‘‘‘residuals :Series (<class ‘pandas.core.series.Series‘>)0 -0.0004221 0.0003102 0.0000423 -0.0001264 0.0002055 -0.0003636 -0.0002317 0.0005018 -0.0000679 0.00046510 0.00059711 -0.00017112 -0.000740Name: lean, dtype: float64‘‘‘
Histogram Of Residuals
- 之前我們用過長條圖(histograms )來顯示資料的分布,現在我們也可以顯示殘差的分布,來確認它是否滿足常態分佈(其實有很多統計測試來檢驗常態分佈):
plt.hist(residuals, bins=5)
- 由於我們的資料集只有13個樣本,因此這樣畫出來的長條圖並沒有太大意義,儘管中間最高的有4個樣本
Sum Of Squares
許多線性迴歸模型的統計測量都依賴於三個平方測量值:Error (SSE), Regression Sum of Squares (RSS)以及Total Sum of Squares (TSS).
-
Error (SSE):真實值與預測值的差的平方和
-
Regression Sum of Squares (RSS) :預測值和真實值的均值的差的平方和,其中的均值是真實值的均值。如果將預測值都設定為觀測值的均值,RSS會非常低,但這並沒有什麼意義。反而是一個大的RSS和一個小的SSE表示一個很好的模型。
-
Total Sum of Squares (TSS):觀測值與觀測值的均值的差的平方和,大概就是訓練集的方差。
TSS=RSS+SSE:資料總量的方差 = 模型的方差+殘差的方差
import numpy as np# sum the (predicted - observed) squaredSSE = np.sum((yhat-y.values)**2)‘‘‘1.9228571428562889e-06‘‘‘# Average yybar = np.mean(y.values)# sum the (mean - predicted) squaredRSS = np.sum((ybar-yhat)**2)‘‘‘0.00015804483516480448‘‘‘# sum the (mean - observed) squaredTSS = np.sum((ybar-y.values)**2)‘‘‘0.00015996769230769499‘‘‘print(TSS-RSS-SSE)‘‘‘3.42158959043e-17‘‘‘
R-Squared
- 線性判定(coefficient of determination)也叫R-Squared,是用來測定線性依賴性的。它是一個數字用來告訴我們資料的總方差中模型的方差的佔比:
- 前面提到一個低的SSE和一個高的RSS表示一個很好的模型擬合,這個R-Squared就表示了這個意思,介於0到1之間。
R2 = RSS/TSSprint(R2)‘‘‘0.987979715684‘‘‘
T-Distribution
統計測驗表明塔的傾斜程度與年份有關係,一個常見的統計顯著性測試是student t-test。這個測試的基礎是T分布。和常態分佈很相似,都是鐘型但是峰值較低。T檢驗是用於小樣本(樣本容量小於30)的兩個平均值差異程度的檢驗方法。它是用T分布理論來推斷差異發生的機率,從而判定兩個平均數的差異是否顯著。
from scipy.stats import t# 100 values between -3 and 3x = np.linspace(-3,3,100)# Compute the pdf with 3 degrees of freedomprint(t.pdf(x=x, df=3))‘‘‘[ 0.02297204 0.02441481 0.02596406 0.02762847 0.0294174 0.031341 0.03341025 0.03563701 0.03803403 0.04061509 0.04339497 0.04638952 0.04961567 0.05309149 0.05683617 0.06086996 0.0652142 0.06989116 0.07492395 0.08033633 0.08615245 0.09239652 0.0990924 0.10626304 0.11392986 0.12211193 0.13082504 0.14008063 0.14988449 0.16023537 0.17112343 0.18252859 0.1944188 0.20674834 0.21945618 0.23246464 0.2456783 0.2589835 0.27224841 0.28532401 0.29804594 0.31023748 0.32171351 0.33228555 0.34176766 0.34998293 0.35677032 0.36199128 0.36553585 0.36732769 0.36732769 0.36553585 0.36199128 0.35677032 0.34998293 0.34176766 0.33228555 0.32171351 0.31023748 0.29804594 0.28532401 0.27224841 0.2589835 0.2456783 0.23246464 0.21945618 0.20674834 0.1944188 0.18252859 0.17112343 0.16023537 0.14988449 0.14008063 0.13082504 0.12211193 0.11392986 0.10626304 0.0990924 0.09239652 0.08615245 0.08033633 0.07492395 0.06989116 0.0652142 0.06086996 0.05683617 0.05309149 0.04961567 0.04638952 0.04339497 0.04061509 0.03803403 0.03563701 0.03341025 0.031341 0.0294174 0.02762847 0.02596406 0.02441481 0.02297204]‘‘‘
比薩鐵塔——統計顯著性檢驗