python多元函数 python多元函数求积分

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

Python怎么做最优化

最优化

为什么要做最优化呢？因为在生活中，人们总是希望幸福值或其它达到一个极值，比如做生意时希望成本最小，收入最大，所以在很多商业情境中，都会遇到求极值的情况。

函数求根

这里「函数的根」也称「方程的根」，或「函数的零点」。

先把我们需要的包加载进来。import numpy as npimport scipy as spimport scipy.optimize as optimport matplotlib.pyplot as plt%matplotlib inline

函数求根和最优化的关系？什么时候函数是最小值或最大值？

两个问题一起回答：最优化就是求函数的最小值或最大值，同时也是极值，在求一个函数最小值或最大值时，它所在的位置肯定是导数为 0 的位置，所以要求一个函数的极值，必然要先求导，使其为 0，所以函数求根就是为了得到最大值最小值。

scipy.optimize 有什么方法可以求根？

可以用 scipy.optimize 中的 bisect 或 brentq 求根。f = lambda x: np.cos(x) - x # 定义一个匿名函数x = np.linspace(-5, 5, 1000) # 先生成 1000 个 xy = f(x) # 对应生成 1000 个 f(x)plt.plot(x, y); # 看一下这个函数长什么样子plt.axhline(0, color='k'); # 画一根横线，位置在 y=0

opt.bisect(f, -5, 5) # 求取函数的根0.7390851332155535plt.plot(x, y)plt.axhline(0, color='k')plt.scatter([_], [0], c='r', s=100); # 这里的 [_] 表示上一个 Cell 中的结果，这里是 x 轴上的位置，0 是 y 上的位置

求根有两种方法，除了上面介绍的 bisect，还有 brentq，后者比前者快很多。%timeit opt.bisect(f, -5, 5)%timeit opt.brentq(f, -5, 5)10000 loops, best of 3: 157 s per loopThe slowest run took 11.65 times longer than the fastest. This could mean that an intermediate result is being cached.10000 loops, best of 3: 35.9 s per loop

函数求最小化

求最小值就是一个最优化问题。求最大值时只需对函数做一个转换，比如加一个负号，或者取倒数，就可转成求最小值问题。所以两者是同一问题。

初始值对最优化的影响是什么？

举例来说，先定义个函数。f = lambda x: 1-np.sin(x)/xx = np.linspace(-20., 20., 1000)y = f(x)

当初始值为 3 值，使用 minimize 函数找到最小值。minimize 函数是在新版的 scipy 里，取代了以前的很多最优化函数，是个通用的接口，背后是很多方法在支撑。x0 = 3xmin = opt.minimize(f, x0).x # x0 是起始点，起始点最好离真正的最小值点不要太远plt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300); # 起始点画出来，用圆圈表示plt.scatter(xmin, f(xmin), marker='v', s=300); # 最小值点画出来，用三角表示plt.xlim(-20, 20);

初始值为 3 时，成功找到最小值。

现在来看看初始值为 10 时，找到的最小值点。x0 = 10xmin = opt.minimize(f, x0).xplt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300)plt.scatter(xmin, f(xmin), marker='v', s=300)plt.xlim(-20, 20);

由上图可见，当初始值为 10 时，函数找到的是局部最小值点，可见 minimize 的默认算法对起始点的依赖性。

那么怎么才能不管初始值在哪个位置，都能找到全局最小值点呢？

如何找到全局最优点？

可以使用 basinhopping 函数找到全局最优点，相关背后算法，可以看帮助文件，有提供论文的索引和出处。

我们设初始值为 10 看是否能找到全局最小值点。x0 = 10from scipy.optimize import basinhoppingxmin = basinhopping(f,x0,stepsize = 5).xplt.plot(x, y);plt.scatter(x0, f(x0), marker='o', s=300);plt.scatter(xmin, f(xmin), marker='v', s=300);plt.xlim(-20, 20);

当起始点在比较远的位置，依然成功找到了全局最小值点。

如何求多元函数最小值？

