线性回归从一元线性回归入门

最新推荐文章于 2024-09-22 19:58:36 发布

chaoping315

最新推荐文章于 2024-09-22 19:58:36 发布

阅读量4.8k

点赞数 3

文章标签：线性回归机器学习梯度下降

本文链接：https://blog.csdn.net/chaoping315/article/details/81463544

版权

本文是从一元线性回归为为基础来理解线性回归，适合于线性回归的入门，让初学者对于线性回归有直观的理解。本文也是我对于线性回归算法入门学习，分享给大家。

线性回归的定义

回归是应用于表达输入变量与输出变量的关系，在输入变量与输出变量之间做一个映射，一般的线性回归公式为

$H(\theta_0,\theta_1\cdots\theta_n) = \theta_0 + \theta_1x_1+\cdots+\theta_nx_n \Longrightarrow H(\theta) = \theta^TX$

其中θ为权重参数，现在给定的数据是给定一些X与其对应的Y值，需要求解一组θ,将X映射到Y，在无误差的情况下，用n组已知的X与Y数据，用解方程的方式就可以解出θ。由于误差的存在，用n组数据解出的一组θ并不满足剩于组的数据，并且会一些矛盾解。现在将问题转为如何求解一组θ,让X最优映射到Y。下面从最简单的一元线性回归来理解说明回归的思想，如果一元线性回归的原理理解了，多维线性回归就理解了。

一元线性回归

一元线性回归，表示只有一个因变量，从数学公式来说就是一条直线，数学公式：

$y = ax + b => H(\theta_0,\theta_1) = \theta_0x_0 + \theta_1x_1 \text { }\text { } (x_0=1)$

我们线性的生成100组数据，给定，a,b的，将产生的值叠加噪声，现在通过这100组值，来求解出 $\theta_0,\theta_1$ 让 $\theta_0,\theta_1$ 与a,b接近。下面是程序生的100数据，后面所有的计算都是基于这100组数据产生的。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import *
sampleSize = 100
mu = 0
sigma = 0.01
np.random.seed(0)
g = np.random.normal(mu, sigma, sampleSize)
x = np.linspace(0,1,100)
a = 0.5
b = 0.5
y = x*a + b
print y
y2 = y + g
print y2
plt.plot(x, y)
plt.scatter(x, y2, color='red')
plt.grid()
plt.title("linear regression")
plt.show()

数据的图像：

从图上可以看出，点是分布在直线的两边，直线相当于是这些点的线性拟合，我们通过点来解求出这样的一点直线。给出这些点，我们可以画出很多直线，这时候我们就需要一个度量标准，来评判那条直线的优劣，这个标准也就是我们平时说的目标函数。我们经常采用的是误差平方和来作为目标函数，误差平方和越小，我们认为求出的 $\theta_0,\theta_1$ 越接近于理想值。数据表达式为： $J(\theta_0,\theta_1) = \frac{1}{2}\{ (\theta_1x^1_1+\theta_0-y^1)^2+ (\theta_1x^2_1+\theta_0-y^2)+\cdots+(\theta_1x^n_1+\theta_0-y^n)\} =\frac{1}{2}\sum_{i=1}^n(\theta_1x^i_1+\theta_0-y^i)^2 =\frac{1}{2}\sum_{i=1}^n\theta_1^2(x^i_1)^2+2\theta_0\theta_1x^i_1-2\theta_1x^i_1y^i-2\theta_0y^i+(y^i)^2$

我们用程序来实现，画出图形便于理解

theta0 = np.linspace(-20,20,100)
theta1 = np.linspace(-20,20,100)
J_theta0_theta1 = np.linspace(-1, 1, 100)
for i, t0 in enumerate(theta0):
    t1 = theta1[i]
    sum = 0
    for j,xj in enumerate(x):
        sum = sum + (t1*xj+t0-y2[j])*(t1*xj+t0-y2[j])
    J_theta0_theta1[i] = sum*0.5
print J_theta0_theta1


figure = plt.figure()

ax = Axes3D(figure)
ax.scatter(theta0, theta1, J_theta0_theta1,color='r')

sum_x = np.sum(x)
sum_x2 = np.dot(x, x.T)
sum_y = np.sum(y2)
sum_y2 = np.dot(y2,y2.T)
sum_xy = np.sum(np.multiply(x,y))
print sum_xy
theta0, theta1 = np.meshgrid(theta0, theta1)
#z = sum_x2*np.dot(theta1, theta1.T) + 2*sum_x*np.dot(theta0,theta1.T) - 2 * sum_x * sum_y * theta1 - 2 * sum_y * theta0 + np.dot(theta0,theta0.T)
#print z
#z1 = np.dot(x, theta1.T) + theta0.T - y2.T
z = np.empty_like(theta0)
def computeCost(X, y2, t0, t1):
    sum = 0
    for j, xj in enumerate(X):
        sum = sum + (t1 * xj + t0 - y2[j]) * (t1 * xj + t0 - y2[j])
    return sum*0.5

