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

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

线性回归的定义

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

                                     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()

结果图:

 

当然我们可以调整不同参数,看运行的效果,图中的黑线是变化情况,最终我们也可得到一个近似的解,可以看出并不是我们上面直接解方程的解。

总结

    理解了一元的线性回归,其实按照一元的线性回归的思想,就比较容易理多元的线性回归,并且我们看到随机梯度算法,批量梯度算法,mini-batch,都是梯度算法改进,梯度下降算法理解了,这些理解起来也比较容易了。本文只是回归的入门知识,如有不足,欢迎指正。

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值