以二元函数为例，使用 minimize 求对应的最小值。def g(X): x,y = X return (x-1)**4 + 5 * (y-1)**2 - 2*x*yX_opt = opt.minimize(g, (8, 3)).x # (8,3) 是起始点print X_opt[ 1.88292611 1.37658521]fig, ax = plt.subplots(figsize=(6, 4)) # 定义画布和图形x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, g((X, Y)), 50) # 等高线图ax.plot(X_opt[0], X_opt[1], 'r*', markersize=15) # 最小点的位置是个元组ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax) # colorbar 表示颜色越深，高度越高fig.tight_layout()

画3D 图。from mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import cmfig = plt.figure()ax = fig.gca(projection='3d')x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)surf = ax.plot_surface(X, Y, g((X,Y)), rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)cset = ax.contour(X, Y, g((X,Y)), zdir='z',offset=-5, cmap=cm.coolwarm)fig.colorbar(surf, shrink=0.5, aspect=5);

曲线拟合

曲线拟合和最优化有什么关系？

曲线拟合的问题是，给定一组数据，它可能是沿着一条线散布的，这时要找到一条最优的曲线来拟合这些数据，也就是要找到最好的线来代表这些点，这里的最优是指这些点和线之间的距离是最小的，这就是为什么要用最优化问题来解决曲线拟合问题。

举例说明，给一些点，找到一条线，来拟合这些点。

先给定一些点：N = 50 # 点的个数m_true = 2 # 斜率b_true = -1 # 截距dy = 2.0 # 误差np.random.seed(0)xdata = 10 * np.random.random(N) # 50 个 x，服从均匀分布ydata = np.random.normal(b_true + m_true * xdata, dy) # dy 是标准差plt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');

上面的点整体上呈现一个线性关系，要找到一条斜线来代表这些点，这就是经典的一元线性回归。目标就是找到最好的线，使点和线的距离最短。要优化的函数是点和线之间的距离，使其最小。点是确定的，而线是可变的，线是由参数值，斜率和截距决定的，这里就是要通过优化距离找到最优的斜率和截距。

点和线的距离定义如下：def chi2(theta, x, y): return np.sum(((y - theta[0] - theta[1] * x)) ** 2)

上式就是误差平方和。

误差平方和是什么？有什么作用？

误差平方和公式为：

误差平方和大，表示真实的点和预测的线之间距离太远，说明拟合得不好，最好的线，应该是使误差平方和最小，即最优的拟合线，这里是条直线。

误差平方和就是要最小化的目标函数。

找到最优的函数，即斜率和截距。theta_guess = [0, 1] # 初始值theta_best = opt.minimize(chi2, theta_guess, args=(xdata, ydata)).xprint(theta_best)[-1.01442005 1.93854656]

上面两个输出即是预测的直线斜率和截距，我们是根据点来反推直线的斜率和截距，那么真实的斜率和截距是多少呢？-1 和 2，很接近了，差的一点是因为有噪音的引入。xfit = np.linspace(0, 10)yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');

最小二乘（Least Square）是什么？

上面用的是 minimize 方法，这个问题的目标函数是误差平方和，这就又有一个特定的解法，即最小二乘。

最小二乘的思想就是要使得观测点和估计点的距离的平方和达到最小，这里的“二乘”指的是用平方来度量观测点与估计点的远近（在古汉语中“平方”称为“二乘”），“最小”指的是参数的估计值要保证各个观测点与估计点的距离的平方和达到最小。

关于最小二乘估计的计算，涉及更多的数学知识，这里不想详述，其一般的过程是用目标函数对各参数求偏导数，并令其等于 0，得到一个线性方程组。具体推导过程可参考斯坦福机器学习讲义第 7 页。def deviations(theta, x, y): return (y - theta[0] - theta[1] * x)theta_best, ier = opt.leastsq(deviations, theta_guess, args=(xdata, ydata))print(theta_best)[-1.01442016 1.93854659]

最小二乘 leastsq 的结果跟 minimize 结果一样。注意 leastsq 的第一个参数不再是误差平方和 chi2，而是误差本身 deviations，即没有平方，也没有和。yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');

非线性最小二乘

