Question

让我们亲自体验一下，看一看我们假想的互联网创业公司MLAAS。它通过HTTP向用户推销机器学习算法服务。但随着公司不断取得成功，要为所有Web访问请求都提供优质服务，就需要具备更好的基础设施。我们并不愿意分配过多的资源，因为这些资源非常昂贵。另一方面，如果没有足够的资源来为所有请求提供服务，我们也将会赔钱。现在的问题是，我们何时会到达目前基础设施的极限。这个极限我们估计是每小时100 000个请求。我们希望事先知道什么时候不得不申请更多的云端服务器来服务于所有请求，同时不必为未使用的服务器承担费用。

1.5.1　读取数据

我们已经收集了上个月的Web统计信息，并把它们汇聚到了ch01/data/web_traffic.tsv（因为tsv包含以Tab字符分割的数字）。它们存储着每小时的访问次数。每一行包含连续的小时信息，以及该小时内的Web访问次数。

文件前面几行如下图所示：

使用SciPy的genfromtxt() 很容易读取数据。

import scipy as sp
data = sp.genfromtxt("web_traffic.tsv", delimiter="\t")

必须使用Tab作为分隔符，确保可以正确读取各个列的数据。

经快速检验显示，我们已经正确地读取了数据。

>>> print(data[:10])
[[ 1.00000000e+00 2.27200000e+03]
[ 2.00000000e+00 nan]
[ 3.00000000e+00 1.38600000e+03]
[ 4.00000000e+00 1.36500000e+03]
[ 5.00000000e+00 1.48800000e+03]
[ 6.00000000e+00 1.33700000e+03]
[ 7.00000000e+00 1.88300000e+03]
[ 8.00000000e+00 2.28300000e+03]
[ 9.00000000e+00 1.33500000e+03]
[ 1.00000000e+01 1.02500000e+03]]
>>> print(data.shape)
(743, 2)

我们有743个二维数据点。

1.5.2　预处理和清洗数据

在SciPy中，为了处理起来更加便利，我们将各维度分成两个向量，其中每个向量的大小是743。第一个向量x 包含小时信息，而另一个向量y 包含某个小时内的Web访问数。这个切分过程是通过我们选择的某些列，由SciPy的特殊索引标记来完成的。

x = data[:,0]
y = data[:,1]

注意 　这里还有很多从SciPy数组中选取数据的方法。可以在http://www.scipy.org/Tentative_NumPy_Tutorial 上查看更多关于索引（indexing）、切割（slicing）和迭代（iterating）的详情。

需要说明的是，y 中仍然有一些项包含了无效值nan 。但问题是，该如何处理这些无效值呢？让我们看一下有多少小时的数据中包含了无效值。

>>> sp.sum(sp.isnan(y))
8

我们看到743个项中只有8个值缺失了，因此把它们删除是可以承受的。记住，我们能够用另一个数组来索引SciPy的数组。sp.isnan(y) 返回一个布尔型的数组，用来表示某个数组项中的内容是否是一个数字。我们可以使用~ 在逻辑上对数组取反，使我们可以在x 和y 中只选择y 值合法的项。

x = x[~sp.isnan(y)]
y = y[~sp.isnan(y)]

为了获得对数据的第一印象，让我们用Matplotlib在散点图上将数据画出来。Matplotlib包含了pyplot 包。它模仿了Matlab的接口——一个非常方便和易用的接口。（你可以在http://matplotlib.org/users/pyplot_tutorial.html 上看到更多画图方面的教程。）

import matplotlib.pyplot as plt
plt.scatter(x,y)
plt.title("Web traffic over the last month")
plt.xlabel("Time")
plt.ylabel("Hits/hour")
plt.xticks([w*7*24 for w in range(10)],
   ['week %i'%w for w in range(10)])
plt.autoscale(tight=True)
plt.grid()
plt.show()

在绘出的图上，可以看到虽然前面几个星期的流量差不多相同，但最后那个星期呈现出显著上升的趋势。

1.5.3　选择正确的模型和学习算法

现在我们已经对数据有了一个初步印象，那么回到起初的问题：服务器要用多长时间来处理进来的Web流量呢？要回答这个问题，必须先做到以下两点：

找到有噪数据背后真正的模型；

使用这个模型预测未来，以便及时找到我们的基础设施必须扩展的地方。

1. 在构建第一个模型之前

谈到模型，你可以把它想象成对复杂现实世界的简化的理论近似。它总会包含一些劣质内容，而这又叫做近似误差。这个误差将指引我们在无数选择中寻找正确的模型。我们用模型预测值到真实值的平方距离来计算这个误差。具体来说，对于一个训练好的模型函数f ，按照下面这样来计算误差：

def error(f, x, y):
    return sp.sum((f(x)-y)**2)

向量x 和y 包含我们之前提取的Web统计数据。这正是Scipy向量化函数（这里采用的是f(x)）的美妙之处。在训练好的模型中，我们假定它把一个向量作为输入，并返回一个相同大小的向量。这样，我们就可以用它来计算与y之间的差距。

