注解
此笔记本可在此处下载： 01_Linear_regression_Linear_model.ipynb
线性回归-线性模型
代码作者：Emile Roux emile.roux@univ-smb.fr
RISE 幻灯片
范围
有限数 \(N\) 有个数据点可供选择：找到最适合的线路，解决这个问题。 \(N\) 点。
https://en.wikipedia.org/wiki/Linear_regression
规则
#setup
%load_ext autoreload
%matplotlib notebook
%autoreload 2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import ipywidgets as ipw
简单两点直线方程介绍
让我们来分析一下如何确定一个线性函数的简单方程 (\(y=ax+b\) ）Kowing线的牵引点。
# the known points
y=np.array([-1,2])
x=np.array([.5,3])
x,y
(array([0.5, 3. ]), array([-1, 2]))
fig = plt.figure()
plt.plot(x,y,'o')
plt.grid()
<IPython.core.display.Javascript object>
使 \(a\) 和 \(b\) 在系数方面，我们必须解决以下系统：
\[y_0=ax_0+b\]
\[y_1=ax_1+b\]
让我们重新组织一下这个系统：
\[b+ax_0=y_0\]
\[B+轴1=Y轴1\]
现在很容易把它写成线性系统：
\[\开始{bmatrix}\]
在这之前 \(a\) 和 \(b\) 系数可以通过反演线性系统来确定：
\[\开始{bmatrix}\]
让我们用python来做这个
# definition of the involved matrix :
X = np.array([[1, x[0]],[1,x[1]]])
print("X=",X)
Y = np.array([y[0],y[1]]).T
print("Y=",Y)
X= [[1. 0.5]
[1. 3. ]]
Y= [-1 2]
# Inversion of the systeme
coef=np.dot(np.linalg.inv(X),Y)
a = coef[1]
b = coef[0]
print("a=",a,)
print("b=",b)
a= 1.2000000000000002
b= -1.6
# draw the obtained line
N = 10
xmin, xmax = 0., 4
xi = np.linspace(xmin, xmax, N)
yi = a * xi + b
fig = plt.figure()
plt.plot(x,y,'o')
plt.plot(xi,yi,'r')
plt.grid()
<IPython.core.display.Javascript object>
如果我们添加第三个点，附加什么？
# the known points
y=np.array([-1,2,0])
x=np.array([.5,3,1.5])
x,y
(array([0.5, 3. , 1.5]), array([-1, 2, 0]))
使 \(a\) 和 \(b\) 在系数方面，我们必须解决以下系统：
\[y_0=ax_0+b\]
\[y_1=ax_1+b\]
\[y_2=ax_2+b\]
这个系统定义过度：我们有两个未知 (\(a\) ，。\(b\) ）和3个方程式。
现在的目标是找到 \(a\) 和 \(b\) 最适合给定数据的系数（即寻找线性方程）。
我们可以用以下线性系统的形式重写方程：
\[\开始{bmatrix}\]
该系统可以用性别表达：
\[x\β=y\]
矩阵 \(X\) 不是方形的，因此不能反转。上述方法不再适用。
诀窍是将系统乘以 \(X\) .
\[x^t x \β=x^t y\]
Therm公司 \(X^t X\) 现在是一个方阵。
让我们检查一下这个确认：
# Construcion of the X matrix
X = np.array([np.ones(x.size), x]).T
X
array([[1. , 0.5],
[1. , 3. ],
[1. , 1.5]])
# Computation of Xt X
np.dot(X.T,X)
array([[ 3. , 5. ],
[ 5. , 11.5]])
得到的矩阵是一个正方形矩阵，可以证明该矩阵是对称的，因而总是可逆的。
从上一个公式：
\[x^t x \β=x^t y\]
这个 \(\beta\) 矢量可以表示为：
\[\β=\左（x^t x\右）^-1 x^t y\]
这个算符是最小二乘估计量。
让我们用python来做这个
# Construcion of the X matrix
X = np.array([np.ones(x.size), x]).T
# Construcion of the Y matrix
Y = y.T
print("X=",X)
print("Y=",Y)
X= [[1. 0.5]
[1. 3. ]
[1. 1.5]]
Y= [-1 2 0]
# Computation of the least square estimator
beta = np.dot(np.linalg.inv(np.dot(X.T,X)),np.dot(X.T,Y))
print("beta=",beta)
print("The fitted linear function is:")
print("y=",beta[1],"x +",beta[0])
beta= [-1.68421053 1.21052632]
The fitted linear function is:
y= 1.210526315789474 x + -1.6842105263157898
现在该画线了
# draw the obtained line
N = 10
xmin, xmax = 0., 4
xi = np.linspace(xmin, xmax, N)
yi = beta[0] + beta[1] * xi
fig = plt.figure()
plt.plot(x,y,'o')
plt.plot(xi,yi,'r')
plt.grid()
<IPython.core.display.Javascript object>
有更大的数据集吗？
数据集：具有一定噪声的综合数据集
N = 100
xmin, xmax = 0., 4
x = np.linspace(xmin, xmax, N) + np.random.rand(N)
y = (3.3 * x -2) + 2*np.random.randn(N)
fig = plt.figure()
plt.plot(x,y,'o',label = "Data")
plt.legend()
plt.grid()
<IPython.core.display.Javascript object>
回归法确定系数
# Construcion of the X matrix
X = np.array([np.ones(x.size), x]).T
# Construcion of the Y matrix
Y = y.T
# Computation of the least square estimator
beta = np.dot(np.linalg.inv(np.dot(X.T,X)),np.dot(X.T,Y))
print("The fitted linear function is:")
print("y=",beta[1],"x +",beta[0])
The fitted linear function is:
y= 3.5040953300898092 x + -2.7670001724748126
绘制获得回归线
N = 10
xmin, xmax = 0., 5
xi = np.linspace(xmin, xmax, N)
yi = beta[0] + beta[1] * xi
fig = plt.figure()
plt.plot(x,y,'o',label = "Data")
plt.plot(xi,yi,'r',label = "Regression")
plt.legend()
plt.grid()
<IPython.core.display.Javascript object>
python知道如何在一行中做到这一点吗？
当然可以！
使用numpy的polyfit函数：
# compute the coefficient using the least square methode
beta=np.polyfit(x, y, 1)
print("beta=",beta)
# creat the associate 1d polynomila function
fit = np.poly1d(beta)
beta= [ 3.50409533 -2.76700017]
绘制obtaine回归线
N = 10
xmin, xmax = 0., 5
xi = np.linspace(xmin, xmax, N)
yi_numpy = fit(xi)
fig = plt.figure()
plt.plot(x,y,'o',label = "Data")
plt.plot(xi,yi,'r',label = "Regression")
plt.plot(xi,yi_numpy,'g',label = "Regression Numpy")
plt.legend()
plt.grid()
<IPython.core.display.Javascript object>
使用numpy linalg函数：
对于这种方法 \(X\) 自动导航系统 \(Y\) 矩阵必须由用户给出。
beta, r, rank, s = np.linalg.lstsq(X, Y,rcond=None)
print("beta=",beta)
beta= [-2.76700017 3.50409533]
幸运的是，所有的方法都给出了相同的结果！
噪声对由此产生的回归有什么影响？
N = 100
xmin, xmax = 0., 4
def get_data(N,sigma):
x = np.linspace(xmin, xmax, N)
y = (3.3 * x -2) + sigma*np.random.randn(N)
return x,y
def compute_reg(x,y):
beta,residuals, rank, singular_values, rcond =np.polyfit(x, y, 1, full=True)
fit= np.poly1d(beta)
SStot = ((y - y.mean())**2).sum()
SSres = residuals[0]
R2 = 1 - SSres/SStot
return beta,fit([xmin,xmax]),R2
x,y = get_data(100, 0.)
beta,fity, R2 = compute_reg(x,y)
fig = plt.figure()
points_plot,=plt.plot(x,y,'o',label = "Data")
line,=plt.plot([xmin,xmax],fity,'-r',label = "Regression")
plt.plot([xmin,xmax],fity,'-g',label = "Non Noisy data")
plt.grid()
plt.legend()
@ipw.interact(N=(5,100,10),sigma=(0.0,10.0,.1))
def update(N,sigma):
x,y = get_data(N,sigma)
beta,fity, R2 = compute_reg(x,y)
points_plot.set_xdata(x)
points_plot.set_ydata(y)
line.set_ydata(fity)
plt.title('$R^2$={:0.2f}'.format(R2))
<IPython.core.display.Javascript object>
interactive(children=(IntSlider(value=45, description='N', min=5, step=10), FloatSlider(value=5.0, description…