上面是给一些点，拟合一条直线，拟合一条曲线也是一样的。def f(x, beta0, beta1, beta2): # 首先定义一个非线性函数，有 3 个参数 return beta0 + beta1 * np.exp(-beta2 * x**2)beta = (0.25, 0.75, 0.5) # 先猜 3 个 betaxdata = np.linspace(0, 5, 50)y = f(xdata, *beta)ydata = y + 0.05 * np.random.randn(len(xdata)) # 给 y 加噪音def g(beta): return ydata - f(xdata, *beta) # 真实 y 和预测值的差，求最优曲线时要用到beta_start = (1, 1, 1)beta_opt, beta_cov = opt.leastsq(g, beta_start)print beta_opt # 求到的 3 个最优的 beta 值[ 0.25525709 0.74270226 0.54966466]

拿估计的 beta_opt 值跟真实的 beta = (0.25, 0.75, 0.5) 值比较，差不多。fig, ax = plt.subplots()ax.scatter(xdata, ydata) # 画点ax.plot(xdata, y, 'r', lw=2) # 真实值的线ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2) # 拟合的线ax.set_xlim(0, 5)ax.set_xlabel(r"$x$", fontsize=18)ax.set_ylabel(r"$f(x, \beta)$", fontsize=18)fig.tight_layout()

除了使用最小二乘，还可以使用曲线拟合的方法，得到的结果是一样的。beta_opt, beta_cov = opt.curve_fit(f, xdata, ydata)print beta_opt[ 0.25525709 0.74270226 0.54966466]

有约束的最小化

有约束的最小化是指，要求函数最小化之外，还要满足约束条件，举例说明。

边界约束def f(X): x, y = X return (x-1)**2 + (y-1)**2 # 这是一个碗状的函数x_opt = opt.minimize(f, (0, 0), method='BFGS').x # 无约束最优化

假设有约束条件，x 和 y 要在一定的范围内，如 x 在 2 到 3 之间，y 在 0 和 2 之间。bnd_x1, bnd_x2 = (2, 3), (0, 2) # 对自变量的约束x_cons_opt = opt.minimize(f, np.array([0, 0]), method='L-BFGS-B', bounds=[bnd_x1, bnd_x2]).x # bounds 矩形约束fig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X,Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 没有约束下的最小值，蓝色五角星ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 有约束下的最小值，红色星星bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]), bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0], facecolor="grey")ax.add_patch(bound_rect)ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()

不等式约束

介绍下相关理论，先来看下存在等式约束的极值问题求法，比如下面的优化问题。

目标函数是 f(w)，下面是等式约束，通常解法是引入拉格朗日算子，这里使用 ββ 来表示算子，得到拉格朗日公式为

l 是等式约束的个数。

然后分别对 w 和ββ 求偏导，使得偏导数等于 0，然后解出 w 和βiβi，至于为什么引入拉格朗日算子可以求出极值，原因是 f(w) 的 dw 变化方向受其他不等式的约束，dw的变化方向与f(w)的梯度垂直时才能获得极值，而且在极值处，f(w) 的梯度与其他等式梯度的线性组合平行，因此他们之间存在线性关系。（参考《最优化与KKT条件》）

对于不等式约束的极值问题

常常利用拉格朗日对偶性将原始问题转换为对偶问题，通过解对偶问题而得到原始问题的解。该方法应用在许多统计学习方法中。有兴趣的可以参阅相关资料，这里不再赘述。def f(X): return (X[0] - 1)**2 + (X[1] - 1)**2def g(X): return X[1] - 1.75 - (X[0] - 0.75)**4x_opt = opt.minimize(f, (0, 0), method='BFGS').xconstraints = [dict(type='ineq', fun=g)] # 约束采用字典定义，约束方式为不等式约束，边界用 g 表示x_cons_opt = opt.minimize(f, (0, 0), method='SLSQP', constraints=constraints).xfig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X, Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 蓝色星星，没有约束下的最小值ax.plot(x_, 1.75 + (x_-0.75)**4, '', markersize=15)ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color="grey")ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 在区域约束下的最小值ax.set_ylim(-1, 3)ax.set_xlabel(r"$x_0$", fontsize=18)ax.set_ylabel(r"$x_1$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()