2. 从一条简单的直线开始

让我们假设另外一个例子，它的模型是一条直线。这里的挑战是如何在图中画出一条最佳的直线，使结果中的近似误差最小。SciPy的polyfit() 函数正是用来解决这个问题的。给定数据x 和y ，以及期望的多项式的阶（直线的阶是1），它可以找到一个模型，能够最小化之前定义的误差函数。

fp1, residuals, rank, sv, rcond = sp.polyfit(x, y, 1, full=True)

polyfit() 函数会把拟合的模型函数所使用的参数返回，即fp1 ；而且通过把full 置成True ，我们还可以获得更多逼近过程的背景信息。在这里面，我们只对残差感兴趣，而这正是近似误差。

>>> print("Model parameters: %s" % fp1)
Model parameters: [ 2.59619213 989.02487106]
>>> print(res)
[ 3.17389767e+08]

这里的意思是说，最优的近似直线如下面这个函数所示：

f(x) = 2.59619213 * x + 989.02487106

然后用poly1d() 根据这些参数创建一个模型函数。

>>> f1 = sp.poly1d(fp1)
>>> print(error(f1, x, y))
317389767.34

我们已经利用full=True 得到了更多的关于逼近过程的细节。正常来说，我们并不需要这个，只需要返回模型参数即可。

注意 　事实上，我们在这里只是做了曲线拟合。更多详细信息，请参考http://en.wikipedia.org/wiki/Curve_fitting 。

现在用f1() 画出第一个训练后的模型。在前述绘图命令之外，我们简单加入如下代码：

fx = sp.linspace(0,x[-1], 1000) # 生成X值用来作图
plt.plot(fx, f1(fx), linewidth=4)
plt.legend(["d=%i" % f1.order], loc="upper left")

下面这个图中显示了我们第一个训练后的模型：

虽然前面4个星期的数据好像并没有偏离太多，但我们仍然可以清楚地看到，最初的直线模型假设是有问题的。此外，实际误差值317 389 767.34到底是好还是坏呢？

我们从来不拿误差的绝对值单独使用。然而当比较两个竞争的模型时，可以利用它们的绝对误差来判断哪一个更好。尽管第一个模型显然不是我们想要的，但它的工作流程却有一个重要作用：我们可以把它当做基线，直到找到更好的模型。无论将来构造出了什么样的模型，我们都会去和当前的基线做比较。

3. 一些高级话题

现在我们要用一个更复杂的模型来做拟合，来看一个阶数为2的多项式，看看它是否可以更好地“理解”我们的数据。

>>> f2p = sp.polyfit(x, y, 2)
>>> print(f2p)
array([ 1.05322215e-02, -5.26545650e+00, 1.97476082e+03])
>>> f2 = sp.poly1d(f2p)
>>> print(error(f2, x, y))
179983507.878

下面这个图表显示了之前训练好的模型（一阶直线），以及我们新训练出的更复杂的二阶模型（虚线）：

这里的误差是179 983 507.878，几乎是直线模型误差的一半。这个效果看起来很不错，然而，它也是有代价的。我们现在得到了一个更复杂的函数，这意味着在polyfit() 中多了一个参数需要调整。近似的多项式如下：

f(x) = 0.0105322215 * x**2 - 5.26545650 * x + 1974.76082

在这里，如果复杂性越大效果越好，那么为什么不进一步增加复杂性呢？让我们试一下阶数为3、10和100的函数。

数据越复杂，曲线对数据逼近得越好。它们的误差值似乎也反映出了同样的结果。

Error d=1: 317,389,767.339778
Error d=2: 179,983,507.878179
Error d=3: 139,350,144.031725
Error d=10: 121,942,326.363461
Error d=100: 109,318,004.475556

然而，如果近距离观察拟合出的曲线，我们就会开始对它们能否捕捉到真实的数据生成过程心生疑虑。换句话说，我们的模型是否真正代表了广大客户访问我们网站的行为呢？看看10阶和100阶的多项式，我们发现了巨大的震荡。似乎这些模型对数据拟合得太过了。它不但捕捉到了背后的数据生成过程，还把噪声也包含进去了，这就叫做过拟合 （overfitting）。

在这里，我们有如下的选择。

选择其中一个拟合出的多项式模型。

换成另外一类更复杂的模型；样条（splines）？

从不同的角度思考数据，然后重新开始。

在上述这5个拟合模型中，1阶模型明显太过简单了，而10阶和100阶的模型显然是过拟合了。只有2阶和3阶模型似乎还比较匹配数据。然而，如果在数据的两个边界上进行预测，我们会发现它们的效果令人抓狂。

换成另外一类更复杂的模型似乎也是一个错误路线。那么什么样的论据会支持哪类模型呢？在这里，我们意识到，也许我们还没有真正理解数据。

4. 以退为进——另眼看数据

