Python与线性回归

By | 2017年7月14日

目录:

  1. 问题阐述
  2. 最小二乘法与leastsq
  3. sklearn
  4. 拟合优度R^2

 

1.问题阐述

已知多个近似分布于同一直线上的点,想拟合出一个直线方程:设该直线方程为y=ax+b,调整参数a和b,使得所有点到该直线的距离的平方和最小,设此时满足要求的a=A,b=B,则直线方程为y=Ax+B。

(我有个问题,目标是使得所有点到该直线的距离的平方和最小,但是下面的error函数却是用 实际的y-求得的y,这个并不是点到直线的距离啊,这样会导致目标函数并不正确啊,这里我没有想通……)

2.最小二乘法与leastsq

scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

这里的 func 就是损失函数,类似于神经网络中的梯度下降法中的损失函数。

如果func函数如下:

func(params) = ydata - f(xdata, params)

那么目标函数就是:

sum((ydata - f(xdata, params))**2, axis=0)

我们要做的就是使目标函数最小化,也就是均方差最小。

举个例子:

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt

Xi = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])
Yi = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])


def func(p, x):
    k, b = p
    return k*x+b


def error(p, x, y, s):
    print s
    return func(p, x)-y


p0 = [100, 2]
#print( error(p0,Xi,Yi) )

s = "Test the number of iteration"
Para = leastsq(error, p0, args=(Xi, Yi, s))
k, b = Para[0]
print"k=", k, '\n', "b=", b

plt.figure(figsize=(15, 10))
plt.scatter(Xi, Yi, color="red", label="Sample Point", linewidth=3)
x = np.linspace(0, 10, 1000)
y = k*x+b
plt.plot(x, y, color="orange", label="Fitting Line", linewidth=2)
plt.legend()
plt.show()

figure_1

3.sklearn

3.1 Pipeline

Pipeline的目的就是当设置不同的参数时组合几个可以一起交叉验证的步骤。所以可以使用组合这几个步骤的名字和它们的属性参数(不过需要在参数前面加_来连接)。

Pipeline的参数是一个list,list里面是一个个(name,transform)格式的tuple。最后一个tuple是估计函数(就是我们训练的模型类型)。而前面的tuple就是交叉验证的步骤。

clf = Pipeline([('poly', PolynomialFeatures(degree=d)), ('linear', LinearRegression(fit_intercept=False))])

 

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from sklearn.pipeline import Pipeline
from scipy.stats import norm
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

degree = [1, 2, 10]   # the degree of parameter


def rmse(y_test, y):
    return sp.sqrt(sp.mean((y_test - y) ** 2))


def R2(y_test, y_true):
    return 1 - ((y_test - y_true)**2).sum() / ((y_true - y_true.mean())**2).sum()

Xi = np.arange(0, 1, 0.002)
Yi = norm.rvs(0, size=500, scale=0.1)
Yi = Yi + Xi**2
'''
Xi = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])
Yi = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])

Xi = np.arange(1, 17, 1)
Yi = np.array([4.00, 6.40, 8.00, 8.80, 9.22, 9.50, 9.70, 9.86, 10.00, 10.20, 10.32, 10.42, 10.50, 10.55, 10.58, 10.60])
'''
plt.figure(figsize=(8, 6))
plt.scatter(Xi, Yi, color="black", label="Sample Point", linewidth=0.3)

for d in degree:
    clf = Pipeline([('poly', PolynomialFeatures(degree=d)), ('linear', LinearRegression(fit_intercept=False))])
    clf.fit(Xi[:, np.newaxis], Yi)
    y_test = clf.predict(Xi[:, np.newaxis])   # 预测估计值
    print (clf.named_steps['linear'].coef_)
    print ('rmse=%.4f, R2=%.4f, clf.score=%.4f' % (rmse(y_test, Yi), R2(y_test, Yi), clf.score(Xi[:, np.newaxis], Yi)))
    plt.plot(Xi, y_test, linewidth=2)
plt.grid()
plt.legend(['1', '2', '10'], loc='upper left')
plt.show()

这里我用了三组数据进行测试,结果分别如下:

第二组数据在一次的条件下效果还行,次数越高越奇怪了。

figure_3

figure_2

figure_1

如果想用训练好的曲线求新的x对应的y值,就可以加类似的下面三句话:

OutofExpectation = 2.3
y = clf.predict(OutofExpectation)
print y[0]

微信截图_20170717163142

4.拟合优度R^2

做回归分析,常用的误差主要有均方误差根(RMSE)和R-平方(R2)。

RMSE是预测值与真实值的误差平方根的均值。这种度量方法很流行(Netflix机器学习比赛的评价方法),是一种定量的权衡方法。

R2方法是将预测值跟只使用均值的情况下相比,看能好多少。其区间通常在(0,1)之间。0表示还不如什么都不预测,直接取均值的情况,而1表示所有预测跟真实结果完美匹配的情况。

上面的代码已经包含了R^2的计算,这里只是单独将其定义式再拿出来:

def R2(y_test, y_true):
    return 1 - ((y_test - y_true)**2).sum() / ((y_true - y_true.mean())**2).sum()

上面三组数据的R^2结果分别如下:

微信截图_20170717154001

微信截图_201707171539142

微信截图_20170717154443

 

 

 

参考文献:

  1. Python闲谈(二)聊聊最小二乘法以及leastsq函数
  2. scipy.optimize.leastsq
  3. 用Python开始机器学习(3:数据拟合与广义线性回归)
  4. sklearn.pipeline.Pipeline类的用法
  5. python曲线拟合

 

发表评论

电子邮件地址不会被公开。 必填项已用*标注