scipy.optimize.minimize 中包括了多种最优化算法，每种算法使用范围不同，详细参考官方文档。

相关问答

Q1: python数据分析与应用第三章代码3-5的数据哪来的

savetxt

import numpy as np

i2 = np.eye(2)

np.savetxt("eye.txt", i2)

3.4 读入CSV文件

# AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800

c,v=np.loadtxt('data.csv', delimiter=',', usecols=(6,7), unpack=True) #index从0开始

3.6.1 算术平均值

np.mean(c) = np.average(c)

3.6.2 加权平均值

t = np.arange(len(c))

np.average(c, weights=t)

3.8 极值

np.min(c)

np.max(c)

np.ptp(c) 最大值与最小值的差值

3.10 统计分析

np.median(c) 中位数

np.msort(c) 升序排序

np.var(c) 方差

3.12 分析股票收益率

np.diff(c) 可以返回一个由相邻数组元素的差

值构成的数组

returns = np.diff( arr ) / arr[ : -1] #diff返回的数组比收盘价数组少一个元素

np.std(c) 标准差

对数收益率

logreturns = np.diff( np.log(c) ) #应检查输入数组以确保其不含有零和负数

where 可以根据指定的条件返回所有满足条件的数

组元素的索引值。

posretindices = np.where(returns 0)

np.sqrt(1./252.) 平方根，浮点数

3.14 分析日期数据

# AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800

dates, close=np.loadtxt('data.csv', delimiter=',', usecols=(1,6), converters={1:datestr2num}, unpack=True)

print "Dates =", dates

def datestr2num(s):

return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()

# 星期一 0

# 星期二 1

# 星期三 2

# 星期四 3

# 星期五 4

# 星期六 5

# 星期日 6

#output

Dates = [ 4. 0. 1. 2. 3. 4. 0. 1. 2. 3. 4. 0. 1. 2. 3. 4. 1. 2. 4. 0. 1. 2. 3. 4. 0.

1. 2. 3. 4.]

averages = np.zeros(5)

for i in range(5):

indices = np.where(dates == i)

prices = np.take(close, indices) #按数组的元素运算,产生一个数组作为输出。

a = [4, 3, 5, 7, 6, 8]

indices = [0, 1, 4]

np.take(a, indices)

array([4, 3, 6])

np.argmax(c) #返回的是数组中最大元素的索引值

np.argmin(c)

3.16 汇总数据

# AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800

#得到第一个星期一和最后一个星期五

first_monday = np.ravel(np.where(dates == 0))[0]

last_friday = np.ravel(np.where(dates == 4))[-1]

#创建一个数组，用于存储三周内每一天的索引值

weeks_indices = np.arange(first_monday, last_friday + 1)

#按照每个子数组5个元素，用split函数切分数组

weeks_indices = np.split(weeks_indices, 5)

#output

[array([1, 2, 3, 4, 5]), array([ 6, 7, 8, 9, 10]), array([11,12, 13, 14, 15])]

weeksummary = np.apply_along_axis(summarize, 1, weeks_indices,open, high, low, close)

def summarize(a, o, h, l, c): #open, high, low, close

monday_open = o[a[0]]

week_high = np.max( np.take(h, a) )

week_low = np.min( np.take(l, a) )

friday_close = c[a[-1]]

return("APPL", monday_open, week_high, week_low, friday_close)

np.savetxt("weeksummary.csv", weeksummary, delimiter=",", fmt="%s") #指定了文件名、需要保存的数组名、分隔符(在这个例子中为英文标点逗号)以及存储浮点数的格式。

0818b9ca8b590ca3270a3433284dd417.png

格式字符串以一个百分号开始。接下来是一个可选的标志字符：-表示结果左对齐，0表示左端补0，+表示输出符号(正号+或负号-)。第三部分为可选的输出宽度参数，表示输出的最小位数。第四部分是精度格式符，以”.”开头，后面跟一个表示精度的整数。最后是一个类型指定字符，在例子中指定为字符串类型。