在此，我们退回去从另一个角度来看数据。似乎在第3周和第4周的数据之间有一个拐点。这让我们可以以3.5周作为分界点把数据分成两份，并训练出两条直线来。我们使用到第3周之前的数据来训练第一条线，用剩下的数据训练第2条线。

inflection = 3.5*7*24 # 计算拐点的小时数
xa = x[:inflection] # 拐点之前的数据
ya = y[:inflection]
xb = x[inflection:] # 之后的数据
yb = y[inflection:]
fa = sp.poly1d(sp.polyfit(xa, ya, 1))
fb = sp.poly1d(sp.polyfit(xb, yb, 1))
fa_error = error(fa, xa, ya)
fb_error = error(fb, xb, yb)
print("Error inflection=%f" % (fa + fb_error))
Error inflection=156,639,407.701523

我们在这两组数据的范围之内画出了这两个模型，如下图所示：

很明显，两条线组合起来似乎比之前我们做的任何模型都能更好地拟合数据。但组合之后的误差仍然高于高阶多项式的误差。我们最后能否相信这个误差呢？

换一个方式来问，相比于其他复杂模型，为什么仅在最后一周数据上更相信拟合的直线模型呢？这是因为我们认为它更符合未来数据。如果在未来时间段上画出模型，就可以看到这是非常正确的（d=1 是我们最初的直线模型）。

10阶和100阶的模型在这里似乎并没有什么光明的未来。它们非常努力地对给定数据正确建模，但它们明显没法推广到将来的数据上。这个叫做过拟合。另一方面，低阶模型似乎也不能恰当地拟合数据。这个叫做欠拟合 （underfitting）。

所以让我们公平地看待2阶或者更高阶的模型，并且试验一下如果只拟合最后一周数据的话，会有什么样的效果。毕竟，我们相信最后一周的数据比之前的数据更符合未来数据的趋势。下面这个有些迷幻的图表中给出了结果。这里更加明显地显示出过拟合问题是如何不好的。

然而，当模型只在3.5周及以后数据上训练时，从模型误差中判断，仍然应该选择最复杂的那个模型。

Error d=1: 22143941.107618
Error d=2: 19768846.989176
Error d=3: 19766452.361027
Error d=10: 18949339.348539
Error d=100: 16915159.603877

5. 训练与测试

如果有一些未来数据能用于模型评估，那么仅从近似误差结果中就应该可以判断出我们选择的模型是好是坏了。

尽管我们看不到未来的数据，但可以从现有数据中拿出一部分，来模拟类似的效果。例如，把一定比例的数据删掉，并使用剩下的数据进行训练，然后在拿出的那部分数据上计算误差。由于模型在训练中看不见拿出的那部分数据，所以就可以对模型的未来行为得到一个较为真实的预估。

只利用拐点时间后的数据训练出来的模型，其测试误差显现出了一个完全不同的境况。

Error d=1: 7,917,335.831122
Error d=2: 6,993,880.348870
Error d=3: 7,137,471.177363
Error d=10: 8,805,551.189738
Error d=100: 10,877,646.621984

结果显示在下图中：

看来最终的胜者已经一目了然。2阶模型的测试误差最低，而这个误差是在模型训练中未使用的那部分数据上评估得到的。这让我们相信，当未来数据到来时，不会遇到糟糕的意外。

6. 回答最初的问题

最终得到了一个模型，我们认为它可以最好地代表数据生成过程；现在，要获悉我们的基础设施何时到达每小时100 000次请求，已经是一个简单的事情了，只需要计算何时我们的模型函数到达100 000这个值即可。

对于2阶模型，我们可以简单地计算出它的逆函数，并得到100 000上的结果。当然，我们还希望有一个可以适用于任何模型函数的方法。

可以这样做：从多项式中减去100 000，得到另一个多项式，然后计算出它的根。如果提供了参数的初始值，SciPy的optimize 模块有一个fsolve 函数可以完成这项工作。假设这个胜出的2阶多项式是fbt2 ：

>>> print(fbt2)
2
0.08844 x - 97.31 x + 2.853e+04
>>> print(fbt2-100000)
2
0.08844 x - 97.31 x - 7.147e+04
>>> from scipy.optimize import fsolve
>>> reached_max = fsolve(fbt2-100000, 800)/(7*24)
>>> print("100,000 hits/hour expected at week %f" % reached_max[0])
100,000 hits/hour expected at week 9.827613

模型告诉我们，鉴于目前的用户行为和我们公司的推进力，还有一个月才会到达访问容量的界限。

当然，加入了我们的预测之后，会出现一定的不确定性。要获知真实的情况，需要更复杂的统计学知识来计算我们在望向更远处时所期望的方差。

对一些用户和潜在用户的动态行为，仍然无法准确地建模。但是，目前的预测对我们来说已经不错了。毕竟，现在可以对所有的耗时行为有所准备。如果能够对Web流量密切监控，我们就可以及时发现何时需要分配新的资源。