Dataset
比薩斜塔是義大利最大的旅遊景點之一。幾百年來這座塔慢慢靠向一邊,最終達到5.5度的傾斜角度，在頂端水平偏離了近3米。年度數據pisa.csv文件記錄了從1975年到1987年測量塔的傾斜，其中lean代表了偏離的角度。在這個任務,我們將嘗試使用線性回歸來估計傾斜率以及解釋其係數和統計數據。
# 讀取數據
import pandas
import matplotlib.pyplot as plt
pisa = 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 year
0 2.9642 1975
1 2.9644 1976
2 2.9656 1977
3 2.9667 1978
4 2.9673 1979
5 2.9688 1980
6 2.9696 1981
7 2.9698 1982
8 2.9713 1983
9 2.9717 1984
10 2.9725 1985
11 2.9742 1986
12 2.9757 1987
'''
plt.scatter(pisa["year"], pisa["lean"])
Fit The Linear Model
從散點圖中我們可以看到用曲線可以很好的擬合該數據。在之前我們利用線性回歸來分析葡萄酒的質量以及股票市場，但在這個任務中，我們將學習如何理解關鍵的統計學概念。Statsmodels是Python中進行嚴格統計分析的一個庫，對於線性模型，Statsmodels提供了足夠多的統計方法以及適當的評估方法。sm.OLS這個類用於擬合線性模型，採取的優化方法是最小二乘法。OLS()不會自動添加一個截距到模型中，但是可以自己添加一列屬性，使其值都是1即可產生截距。
import statsmodels.api as sm
y = pisa.lean # target
X = pisa.year # features
# add a column of 1's as the constant term
X = sm.add_constant(X)
# OLS -- Ordinary Least Squares Fit
linear = sm.OLS(y, X)
# fit model
linearfit = linear.fit()
print(linearfit.summary())
'''
OLS Regression Results
==============================================================================
Dep. Variable: lean R-squared: 0.988
Model: OLS Adj. R-squared: 0.987
Method: Least Squares F-statistic: 904.1
Date: Mon, 25 Apr 2016 Prob (F-statistic): 6.50e-12
Time: 13:30:20 Log-Likelihood: 83.777
No. Observations: 13 AIC: -163.6
Df Residuals: 11 BIC: -162.4
Df 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.258
year 0.0009 3.1e-05 30.069 0.000 0.001 0.001
==============================================================================
Omnibus: 0.310 Durbin-Watson: 1.642
Prob(Omnibus): 0.856 Jarque-Bera (JB): 0.450
Skew: 0.094 Prob(JB): 0.799
Kurtosis: 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 are
strong multicollinearity or other numerical problems.
'''
Define A Basic Linear Model
列印summary時發現有很多關於模型的信息，為了弄清楚這些統計指標，我們需要仔細研究一下標準的線性回歸模型，下面模型中的ei是預測值和真實值的差，我們默認模型的誤差是正態分布的，均值為0.：
計算模型的殘差(residuals)：
# Our predicted values of y
yhat = 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 ()
0 -0.000422
1 0.000310
2 0.000042
3 -0.000126
4 0.000205
5 -0.000363
6 -0.000231
7 0.000501
8 -0.000067
9 0.000465
10 0.000597
11 -0.000171
12 -0.000740
Name: lean, dtype: float64
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
Histogram Of Residuals
之前我們用過直方圖(histograms )來顯示數據的分布，現在我們也可以顯示殘差的分布，來確認它是否滿足正態分布(其實有很多統計測試來檢驗正態分布)：
plt.hist(residuals, bins=5)
1
由於我們的數據集只有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) squared
SSE = np.sum((yhat-y.values)**2)
'''
1.9228571428562889e-06
'''
# Average y
ybar = np.mean(y.values)
# sum the (mean - predicted) squared
RSS = np.sum((ybar-yhat)**2)
'''
0.00015804483516480448
'''
# sum the (mean - observed) squared
TSS = 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越大越好，1最好。
R2 = RSS/TSS
print(R2)
'''
0.987979715684
'''
T-Distribution
統計測驗表明塔的傾斜程度與年份有關係，一個常見的統計顯著性測試是student t-test。這個測試的基礎是T分布，和正態分布很相似，都是鍾型但是峰值較低。顯著性檢驗利用了T分布的概率密度函數*probability density function (pdf)，pdf描述了兩個連續隨機變量的相關性。這個隨機變量的累計密度函數cumulative density function (cdf)小於或等於某一個點。自由度degrees of freedom (df) 描述的是觀測樣本的數量(一般將其設置為樣本數-2)。T檢驗是用於小樣本(樣本容量小於30)的兩個平均值差異程度*的檢驗方法。它是用T分布理論來推斷差異發生的概率，從而判定兩個均值的差異是否顯著。
from scipy.stats import t
# 100 values between -3 and 3
x = np.linspace(-3,3,100)
# Compute the pdf with 3 degrees of freedom
print(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]
'''
Statistical Significance Of Coefficients
統計顯著性首先：我們要測試這個斜塔的傾斜程度是否與年份有關，我們將null hypothesis(H0)：沒有關係，也就是年份,這個屬性的係數(coefficient )β1=0，相反β1≠0。
1.提出一個假設：H0:β1=0 H1:β1≠0
然後測試null hypothesis，我們需要利用T分布計算一個統計測量值：
2.計算t-statistic的計算公式如下：
tstat = linearfit.params["year"] / np.sqrt(s2b1)
'''
30.068584687652137
'''
現在我們得到了t值，我們需要計算這個β1與0不同的概率，我們設置顯著性為95%，意思就是β1與0不同的概率大於95%，我們才認為β1≠0。這需要用到t分布，給定一個p值和自由度計算這個t分布的累積密度函數，我們就可以獲得一個概率：
在雙邊測試中，我們還需要觀察這個相關性是正相關還是負相關。T分布式關於原點對稱的，由於我們此處是看是否相關，因此兩邊的和是5%，此處的這個顯著性值要>97.5%，因為是雙邊的。
If Tcdf(|t|,df) < 0.975 ： Accept H0:β1=0
Else ：Accept H1
# At the 95% confidence interval for a two-sided t-test we must use a p-value of 0.975
pval = 0.975
# The degrees of freedom
df = pisa.shape[0] - 2
# The probability to test against
p = t.cdf(tstat, df=df)
'''
0.99999999999674838
'''
beta1_test = p > pval
'''
True
'''
最終p值>0.975因此接受β1≠0這個假設，也就是斜塔的傾斜度與年份顯著性相關。
Conclusion
此文中給出了如何計算截距的方法：添加一列值為1的屬性
R-squared是一個很強大的度量值，但是它經常被過度使用，一個很低的R-squared不代表這兩個變量之間沒有依賴性，比如y=sin(x)的R-squared的值為0，但是很明顯x和y是相關的。另一方面，一個很高的R-squared值也不代表這個模型能很好的預測未來的事件，因為它沒有計算觀測樣本的數量。
Student t-tests適用於多變量回歸，它可以幫助我們剔除掉一些沒有相關性的屬性。