numpy.apply_along_axis(func1d, axis, arr, *args, **kwargs)

def my_func(a):

... """Average first and last element of a 1-D array"""

... return (a[0] + a[-1]) * 0.5

b = np.array([[1,2,3], [4,5,6], [7,8,9]])

np.apply_along_axis(my_func, 0, b) #沿着X轴运动，取列切片

array([ 4., 5., 6.])

np.apply_along_axis(my_func, 1, b) #沿着y轴运动，取行切片

array([ 2., 5., 8.])

b = np.array([[8,1,7], [4,3,9], [5,2,6]])

np.apply_along_axis(sorted, 1, b)

array([[1, 7, 8],

[3, 4, 9],

[2, 5, 6]])

3.20 计算简单移动平均线

(1) 使用ones函数创建一个长度为N的元素均初始化为1的数组，然后对整个数组除以N，即可得到权重。如下所示：

N = int(sys.argv[1])

weights = np.ones(N) / N

print "Weights", weights

在N = 5时，输出结果如下：

Weights [ 0.2 0.2 0.2 0.2 0.2] #权重相等

(2) 使用这些权重值，调用convolve函数：

c = np.loadtxt('data.csv', delimiter=',', usecols=(6,),unpack=True)

sma = np.convolve(weights, c)[N-1:-N+1] #卷积是分析数学中一种重要的运算，定义为一个函数与经过翻转和平移的另一个函数的乘积的积分。

t = np.arange(N - 1, len(c)) #作图

plot(t, c[N-1:], lw=1.0)

plot(t, sma, lw=2.0)

show()

3.22 计算指数移动平均线

指数移动平均线(exponential moving average)。指数移动平均线使用的权重是指数衰减的。对历史上的数据点赋予的权重以指数速度减小，但永远不会到达0。

x = np.arange(5)

print "Exp", np.exp(x)

#output

Exp [ 1. 2.71828183 7.3890561 20.08553692 54.59815003]

Linspace 返回一个元素值在指定的范围内均匀分布的数组。

print "Linspace", np.linspace(-1, 0, 5) #起始值、终止值、可选的元素个数

#output

Linspace [-1. -0.75 -0.5 -0.25 0. ]

(1)权重计算

N = int(sys.argv[1])

weights = np.exp(np.linspace(-1. , 0. , N))

(2)权重归一化处理

weights /= weights.sum()

print "Weights", weights

#output

Weights [ 0.11405072 0.14644403 0.18803785 0.24144538 0.31002201]

(3)计算及作图

c = np.loadtxt('data.csv', delimiter=',', usecols=(6,),unpack=True)

ema = np.convolve(weights, c)[N-1:-N+1]

t = np.arange(N - 1, len(c))

plot(t, c[N-1:], lw=1.0)

plot(t, ema, lw=2.0)

show()

3.26 用线性模型预测价格

(x, residuals, rank, s) = np.linalg.lstsq(A, b) #系数向量x、一个残差数组、A的秩以及A的奇异值

print x, residuals, rank, s

#计算下一个预测值

print np.dot(b, x)

3.28 绘制趋势线

x = np.arange(6)

x = x.reshape((2, 3))

array([[0, 1, 2], [3, 4, 5]])

np.ones_like(x) #用1填充数组

array([[1, 1, 1], [1, 1, 1]])

类似函数

zeros_like

empty_like

zeros

ones

empty

3.30 数组的修剪和压缩

a = np.arange(5)

print "a =", a

print "Clipped", a.clip(1, 2) #将所有比给定最大值还大的元素全部设为给定的最大值，而所有比给定最小值还小的元素全部设为给定的最小值

#output

a = [0 1 2 3 4]

Clipped [1 1 2 2 2]

a = np.arange(4)

print a

print "Compressed", a.compress(a 2) #返回一个根据给定条件筛选后的数组

#output

[0 1 2 3]

Compressed [3]

b = np.arange(1, 9)

print "b =", b

print "Factorial", b.prod() #输出数组元素阶乘结果