for (i, j) ,v in np.ndenumerate(z):
    z[i,j] = computeCost(x, y2, theta0[i, j], theta1[i, j])
print z
ax.plot_surface(theta0, theta1, z, rstride=1, cstride=1, cmap='rainbow')
ax.set_zlabel('J(theta0,theta1)')
ax.set_ylabel('theta1')
ax.set_xlabel('theta0')
plt.contour(theta0, theta1, z, 20, alpha=.75, cmap=plt.cm.hot)
#plt.contourf(theta0, theta1, z)
plt.show()

得到结果：

目录函数据我们画了对角线上的点跟等高线，从这个图上可以看出，有一个极值点，极值点也就是最小点。求极值的方式就是让

求导让导数为0的点:

$\begin{array}{l} J(\theta_0,\theta_1) = \frac{1}{2}\sum_{i=1}^n(\theta_1x^i_1+\theta_0-y^i)^2 \\ \frac{\partial J(\theta_0,\theta_1)}{\partial \theta_0} = \sum_{i=1}^n\theta_1x^i_1+\theta_0-y^i = 0 \\ \frac{\partial J(\theta_0,\theta_1)}{\partial \theta_1} = \sum_{i=1}^n(\theta_1x^i_1+\theta_0-y^i)x_1^i = 0 \end{array}$

接上面程序，求解一元线性方程：

from scipy.linalg import solve

print "theta1,theta0", solve([[sum_x,100],[sum_x2,sum_x]],[sum_y,sum_xy])

解出结果为：theta1,theta0 [ 0.49648258 0.50235679] ，这样就能通过直接求的方式得到所需的参数。

目标函数的矩阵形式

结合上面部分，结合矩阵知识对公式进一步推导：

$\begin{array}{l} \theta=[\theta_0,\theta_1] \\ X = [1,x_1] \\ H(\theta) = \theta^TX \\ J(\theta) = \sum_{i=1}^n(H_\theta(x_{i}) - y^i)^2 \\ = (Y-X\theta)^T(Y-X\theta) \\ =(Y^T-\theta^TX^T)(Y-X\theta) \\ =Y^TY - Y^TX\theta-\theta^TX^TY+\theta^TX^TX\theta \\ \frac{\partial J(\theta)}{\partial \theta} = 0 - (Y^TX)^T - X^TY + 2X^TX\theta \\ = 2X^TX\theta - 2X^TY \\ = 0 \\ \theta = (X^TX)^{-1}X^TY \end{array}$

从公式上看，我们可以直接求出θ的值，当然为了防止过拟合，会加入扰动因子，这部分就不在这讨论了。

梯度下降算法的实现

上面的一元线性回归中，对于目标函数，对于任意给定的 $\theta_0,\theta_1$ ,如果沿着它们的负梯度方向衰减，衰减的步长足够小，就可以得到接近函数极值的点。对于一元线性回归的部分

$\begin{array}{l} x_0 = 1 \\ \theta_0^\prime = \theta_0 - \alpha\sum_{i=1}^n(\theta_0x_0+\theta_1x_1-y^i)x_0^i \\ \theta_1^\prime = \theta_1 - \alpha\sum_{i=1}^n(\theta_0x_0+\theta_1x_1-y^i)x_1^i \end{array}$

程序部分：

# 利用梯度下降算法求解 theta0 theta1 初始值为1

# 迭代阀值，当两次迭代损失函数之差小于该阀值时停止迭代
epsilon = 0.0001

# 学习率
alpha = 0.001

theta0 = 20
theta1 = 1

n = len(x)
print "n=",n
J_theta = [0, 0]
iter_times = 0
g_x = []
g_y = []
g_z = []
while True:
    diff = [0, 0]

    for i in range(n):
        diff[0] += (theta0*1 + theta1*x[i] - y2[i])*1
        diff[1] += (theta0*1 + theta1*x[i] - y2[i])*x[i]

    theta0 -= alpha * diff[0]
    theta1 -= alpha * diff[1]
    g_x.append(theta0)
    g_y.append(theta1)

    J_theta[1] = 0
    for i in range(n):
        J_theta[1] += (theta0*1 + theta1*x[i] - y2[i])**2

    J_theta[1] = 0.5 * J_theta[1]
    g_z.append(J_theta[1])
    print "iter times %d ,theta0 = %s, theta1= %s, epsilon = %s " % (iter_times,theta0,theta1,abs(J_theta[0] - J_theta[1]))
    if abs(J_theta[0] - J_theta[1]) < epsilon:
        break
    else:
        J_theta[0] = J_theta[1]
    iter_times += 1
print "theta0,theta1", theta0,theta1
ax.plot(g_x,g_y,g_z,c='black')
plt.show()

结果图：