一元回归分析

用户3577892

发布于 2021-01-14 10:29:41

1.1K0

发布于 2021-01-14 10:29:41

文章被收录于专栏：数据科学CLUB

理论部分
给出样本数据
计算相应指标
可视化

理论部分

问题

\quad

考察两个变量

与

之间是否存在线性相关关系,其中

是一般 ( 可控) 变量,

是随机变量,其线性相关关系可表示如下 ( 可用散点图显示) :

y=\beta_{0}+\beta_{1} x+\varepsilon

其中

\beta_{0}

为截距,

\beta,

为斜率

, \varepsilon

为随机误差,常假设

\varepsilon \sim N\left(0, \sigma^{2}\right) .

这里

\beta_{0}, \beta_{1}, \sigma^{2}

是三个待估参数. 上式表明,

与

之间有线性关系,但受到随机误差的干扰.

数据

\quad

对

与

通过试验或观察可得

对数据（注 : 数据是成对的,不允许错位). 在

与

之间存在线性关系的假设下,有如下统计模型：

\left\{\begin{array}{l} y_{i}=\beta_{0}+\beta_{1} x_{i}+\varepsilon_{i}, i=1,2, \cdots, n \\ \text { 各 } \varepsilon_{i} \text { 独立同分布,其分布为 } N\left(0, \sigma^{2}\right) . \end{array}\right.

利用成对数据可获得

\boldsymbol{\beta}_{0}

与

\boldsymbol{\beta}_{1}

的估计,设估计分别为

\hat{\boldsymbol{\beta}}_{0}

与

\hat{\boldsymbol{\beta}}_{1},

则称

\hat{y}=\hat{\boldsymbol{\beta}}_{0}+\hat{\boldsymbol{\beta}}_{1} x

为回归方程,其图形称为回归直线.

参数估计

\quad

用最小二乘法可得

\beta_{0}

与

\beta_{1}

的无偏估计

\left\{\begin{array}{l} \hat{\beta}_{1}=l_{x y} / l_{x x} \\ \hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1} \bar{x} \end{array}\right.

其中

\bar{x}=\frac{1}{n} \sum x_{i}, \bar{y}=\frac{1}{n} \sum y_{i}\left(\right.

此处

\sum

表示

\sum_{i=1}^{n},

下同

\begin{array}{l} l_{x y}=\sum\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)=\sum x_{i} y_{i}-n \bar{x} \cdot \bar{y}=\sum x_{i} y_{i}-\frac{1}{n} \sum x_{i} \sum y_{i} \\ l_{s x}=\sum\left(x_{i}-\bar{x}\right)^{2}=\sum x_{i}^{2}-n \bar{x}^{2}=\sum x_{i}^{2}-\frac{1}{n}\left(\sum x_{i}\right)^{2} \\ l_{y}=\sum\left(y_{i}-\bar{y}\right)^{2}=\sum y_{i}^{2}-n \bar{y}^{2}=\sum y_{i}^{2}-\frac{1}{n}\left(\sum y_{i}\right)^{2} \end{array}

回归方程的显著性检验

\quad

回占方程的显著性检验就是要对如下一对假设作出判断：

\boldsymbol{H}_{0}: \boldsymbol{\beta}_{1}=0 \quad \text { vs } \quad H_{1}: \boldsymbol{\beta}_{1} \neq 0

检验方法如下.

检验如下的平方和分解式是非常重要的,它在许多统计领域得到应用 :

S_{r}=S_{R}+S_{e}, f_{T}=f_{R}+f_{e}

其中

S_{T}=\sum\left(y_{i}-\bar{y}\right)^{2}=\sum y_{i}^{2}-\frac{1}{n}\left(\sum y_{i}\right)^{2}=l_{y y}

是总平方和

其自由度

f_{T}=n-1

S_{R}=\sum\left(\hat{y}_{i}-\bar{y}\right)^{2}=\sum\left(\hat{\beta}_{0}+\hat{\beta}_{1} x_{i}-\bar{y}\right)^{2}=\hat{\beta}_{1} l_{x y}=\hat{\beta}_{1}^{2} l_{x x}

是回归平方和,其自由度

f_{R}=1

S_{e}=\sum\left(y_{i}-\hat{y}_{i}\right)^{2}=\sum\left(y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1} x_{i}\right)^{2}

是残差平方和,其自由度

f_{e}=n-2 .

而

\hat{y}_{i}=\hat{\beta}_{0}+\hat{\beta}_{1} x_{i}

是在

x=x_{i}

的回归值(拟合值),它与实测值

y_{i}

通常是不相等的. 在原假设

H_{0}

成立的条件下,检验统计胃

F=\frac{S_{R}}{S_{e} /(n-2)} \sim F(1, n-2)

,拒绝域为

W=\left\{F \geqslant F_{1-\alpha}(1, n-2)\right\}