#output

b = [1 2 3 4 5 6 7 8]

Factorial 40320

print "Factorials", b.cumprod()

#output

Q2: python中scipy包中的optimize里面的函数具体怎么用

from scipy.optimize import fmin

def myfunc(x):

return x**2-4*x+8

print fmin(myfunc, 0)

def myfunc(p):

x, y = p

return x**2+y**2+8

print fmin(myfunc, (1, 1))

复制代码

fmin的第一个参数是一个函数，这个函数的参数是一个数组，数组中每个元素是一个变量，因此对于多元函数，需要在myfunc内部将数组的内容展开。

Q3: python多元线性回归怎么计算

1、什么是多元线性回归模型？

当y值的影响因素不唯一时,采用多元线性回归模型。

y =y=β0+β1x1+β2x2+...+βnxn

例如商品的销售额可能不电视广告投入,收音机广告投入,报纸广告投入有关系,可以有 sales =β0+β1*TV+β2* radio+β3*newspaper.

2、使用pandas来读取数据

pandas 是一个用于数据探索、数据分析和数据处理的python库

[python] view plain copy

import pandas as pd

[html] view plain copy

pre name="code" class="python"# read csv file directly from a URL and save the results

data = pd.read_csv('/home/lulei/Advertising.csv')

# display the first 5 rows

data.head()

上面代码的运行结果：

TV Radio Newspaper Sales

0 230.1 37.8 69.2 22.1

1 44.5 39.3 45.1 10.4

2 17.2 45.9 69.3 9.3

3 151.5 41.3 58.5 18.5

4 180.8 10.8 58.4 12.9

上面显示的结果类似一个电子表格，这个结构称为Pandas的数据帧(data frame)，类型全称：pandas.core.frame.DataFrame.

pandas的两个主要数据结构：Series和DataFrame：

Series类似于一维数组，它有一组数据以及一组与之相关的数据标签(即索引)组成。

DataFrame是一个表格型的数据结构，它含有一组有序的列，每列可以是不同的值类型。DataFrame既有行索引也有列索引，它可以被看做由Series组成的字典。

[python] view plain copy

# display the last 5 rows

data.tail()

只显示结果的末尾5行

TV Radio Newspaper Sales

195 38.2 3.7 13.8 7.6

196 94.2 4.9 8.1 9.7

197 177.0 9.3 6.4 12.8

198 283.6 42.0 66.2 25.5

199 232.1 8.6 8.7 13.4

[html] view plain copy

# check the shape of the DataFrame(rows, colums)

data.shape

查看DataFrame的形状,注意第一列的叫索引，和数据库某个表中的第一列类似。

(200,4)

3、分析数据

特征：

TV：对于一个给定市场中单一产品，用于电视上的广告费用（以千为单位）

Radio：在广播媒体上投资的广告费用

Newspaper：用于报纸媒体的广告费用

响应：

Sales：对应产品的销量

在这个案例中，我们通过不同的广告投入，预测产品销量。因为响应变量是一个连续的值，所以这个问题是一个回归问题。数据集一共有200个观测值，每一组观测对应一个市场的情况。

注意：这里推荐使用的是seaborn包。网上说这个包的数据可视化效果比较好看。其实seaborn也应该属于matplotlib的内部包。只是需要再次的单独安装。

[python] view plain copy

import seaborn as sns

import matplotlib.pyplot as plt

# visualize the relationship between the features and the response using scatterplots

sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8)

plt.show()#注意必须加上这一句，否则无法显示。

[html] view plain copy

这里选择TV、Radio、Newspaper 作为特征，Sales作为观测值

[html] view plain copy

返回的结果：

seaborn的pairplot函数绘制X的每一维度和对应Y的散点图。通过设置size和aspect参数来调节显示的大小和比例。可以从图中看出，TV特征和销量是有比较强的线性关系的，而Radio和Sales线性关系弱一些，Newspaper和Sales线性关系更弱。通过加入一个参数kind='reg'，seaborn可以添加一条最佳拟合直线和95%的置信带。

[python] view plain copy

sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8, kind='reg')

plt.show()

结果显示如下：

4、线性回归模型

优点：快速；没有调节参数；可轻易解释；可理解。

缺点：相比其他复杂一些的模型，其预测准确率不是太高，因为它假设特征和响应之间存在确定的线性关系，这种假设对于非线性的关系，线性回归模型显然不能很好的对这种数据建模。

线性模型表达式： y=β0+β1x1+β2x2+...+βnxn 其中

y是响应

β0是截距

β1是x1的系数，以此类推

在这个案例中： y=β0+β1∗TV+β2∗Radio+...+βn∗Newspaper

(1)、使用pandas来构建X(特征向量)和y(标签列)

scikit-learn要求X是一个特征矩阵，y是一个NumPy向量。

pandas构建在NumPy之上。

因此，X可以是pandas的DataFrame，y可以是pandas的Series，scikit-learn可以理解这种结构。

[python] view plain copy

#create a python list of feature names

feature_cols = ['TV', 'Radio', 'Newspaper']

# use the list to select a subset of the original DataFrame

X = data[feature_cols]

# equivalent command to do this in one line

X = data[['TV', 'Radio', 'Newspaper']]

# print the first 5 rows

print X.head()

# check the type and shape of X

print type(X)

print X.shape

输出结果如下：

TV Radio Newspaper

0 230.1 37.8 69.2

1 44.5 39.3 45.1

2 17.2 45.9 69.3

3 151.5 41.3 58.5

4 180.8 10.8 58.4

class 'pandas.core.frame.DataFrame'

(200, 3)

[python] view plain copy

# select a Series from the DataFrame

y = data['Sales']

# equivalent command that works if there are no spaces in the column name

y = data.Sales

# print the first 5 values

print y.head()

输出的结果如下：

0 22.1

1 10.4

2 9.3

3 18.5

4 12.9

Name: Sales

（2）、构建训练集与测试集

[html] view plain copy

pre name="code" class="python"span style="font-size:14px;"##构造训练集和测试集

from sklearn.cross_validation import train_test_split #这里是引用了交叉验证

X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)

#default split is 75% for training and 25% for testing

[html] view plain copy

print X_train.shape

print y_train.shape

print X_test.shape

print y_test.shape

输出结果如下：

(150, 3)

(150,)

(50, 3)

(50,)

注：上面的结果是由train_test_spilit()得到的，但是我不知道为什么我的版本的sklearn包中居然报错：

ImportError Traceback (most recent call last)ipython-input-182-3eee51fcba5a in module() 1 ###构造训练集和测试集---- 2 from sklearn.cross_validation import train_test_split 3 #import sklearn.cross_validation 4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1) 5 # default split is 75% for training and 25% for testingImportError: cannot import name train_test_split

处理方法：1、我后来重新安装sklearn包。再一次调用时就没有错误了。

2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)

（3）sklearn的线性回归

[html] view plain copy

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()

model=linreg.fit(X_train, y_train)

print model

print linreg.intercept_

print linreg.coef_

输出的结果如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)

2.66816623043

[ 0.04641001 0.19272538 -0.00349015]

[html] view plain copy

# pair the feature names with the coefficients

zip(feature_cols, linreg.coef_)

输出如下：

[('TV', 0.046410010869663267),

('Radio', 0.19272538367491721),

('Newspaper', -0.0034901506098328305)]

y=2.668+0.0464∗TV+0.192∗Radio-0.00349∗Newspaper

如何解释各个特征对应的系数的意义？

对于给定了Radio和Newspaper的广告投入，如果在TV广告上每多投入1个单位，对应销量将增加0.0466个单位。就是加入其它两个媒体投入固定，在TV广告上每增加1000美元（因为单位是1000美元），销量将增加46.6（因为单位是1000）。但是大家注意这里的newspaper的系数居然是负数，所以我们可以考虑不使用newspaper这个特征。这是后话，后面会提到的。

（4）、预测

[python] view plain copy

y_pred = linreg.predict(X_test)

print y_pred

[python] view plain copy

print type(y_pred)

输出结果如下：

