一、线性回归理论概述
概述大致引用自机器学习之线性回归-AndrewNg学习笔记
注意，此文章的意义在于用符号计算的方式体现线性回归的每一个步骤，如何最快的实现可以查看其他线性回归教程。
回归问题是确定两种或两种以上变量相互依赖的定量关系的一种统计分析方法。统计学中，线性回归(Linear Regression)是利用称为线性回归方程的最小平方函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。
为了运行监督学习算法，首先要决定的是如何在计算机中表示假设函数(Hypothesis Function)h。在这里我们决定把房价表示为特征的线性函数：
假设函数
其中θ为参数(或者可以理解为权重)，依照惯例，x0 = 1，可以得到向量表示的形式：
假设函数向量表示
现在对于给定的训练集合，问题是我们要如何选取(学习)参数θ。一个合理的想法是使得h(x)的值尽量逼近真实标记y(起码在给定的训练集合上努力做到这一点)。我们定义如下的代价函数(它描述的是在任一θ下，h(x)和相应的y的接近程度)，代价函数如下：(那么现在我们的目标就是要使得如下的cost function的值最小)
代价函数
LMS算法：
我们的目标是找到使得J(θ)最小的参数向量θ，为了达到这样的目的，我们需要使用一种搜索算法：给定初始的θ值，并且不断改变θ的值是的J(θ)更小直到θ收敛。在这里我们使用的梯度下降算法：从初始的θ值开始按照以下的形式不断更新：
alpha称为学习速率或者步长。这个算法的基本思想是：要选择θ使得J(θ)达到最小，那么我每次沿着使J下降最快的方向走一步，而这个下降最快的方向就是梯度的方向。接下来是要进行偏导的求解。
LMS(least mean squares)更新规则为：
这就是原始的代价函数的梯度下降算法。需要注意的是：该算法在更新每一个θ时，都要遍历整个训练集。所以该算法也叫做批梯度下降算法。（需要说明的是，尽管梯度下降算法一般会受到局部最小值的影响，但我们在求解线性回归时做面临的最优化问题仅有一个全局最优解）
二、Python实现
需求第三方包：pandas、matplotlib、sympy
需求数据：吴恩达线性回归课后作业ex1data1.txt（可在Coursera的课后作业中下载）
首先读入数据
with open('ex1data1.txt',encoding='ascii') as f:
data = [list(map(float, line[:-1].split(','))) for line in f]
data[:5]
因为源数据没有题头，不能直接用pandas的read_csv读入，所以用python自带的io方法读取，line[:-1]是因为字符串line是以'\n'结尾的。将一行按','分开，并用map将每一个转换为float类型，结果返回map对象，再用list转换为list类型，最终data为二维列表。
import pandas as pd
df = pd.DataFrame(data, columns=['城市人口', '利润'])
df.head()
%matplotlib inline
import matplotlib.pyplot as plt
plt.scatter(*list(zip(*data)), s=30, marker="x", color='#ee0000')
引入pandas库，并将所有点画出来，注意，这个程序是运行在jupyter notebook中的，所以有%matplotlib inline这个魔术命令，如果是在Pycharm或其他IDE中，请去掉这一条，并在最后加上plt.show()
plt.scatter是散点图绘制函数，第一个参数为x轴所有点的列表，第二个参数为y轴所有点的列表，在程序中整体表现为*list(zip(*data))，s是大小，marker为标记的类型，color是用‘rgb’的方式表示的，详情可参阅matplotlib官方文档。
scatter中的*list(zip(*data))涉及到解包和压缩（转置）操作。首先zip函数常用来做转置操作，能把m*n的二维列表变成n*m的二维列表，例如zip([1,2,3,4],[5,6,7,8])，结果就是[(1,5),(2,6),(3,7),(4,8)]，而在这个程序中可以把n*2的二维列表，变成由x组成的列表1和y组成得到列表2满足plt.scatter前两个参数的需求。两个*的解包操作呢？首先，zip函数要求它的参数是几个长度相同的可迭代类型对象（列表、元组等），而传入的data是n个长度为2的列表没错，但他们整体被包含在一个列表中，会被认为是一个参数，因此用*操作符解开列表，才能用于zip；而前面的*号同理，不用解包操作会被误认为是一个参数交给第一个变量，而解包后，就会分别赋值给scatter前两个变量。
import sympy as sp
sp.init_printing(use_unicode=True)
alpha = 0.001
t0 = 0
t1 = 1
现在引用符号运算库sympy，这个库可以定义符号，求表达式值、求微积分等等，init_printing可以在jupyter中输出格式化的数学公式方便看
alpha是梯度下降步长，t0、t1为θ的初始值，可以在之后的梯度下降中多试验几次
stheta0 = sp.Symbol('theta0')
stheta1 = sp.Symbol('theta1')
sx = sp.Symbol('x')
def H(x=None):
res = stheta0 + sx * stheta1
if x is not None:
res = res.subs(sx, x)
return res
def J(df):
sum = 0
for i in range(len(df)):
row = df.ix[i]
sum += (H(x=row['城市人口'])-row['利润'])**2
return (1/(2*len(df))) * sum
定义θ0、θ1、x的符号并组成假设函数和代价函数，输出J(θ)可以看到没有用叠加符号的数学公式。
对于sp.Symbol('theta0')和它们乘积、加和的结果，可以看：
它们的类型为：
这些对象（Symbol、Add）都可以用一连串的subs(symbol, value)的方式或者字典{symbol:value}的方式代入值：
注意，当所有未知数都消去时（如In [13]所示），类型会变成sympy包中的一个整数（Integer）类型，如果需要转化，直接用python内置的int转化就好，同样的方法适用于numpy的int32等类型（之前刷leetcode题时用python，线下调试正确，上传却失败，明明支持numpy等第三方库，最后发现返回的类型不是python的int而是numpy的int32）
import matplotlib
#开启动态（交互式）
is_ipython = 'inline' in matplotlib.get_backend()
if is_ipython:
from IPython import display
plt.ion()
for _ in range(int(1/alpha)):
plt.cla()
plt.scatter(*list(zip(*data)), s=30, marker="x", color='#ee0000')
plt.plot([5,23], [t0+t1*5, t0+t1*23])
if is_ipython:
display.clear_output(wait=True)
display.display(plt.gcf())
temp_t0 = float(t0-alpha * (sp.diff(J(df), stheta0).subs(stheta0, t0).subs(stheta1, t1)))
temp_t1 = float(t1-alpha * (sp.diff(J(df), stheta1).subs(stheta0, t0).subs(stheta1, t1)))
t0 = temp_t0
t1 = temp_t1
前面是在jupyter开启动态展示的代码，从for开始迭代，我设定的迭代次数为步长的倒数，可以任意更改。
plt.cla()是清理画布，改变上次的绘图结果，显示最新的绘图结果
之后画点，并画出以当前θ0和θ1下5-23的线【为何取5-23？因为这是这组数据的最大、最小值，可根据数据自行更改】
画完当前的图后，开始梯度下降，先代入x值，求偏微分再代入θ0和θ1值；这里遵循了之前的重要原则，当把所有新的梯度下降值都求出来后再批量赋值
sp.diff是求微分的函数，而sp.diff(J(df), stheta1)是对θ1求偏微分的意思
会看到线一点点变化，变化越来越小看起来就像没动一样
三、总结
本程序是实现最基本的方法，正常运算时不会现求微分，而是直接套用LMS的更新规则，本例的目的是明确梯度下降法线性回归的步骤。
网友评论