上述检验过程一般用如下方差分析表列出 :

估计与预测

\bullet

当

x=x_{0}

时

, \hat{y}_{0}=\hat{\beta}_{0}+\hat{\beta}_{1} x_{0}

是

E\left(y_{0}\right)=\beta_{0}+\beta_{1} x_{0}

的点估计

;

\bullet

当

x=x_{0}

时

, E\left(y_{0}\right)=\boldsymbol{\beta}_{0}+\boldsymbol{\beta}_{1} x_{0}

的置信水平由

1-\alpha

的置信区间是

\left[\hat{y}_{0}-\delta_{0},\right.

\left.\hat{y}_{0}+\delta_{0}\right],

其中

\delta_{0}=t_{1-\alpha / 2}(n-2) \hat{\sigma} \sqrt{\frac{1}{n}+\frac{\left(x_{0}-\bar{x}\right)^{2}}{l_{x x}}}, \quad \hat{\sigma}=\sqrt{M S_{e}}

\bullet

当

x=x_{0}

时

, y_{0}=\beta_{0}+\beta_{1} x_{0}+\varepsilon_{0}

的

1-\alpha

预测区间是

\left[\hat{y}_{0}-\delta, \hat{y}_{0}+\delta\right]

,其中

\delta=\delta\left(x_{0}\right)=t_{1-\alpha / 2}(n-2) \hat{\sigma} \sqrt{1+\frac{1}{n}+\frac{\left(x_{0}-\bar{x}\right)^{2}}{l_{x x}}}

注

: E\left(y_{0}\right)

是未知参数,而

y_{0}

是随机变量. 对

E\left(y_{0}\right)

谈论的是置信区间,对

y_{0}

谈论的是预测区间,两者是不同的,显然,预测区间要比置信区间宽很多. 要提高预测区间(置信区间也一样) 的精度,即要使

\delta\left(\right.

或

\left.\delta_{0}\right)

较小,这要求 : (1) 增大样本量

n ;(2)

增大

l_{x x},

即要求

x_{1}, x_{2}, \cdots, x_{n}

较为分散

;(3)

使

x_{0}

靠近

\bar{x} .

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

给出样本数据

x = np.array([15.3, 10.8, 8.1, 19.5, 7.2, 5.3, 9.3, 11.1, 7.5, 12.2,
                6.7, 5.2, 19.0, 15.1, 6.7, 8.6, 4.2, 10.3, 12.5, 16.1, 
                13.3, 4.9, 8.8, 9.5])
y = np.array([1.76, 1.34, 1.27, 1.47, 1.27, 1.49, 1.31, 1.09, 1.18, 
                1.22, 1.25, 1.19, 1.95, 1.28, 1.52, 1.5, 1.12, 1.37, 
                1.19, 1.05, 1.32, 1.03, 1.12, 1.70])

计算相应指标

n = len(x)   
Lxx = np.sum(x**2) - np.sum(x)**2/n
Lyy = np.sum(y**2) - np.sum(y)**2/n    
Lxy = np.sum(x*y) - np.sum(x)*np.sum(y)/n    
mean_x = np.mean(x)
mean_y = np.mean(y)

# 斜率和截距的最小二乘估计和MLE是一样的
b = Lxy/Lxx
a = mean_y - b*mean_x
fit = lambda xx: a + b*xx  

alpha = 0.05
# 残差
residuals = y - fit(x)
# MSE
var_res = np.sum(residuals**2)/(n-2)
sd_res = np.sqrt(var_res)

# 残差自由度
df = n-2     
# t值
tval = stats.t.isf(alpha/2., df)  

# 置信区间
se_fit     = lambda x: sd_res * np.sqrt(  1./n + (x-mean_x)**2/Lxx)
# 预测区间
se_predict = lambda x: sd_res * np.sqrt(1+1./n + (x-mean_x)**2/Lxx)

可视化

plt.figure()      
plt.plot(x, fit(x),'k', label='Regression line')
plt.plot(x,y,'k.')
        
x.sort()
# 置信度
limit = (1-alpha)*100

# 置信区间范围
plt.plot(x, fit(x)+tval*se_fit(x), 'r--', lw=2, label='Confidence limit ({}%)'.format(limit))
plt.plot(x, fit(x)-tval*se_fit(x), 'r--', lw=2 )
        
# 预测区间范围
plt.plot(x, fit(x)+tval*se_predict(x), 'b*-', lw=2,  label='Prediction limit ({}%)'.format(limit))
plt.plot(x, fit(x)-tval*se_predict(x), 'b*-', lw=2)

plt.xlabel('X values')
plt.ylabel('Y values')
plt.title('Linear regression and confidence limits')


plt.legend(loc=0)
plt.show()

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-01-06，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自数据科学CLUB 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度