[ 14.58678373 7.92397999 16.9497993 19.35791038 7.36360284

7.35359269 16.08342325 9.16533046 20.35507374 12.63160058

22.83356472 9.66291461 4.18055603 13.70368584 11.4533557

4.16940565 10.31271413 23.06786868 17.80464565 14.53070132

15.19656684 14.22969609 7.54691167 13.47210324 15.00625898

19.28532444 20.7319878 19.70408833 18.21640853 8.50112687

9.8493781 9.51425763 9.73270043 18.13782015 15.41731544

5.07416787 12.20575251 14.05507493 10.6699926 7.16006245

11.80728836 24.79748121 10.40809168 24.05228404 18.44737314

20.80572631 9.45424805 17.00481708 5.78634105 5.10594849]

type 'numpy.ndarray'

5、回归问题的评价测度

(1) 评价测度

对于分类问题，评价测度是准确率，但这种方法不适用于回归问题。我们使用针对连续数值的评价测度(evaluation metrics)。

这里介绍3种常用的针对线性回归的测度。

1)平均绝对误差(Mean Absolute Error, MAE)

(2)均方误差(Mean Squared Error, MSE)

(3)均方根误差(Root Mean Squared Error, RMSE)

这里我使用RMES。

[python] view plain copy

pre name="code" class="python"#计算Sales预测的RMSE

print type(y_pred),type(y_test)

print len(y_pred),len(y_test)

print y_pred.shape,y_test.shape

from sklearn import metrics

import numpy as np

sum_mean=0

for i in range(len(y_pred)):

sum_mean+=(y_pred[i]-y_test.values[i])**2

sum_erro=np.sqrt(sum_mean/50)

# calculate RMSE by hand

print "RMSE by hand:",sum_erro

最后的结果如下：

type 'numpy.ndarray' class 'pandas.core.series.Series'

50 50

(50,) (50,)

RMSE by hand: 1.42998147691

（2）做ROC曲线

[python] view plain copy

import matplotlib.pyplot as plt

plt.figure()

plt.plot(range(len(y_pred)),y_pred,'b',label="predict")

plt.plot(range(len(y_pred)),y_test,'r',label="test")

plt.legend(loc="upper right") #显示图中的标签

plt.xlabel("the number of sales")

plt.ylabel('value of sales')

plt.show()

显示结果如下：（红色的线是真实的值曲线，蓝色的是预测值曲线）

直到这里整个的一次多元线性回归的预测就结束了。

6、改进特征的选择

在之前展示的数据中，我们看到Newspaper和销量之间的线性关系竟是负关系（不用惊讶，这是随机特征抽样的结果。换一批抽样的数据就可能为正了），现在我们移除这个特征，看看线性回归预测的结果的RMSE如何？

依然使用我上面的代码，但只需修改下面代码中的一句即可：

[python] view plain copy

#create a python list of feature names

feature_cols = ['TV', 'Radio', 'Newspaper']

# use the list to select a subset of the original DataFrame

X = data[feature_cols]

# equivalent command to do this in one line

#X = data[['TV', 'Radio', 'Newspaper']]#只需修改这里即可pre name="code" class="python" style="font-size: 15px; line-height: 35px;"X = data[['TV', 'Radio']] #去掉newspaper其他的代码不变

# print the first 5 rowsprint X.head()# check the type and shape of Xprint type(X)print X.shape

最后的到的系数与测度如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)

2.81843904823

[ 0.04588771 0.18721008]

RMSE by hand: 1.28208957507

然后再次使用ROC曲线来观测曲线的整体情况。我们在将Newspaper这个特征移除之后，得到RMSE变小了，说明Newspaper特征可能不适合作为预测销量的特征，于是，我们得到了新的模型。我们还可以通过不同的特征组合得到新的模型，看看最终的误差是如何的。

备注：

之前我提到了这种错误：

注：上面的结果是由train_test_spilit()得到的，但是我不知道为什么我的版本的sklearn包中居然报错：

处理方法：1、我后来重新安装sklearn包。再一次调用时就没有错误了。

2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)