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适用于多变量回归，它可以帮助我们剔除掉一些没有相关性的属性。