Python-数据分析秘籍-三-

Python 数据分析秘籍（三）

原文：Mastering Python Data Analysis

协议：CC BY-NC-SA 4.0

七、使用金融数据分析选股

在本章中，我们将介绍以下食谱:

• 计算简单和日志返回
• 用夏普比率和流动性对股票进行排名
• 用卡尔马尔比率和索蒂诺比率对股票进行排名
• 分析退货统计数据
• 将个股与大盘联系起来
• 探索风险和回报
• 用非参数游程检验检验市场
• 随机漫步测试
• 用自回归模型确定市场效率
• 为股票价格数据库创建表格
• 填充股票价格数据库
• 优化等权双资产投资组合

简介

金融学涉及许多学科，如货币、储蓄、投资和保险。在这一章中，我们将集中讨论股票投资，因为股价数据非常丰富。根据学术理论，一般投资者不应该投资个股，而应该投资整个市场，例如，代表一个国家内大公司的一篮子股票。经济学家为这一理论提出了几个这样的论点。首先，金融市场是随机的；因此，通过选股击败一个平均篮子是非常困难的。其次，个股波动剧烈，价格波动剧烈。这些价格变动在一个篮子里得到平均，这使得投资一组股票的风险降低。

我们将分析股票价格，但没有什么能阻止你重用食谱来分析共同基金和交易所交易基金或其他金融资产。为了简化分析，我将选择范围限制在六只知名美国公司的股票，这些公司也在标准普尔 500 股票指数中有代表。

计算简单和日志返回

收益衡量(股票)价格的变化率。使用收益的好处是收益是无量纲的，所以我们可以很容易地比较不同金融证券的收益。相比之下，仅金融资产的价格并不能告诉我们太多。在本章中，我们计算每日回报，因为我们的数据是每天采样的。通过小的调整，您应该能够在不同的时间范围内应用相同的分析。

事实上，有各种类型的回报。出于基本分析的目的，我们只需要了解 simple (7.1)和 log(arithmic)返回(7.2)，如下式所示:

实际上，这些类型的返回可以很容易地转换，从简单的返回到日志返回。如果给你选择，日志返回是你应该喜欢的，因为它们更容易计算。

怎么做...

1. 进口情况如下:

   import dautil as dl
   import ch7util
   import matplotlib.pyplot as plt
   

2. 下载标准普尔 500 指数数据:

   ohlc = dl.data.OHLC()
   sp500 = ohlc.get('^GSPC')['Adj Close']
   rets = sp500[1:]/sp500[:-1] - 1
   

3. 绘制简单和日志返回:

   _, ax = plt.subplots()
   cp = dl.plotting.CyclePlotter(ax)
   cp.plot(sp500.index, rets, label='Simple')
   cp.plot(sp500.index[1:], ch7util.log_rets(sp500), label='Log')
   ax.set_title('Simple and Log Returns')
   ax.set_xlabel('Date')
   ax.set_ylabel('Return')
   ax.legend(loc='best')
   

最终结果参考以下截图(简单和日志返回的值非常接近):

本食谱的代码在本书代码包的simple_log_rets.ipynb文件中。

另见

• 相关维基百科页面在en.wikipedia.org/wiki/Rate_o…(2015 年 10 月检索)

用夏普比率和流动性对股票进行排名

由 T2 威廉夏普定义的夏普比率是一个基本的投资指标。比例如下:

这个比率取决于资产的回报和基准的回报。我们将使用标准普尔 500 指数作为基准。这个比率应该代表一个报酬与风险的比率。我们希望在最小化风险的同时最大化回报，这与最大化夏普比率相对应。

另一个重要的投资变量是流动性。现金是最终的流动资产，但大多数其他资产的流动性较低，这意味着当我们试图出售或购买它们时，它们会改变价值。我们将使用本食谱中的交易量作为流动性的衡量标准。(交易量对应于金融资产的交易数量。流动性衡量一项资产的流动性——买卖它有多容易。)

怎么做...

你可以在本书的代码包sharpe_liquidity.ipynb文件中找到代码:

1. 进口情况如下:

   import numpy as np
   import dautil as dl
   import matplotlib.pyplot as plt
   import ch7util
   

2. 定义以下函数计算平均交易量的比值和对数:

   def calc_metrics(ticker, ohlc):
       stock = ohlc.get(ticker)
       sp500 = ohlc.get('^GSPC')
       merged = ch7util.merge_sp500(stock, sp500)
       rets_stock = ch7util.log_rets(merged['Adj Close_stock'])
       rets_sp500 = ch7util.log_rets(merged['Adj Close_sp500'])
       stock_sp500 = rets_stock - rets_sp500
       sharpe_stock = stock_sp500.mean()/stock_sp500.std()
       avg_vol = np.log(merged['Volume_stock'].mean())
   
       return (sharpe_stock, avg_vol)
   

3. 从ch7util模块计算我们一篮子股票的指标:

   dfb = dl.report.DFBuilder(cols=['Ticker', 'Sharpe', 'Log(Average Volume)'])
   
   ohlc = dl.data.OHLC()
   
   for symbol in ch7util.STOCKS:
       sharpe, vol = calc_metrics(symbol, ohlc)
       dfb.row([symbol, sharpe, vol])
   
   df = dfb.build(index=ch7util.STOCKS)
   

4. 绘制股票的比率和对数平均成交量:

   _, ax = plt.subplots()
   ax.scatter(df['Log(Average Volume)'], df['Sharpe'])
   dl.plotting.plot_polyfit(ax, df['Log(Average Volume)'], df['Sharpe'])
   
   dl.plotting.plot_text(ax, df['Log(Average Volume)'],
                         df['Sharpe'], ch7util.STOCKS)
   ax.set_xlabel('Log(Average Volume)')
   ax.set_ylabel('Sharpe')
   ax.set_title('Sharpe Ratio & Liquidity')
   

有关最终结果，请参考以下屏幕截图:

另见

• 相关维基百科页面在en.wikipedia.org/wiki/Sharpe…(2015 年 10 月检索)

用卡尔马尔比率和索蒂诺比率对股票进行排名

Sortino 和 Calmar 比率是与夏普比率相当的业绩比率(参考用夏普比率和流动性配方对股票进行排名)。比率甚至更多；然而，夏普比率一直在附近最长，因此应用非常广泛。

Sortino 比率以 Frank Sortino 命名，但它是由 Brian Rom 定义的。比率将风险定义为低于基准的下行方差。基准可以是指数，也可以是固定收益，如零。该比率定义如下:

R 为资产回报， T 为目标基准， DR 为下行风险。卡尔马尔比率是由特里·杨发明的，并以他的公司和时事通讯命名。该比率将风险定义为资产的最大提取(价格从峰值跌至底部)。

怎么做...

以下是本书代码包中calmar_sortino.ipynb文件的分解:

1. 进口情况如下:

   import numpy as np
   import dautil as dl
   import ch7util
   from scipy.signal import argrelmin
   from scipy.signal import argrelmax
   import matplotlib.pyplot as plt
   

2. 定义以下函数计算排序比:

   def calc_sortino(rets):
       # Returns below target
       semi_var = rets[rets < 0] ** 2
       semi_var = semi_var.sum()/len(rets)
       sortino = np.sqrt(semi_var)
   
       return rets.mean()/sortino
   

3. 定义以下函数计算卡尔马尔比:

   def calc_calmar(rets):
       # Peaks and bottoms indexes in sequence
       mins = np.ravel(argrelmin(rets))
       maxs = np.ravel(argrelmax(rets))
       extrema = np.concatenate((mins, maxs))
       extrema.sort()
   
       return -rets.mean()/np.diff(rets[extrema]).min()
   

4. 为我们的股票列表计算卡尔马尔和索蒂诺比率:

   ohlc = dl.data.OHLC()
   dfb = dl.report.DFBuilder(cols=['Ticker', 'Sortino', 'Calmar'])
   
   for symbol in ch7util.STOCKS:
       stock = ohlc.get(symbol)
       rets = ch7util.log_rets(stock['Adj Close'])
       sortino = calc_sortino(rets)
       calmar = calc_calmar(rets)
       dfb.row([symbol, sortino, calmar])
   
   df = dfb.build(index=ch7util.STOCKS).dropna()
   

5. 绘制股票的 Sortino 和 Calmar 比率:

   _, ax = plt.subplots()
   ax.scatter(df['Sortino'], df['Calmar'])
   dl.plotting.plot_polyfit(ax, df['Sortino'], df['Calmar'])
   dl.plotting.plot_text(ax, df['Sortino'], df['Calmar'], ch7util.STOCKS)
   ax.set_xlabel('Sortino')
   ax.set_ylabel('Calmar')
   ax.set_title('Sortino & Calmar Ratios')
   

有关最终结果，请参考以下屏幕截图:

另见

• 维基百科关于 en.wikipedia.org/wiki/Sortin… 索特诺比率的页面(2015 年 10 月检索)
• 维基百科关于 en.wikipedia.org/wiki/Calmar… 卡尔马尔比率的页面(2015 年 10 月检索)

分析退货统计

回报，尤其是股票指数的回报，已经被广泛研究。过去，人们认为回报是正态分布的。然而，现在很明显，回报分布有肥尾(比正态分布更肥)。更多信息可在en.wikipedia.org/wiki/Fat-ta…获得(2015 年 10 月检索)。检查数据是否符合正态分布非常容易。我们只需要样本的平均值和标准差。

在本食谱中，我们将探讨许多主题:

• 股票收益的偏斜度和峰度是一个值得研究的问题。偏斜度在股票期权模型中尤其重要。分析师通常将自己限制在均值和标准差范围内，假设它们分别对应于报酬和风险。
• 如果我们对趋势的存在感兴趣，那么我们应该看一看自相关图。这是一个自相关图——即信号和某个滞后的信号之间的相关性(也在 Python 数据分析、 Packt Publishing 中解释)。
• 我们还将绘制负回报(其绝对值)和对数-对数标度上的相应计数，因为这些似乎近似遵循幂律(尤其是尾值)。

怎么做...

分析可以在本书代码包的rets_stats.ipynb文件中找到:

1. 进口情况如下:

   import dautil as dl
   import ch7util
   import matplotlib.pyplot as plt
   from scipy.stats import skew
   from scipy.stats import kurtosis
   from pandas.tools.plotting import autocorrelation_plot
   import numpy as np
   from scipy.stats import norm
   from IPython.display import HTML
   

2. 计算我们股票的回报:

   ohlc = dl.data.OHLC()
   rets_dict = {}
   
   for i, symbol in enumerate(ch7util.STOCKS):
       rets = ch7util.log_rets(ohlc.get(symbol)['Adj Close'])
       rets_dict[symbol] = rets
   
   sp500 = ch7util.log_rets(ohlc.get('^GSPC')['Adj Close'])
   

3. 绘制标准普尔 500 回报和相应理论正态分布的直方图:

   sp = dl.plotting.Subplotter(2, 2, context)
   sp.ax.set_xlim(-0.05, 0.05)
   _, bins, _ = sp.ax.hist(sp500, bins=dl.stats.sqrt_bins(sp500),
                        alpha=0.6, normed=True)
   sp.ax.plot(bins, norm.pdf(bins, sp500.mean(), sp500.std()), lw=2)
   

4. 绘制收益的偏斜度和峰度:

   skews = [skew(rets_dict[s]) for s in ch7util.STOCKS]
   kurts = [kurtosis(rets_dict[s]) for s in ch7util.STOCKS]
   sp.label()
   
   sp.next_ax().scatter(skews, kurts)
   dl.plotting.plot_text(sp.ax, skews, kurts, ch7util.STOCKS)
   sp.label()
   

5. 绘制标准普尔 500 回报的自相关图:

   autocorrelation_plot(sp500, ax=sp.next_ax())
   sp.label()
   

6. 绘制负回报(绝对值)和计数的对数-对数图:

   # Negative returns
   counts, neg_rets = np.histogram(sp500[sp500 < 0])
   neg_rets = neg_rets[:-1] + (neg_rets[1] - neg_rets[0])/2
   # Adding 1 to avoid log(0)
   dl.plotting.plot_polyfit(sp.next_ax(), np.log(np.abs(neg_rets)),
                            np.log(counts + 1), plot_points=True)
   sp.label()
   
   HTML(sp.exit())
   

有关最终结果，请参考以下屏幕截图:

将个股与大盘关联起来

当我们定义一个股市或指数时，我们通常会选择在某些方面相似的股票。例如，股票可能在同一个国家或大陆。鸟类的位置可以从它们所属的鸟群的位置来粗略估计。同样，我们预计股票回报将与其市场相关，尽管不一定完全相关。

我们将探讨以下指标:

• 最明显的指标据称是个股回报率和标准普尔 500 指数的相关系数。
• 另一个指标是从线性回归而不是相关性中获得的斜率。
• 我们也可以分析收益的平方差，有点类似于回归诊断中的平方误差。
• 我们也可以将交易量和波动性联系起来，而不是将回报联系起来。为了衡量波动性，我们将使用有点不寻常的高低价格差的平方值。实际上，我们应该把这个值除以一个常数；然而，这对于相关系数计算来说不是必需的。

怎么做...

分析在本书代码包的correlating_market.ipynb文件中:

1. 进口情况如下:

   import ch7util
   import dautil as dl
   import numpy as np
   import matplotlib.pyplot as plt
   from IPython.display import HTML
   

2. 定义以下函数来计算波动率:

   def hl2(df, suffix):
       high = df['High_' + suffix]
       low = df['Low_' + suffix]
   
       return (high - low) ** 2
   

3. 定义以下函数来关联标准普尔 500 和我们的股票:

   def correlate(stock, sp500):
       merged = ch7util.merge_sp500(stock, sp500)
       rets = ch7util.log_rets(merged['Adj Close_stock'])
       sp500_rets = ch7util.log_rets(merged['Adj Close_sp500'])
       result = {}
   
       result['corrcoef'] = np.corrcoef(rets, sp500_rets)[0][1]
       slope, _ = np.polyfit(sp500_rets, rets, 1)
       result['slope'] = slope
   
       srd = (sp500_rets - rets) ** 2
       result['msrd'] = srd.mean()
       result['std_srd'] = srd.std()
   
       result['vols'] = np.corrcoef(merged['Volume_stock'],
                                    merged['Volume_sp500'])[0][1]
   
       result['hl2'] = np.corrcoef(hl2(merged, 'stock'),
                                   hl2(merged, 'sp500'))[0][1]
   
       return result
   

4. 将我们的一组股票与标准普尔 500 指数联系起来

5. 绘制股票的相关系数:

   sp = dl.plotting.Subplotter(2, 2, context)
   dl.plotting.bar(sp.ax, ch7util.STOCKS,
                   [corr['corrcoef'] for corr in corrs])
   sp.label()
   
   dl.plotting.bar(sp.next_ax(), ch7util.STOCKS,
                   [corr['slope'] for corr in corrs])
   sp.label()
   

6. 绘制差异统计的平方:

   sp.next_ax().set_xlim([0, 0.001])
   dl.plotting.plot_text(sp.ax, [corr['msrd'] for corr in corrs],
                         [corr['std_srd'] for corr in corrs],
                         ch7util.STOCKS, add_scatter=True,
                         fontsize=9, alpha=0.6)
   sp.label()
   

7. 绘制成交量和波动率相关系数:

   dl.plotting.plot_text(sp.next_ax(), [corr['vols'] for corr in corrs],
                         [corr['hl2'] for corr in corrs],
                         ch7util.STOCKS, add_scatter=True,
                         fontsize=9, alpha=0.6)
   sp.label()
   
   HTML(sp.exit())
   

最终结果参见以下截图:

探索风险与回报

金融中的贝塔系数是一个线性回归模型的斜率(T2 ),该模型包含了 T3 资产的收益和一个基准的收益，例如 T4 500 指数。该模型定义如下:

根据 资本资产定价模型 ( 资本资产定价模型)的说法，贝塔是衡量风险的指标。预期收益由收益的平均值给出。如果我们绘制各种证券的贝塔值和预期收益，我们就获得了相应市场的证券市场线 ( SML )。SML 的截距给出了无风险率，理论上我们应该在不承担任何风险的情况下获得回报。一般来说，如果一项资产不在 SML，那么根据资本资产定价模型，它被错误定价。

怎么做...

程序在本书代码包的capm.ipynb文件中:

1. 进口情况如下:

   import dautil as dl
   import numpy as np
   import pandas as pd
   import ch7util
   import matplotlib.pyplot as plt
   

2. 定义以下函数来计算β值:

   def calc_beta(symbol):
       ohlc = dl.data.OHLC()
       sp500 = ohlc.get('^GSPC')['Adj Close']
       stock = ohlc.get(symbol)['Adj Close']
       df = pd.DataFrame({'SP500': sp500, symbol: stock}).dropna()
       sp500_rets = ch7util.log_rets(df['SP500'])
       rets = ch7util.log_rets(df[symbol])
       beta, _ = np.polyfit(sp500_rets, rets, 1)
   
       # annualize & percentify
       return beta, 252 * rets.mean() * 100
   

3. 计算我们股票的贝塔值和平均回报:

   betas = []
   means = []
   
   for symbol in ch7util.STOCKS:
       beta, ret_mean = calc_beta(symbol)
       betas.append(beta)
       means.append(ret_mean)
   

4. 绘制结果和市场安全线:

   _, ax = plt.subplots()
   dl.plotting.plot_text(ax, betas, means, ch7util.STOCKS, add_scatter=True)
   dl.plotting.plot_polyfit(ax, betas, means)
   ax.set_title('Capital Asset Pricing Model')
   ax.set_xlabel('Beta')
   ax.set_ylabel('Mean annual return (%)')
   

有关最终结果，请参考以下屏幕截图:

另见

• 关于 en.wikipedia.org/wiki/Beta_%… 测试版的维基百科页面(2015 年 10 月检索)
• 维基百科关于 en.wikipedia.org/wiki/Capita… 资本资产定价模型的页面(2015 年 10 月检索)

用非参数游程检验检验市场

有效市场假说(【EMH】)规定你平均不能通过挑选更好的股票或把握市场时机来“跑赢市场”。根据《EMH》的说法，所有关于市场的信息都可以以这样或那样的形式立即提供给每个市场参与者，并立即反映在资产价格中，因此投资就像是在玩一场纸牌游戏，所有的牌都露出来了。你能赢的唯一方法是在风险很大的股票上下注，然后走运。

法国数学家学士在 1900 年左右为 EMH 开发了一个测试。该测试考察了连续出现的正负价格变化。我们不计算价格没有变化的事件，只使用它们来结束运行。无论如何，对于流动性市场来说，这类事件相对较少。

统计测试本身在金融界之外为人所知，被称为“T2”沃尔德-沃尔福威茨运行测试“T3”。如果我们用“+”表示正变化，用“-”表示负变化，我们可以用 5 次运行得到序列“++++++++++++”。以下计算运行次数的平均值μ (7.6)、标准偏差σ (7.7)和 Z 评分 Z (7.8)的公式 R 还要求负变化次数 N- ，正变化次数 N+ ，以及总变化次数 N :

我们假设运行次数遵循正态分布，这给了我们一种方法，在我们选择的置信水平下，潜在地拒绝运行的随机性。

怎么做...

看看这本书代码包里的non_parametric.ipynb文件。

1. 进口情况如下:

   import dautil as dl
   import numpy as np
   import pandas as pd
   import ch7util
   import matplotlib.pyplot as plt
   from scipy.stats import norm
   from IPython.display import HTML
   

2. 定义以下函数来计算运行次数:

   def count_runs(signs):
       nruns = 0
       prev = None
   
       for s in signs:
           if s != 0 and s != prev:
               nruns += 1
   
           prev = s
   
       return nruns
   

3. 定义以下函数来计算平均值、标准偏差和 z 分数:

   def proc_runs(symbol):
       ohlc = dl.data.OHLC()
       close = ohlc.get(symbol)['Adj Close'].values
       diffs = np.diff(close)
       nplus = (diffs > 0).sum()
       nmin = (diffs < 0).sum()
       n = nplus + nmin
       mean = (2 * (nplus * nmin) / n) + 1
       var = (mean - 1) * (mean - 2) / (n - 1)
       std = np.sqrt(var)
       signs = np.sign(diffs)
       nruns = count_runs(np.diff(signs))
   
       return mean, std, (nruns - mean) / std
   

4. 计算我们股票的指标:

   means = []
   stds = []
   zscores = []
   
   for symbol in ch7util.STOCKS:
       mean, std, zscore = proc_runs(symbol)
       means.append(mean)
       stds.append(std)
       zscores.append(zscore)
   

5. 用表示 95%置信水平的线绘制 z 分数:

   sp = dl.plotting.Subplotter(2, 1, context)
   dl.plotting.plot_text(sp.ax, means, stds, ch7util.STOCKS, add_scatter=True)
   sp.label()
   
   dl.plotting.bar(sp.next_ax(), ch7util.STOCKS, zscores)
   sp.ax.axhline(norm.ppf(0.95), label='95 % confidence level')
   sp.label()
   HTML(sp.exit())
   

最终结果参见以下截图:

另见

• 维基百科关于沃尔德-沃尔福威茨运行测试的页面位于(2015 年 10 月检索)
• 关于 en.wikipedia.org/wiki/Effici… EMH 的维基百科页面(2015 年 10 月检索)

随机行走测试

随机游走假说 ( RWH )就像有效市场假说(参考用非参数游程检验食谱检验市场)一样，声称市场是打不过的。 RWH 规定，资产价格执行随机游走。事实上，只要反复抛硬币，你就能生成相当令人信服的股价图。

1988 年，金融学教授罗和麦克莱恩斯用资产价格的自然对数作为数据为构建了一个测试。该测试指定原木价格在平均值(7.9)附近漂移。我们预计不同频率(例如，一天和两天)的价格变化是随机的。此外，两个不同频率下的方差(7.10 和 7.11)是相关的，根据以下等式，相应的比值(7.12)通常分布在零附近:

怎么做...

代码在本书代码包的random_walk.ipynb文件中:

1. 进口情况如下:

   import dautil as dl
   import numpy as np
   import matplotlib.pyplot as plt
   import ch7util
   

2. 计算我们股票的比率:

   ratios = []
   
   for symbol in ch7util.STOCKS:
       ohlc = dl.data.OHLC()
       P = ohlc.get(symbol)['Adj Close'].values
       N = len(P)
       mu = (np.log(P[-1]) - np.log(P[0]))/N
       var_a = 0
       var_b = 0
   
       for k in range(1, N):
           var_a = (np.log(P[k]) - np.log(P[k - 1]) - mu) ** 2
           var_a = var_a / N
   
       for k in range(1, N//2):
           var_b = (np.log(P[2 * k]) - np.log(P[2 * k - 2]) - 2 * mu) ** 2
           var_b = var_b / N
   
       ratios.append(np.sqrt(N) * (var_b/var_a - 1))
   

3. 绘制比率图，我们预计比率接近于零(7.12):

   _, ax = plt.subplots()
   dl.plotting.bar(ax, ch7util.STOCKS, ratios)
   ax.set_title('Random Walk Test')
   ax.set_ylabel('Ratio')
   

有关最终结果，请参考以下屏幕截图:

另见

• 维基百科页面关于 en.wikipedia.org/wiki/Random… 的随机行走假说(2015 年 10 月检索)
• 菲南德牧师。种马。(1988) 1 (1): 41-66.doi:10 . 1093/rfs/1 . 1 . 41rfs.oxfordjournals.org/content/1/1…(2015 年 10 月检索)

用自回归模型确定市场效率

根据有效市场假说(参考用非参数游程检验方法检验市场)，关于资产的所有信息立即反映在资产的价格中。这意味着以前的价格不会影响现在的价格。以下方程指定了自回归模型(7。13)和受限模型(7。14)将所有系数设置为零:

如果我们相信市场是有效的，我们会期望非限制模型比限制模型没有什么增加，因此，比率(7。15)的相应 R 平方系数应该接近 1。

怎么做...

脚本在本书代码包的autoregressive_test.ipynb文件中:

1. 进口情况如下:

   import dautil as dl
   import ch7util
   import numpy as np
   import matplotlib.pyplot as plt
   import statsmodels.api as sm
   from IPython.display import HTML
   

2. 使用(7.13)和(7.14)拟合模型，然后使用(7.15)为我们的股票列表计算市场效率:

   ohlc = dl.data.OHLC()
   efficiencies = []
   restricted_r2 = []
   unrestricted_r2 = []
   
   for stock in ch7util.STOCKS:
       rets = ch7util.log_rets(ohlc.get(stock)['Adj Close'])
       restricted = sm.OLS(rets, rets.mean() * np.ones_like(rets)).fit()
       rets_1 = rets[3:-1]
       rets_2 = rets[2:-2]
       rets_3 = rets[1:-3]
       rets_4 = rets[:-4]
       x = np.vstack((rets_1, rets_2, rets_3, rets_4)).T
       x = sm.add_constant(x)
       y = rets[4:]
       unrestricted = sm.OLS(y, x).fit()
       restricted_r2.append(restricted.rsquared)
       unrestricted_r2.append(unrestricted.rsquared)
       efficiencies.append(1 - restricted.rsquared/unrestricted.rsquared)
   

3. 将市场效率和 R 平方值绘制如下:

   sp = dl.plotting.Subplotter(2, 1, context)
   dl.plotting.bar(sp.ax, ch7util.STOCKS, efficiencies)
   sp.label()
   dl.plotting.plot_text(sp.next_ax(), unrestricted_r2, np.array(restricted_r2)/10 ** -16,
                         ch7util.STOCKS, add_scatter=True)
   sp.label()
   HTML(sp.exit())
   

有关最终结果，请参考以下屏幕截图:

另见

• en.wikipedia.org/wiki/Autore… T2 的自回归模型的维基百科页面(2015 年 10 月检索)

为股票价格数据库创建表格

只存储股票价格一般来说不是很有用。我们通常希望存储关于公司和相关衍生品的额外静态信息，如股票期权和期货。经济学理论告诉我们，在历史价格数据中寻找周期和趋势或多或少是浪费时间；因此，创建数据库似乎更没有意义。当然你不必相信这个理论，无论如何创建一个股票价格数据库是一个有趣的技术挑战。数据库对于投资组合优化也很有用(参见食谱优化同等权重的 2 资产投资组合)。

我们将基于中的星型模式进行设计，用事实和维度表实现星型模式。事实表将保存价格，包括日期维度表、资产维度表和来源维度表，如下图所示:

显然，模式会随着时间的推移而发展，根据需要添加或删除表、索引和列。我们将使用中的模式填充股票价格数据库配方。

怎么做...

该模式在本书的代码包中的database_tables.py文件中定义:

1. 进口情况如下:

   from sqlalchemy import Column
   from sqlalchemy.ext.declarative import declarative_base
   from sqlalchemy import Date
   from sqlalchemy import ForeignKey
   from sqlalchemy import Integer
   from sqlalchemy import String
   
   Base = declarative_base()
   

2. 为股价事实表定义以下类别:

   class StockPrice(Base):
       __tablename__ = 'stock_price'
       id = Column(Integer, primary_key=True)
       date_id = Column(Integer, ForeignKey('date_dim.id'),
                        primary_key=True)
       asset_id = Column(Integer, ForeignKey('asset_dim.id'),
                         primary_key=True)
       source_id = Column(Integer, ForeignKey('source_dim.id'),
                          primary_key=True)
       open_price = Column(Integer)
       high_price = Column(Integer)
       low_price = Column(Integer)
       close_price = Column(Integer)
       adjusted_close = Column(Integer)
       volume = Column(Integer)
   

3. 为日期维度表定义以下类:

   class DateDim(Base):
       __tablename__ = 'date_dim'
       id = Column(Integer, primary_key=True)
       date = Column(Date, nullable=False, unique=True)
       day_of_month = Column(Integer, nullable=False)
       day_of_week = Column(Integer, nullable=False)
       month = Column(Integer, nullable=False)
       quarter = Column(Integer, nullable=False)
       year = Column(Integer, nullable=False)
   

4. 定义以下类来保存股票信息:

   class AssetDim(Base):
       __tablename__ = 'asset_dim'
       id = Column(Integer, primary_key=True)
       symbol = Column(String, nullable=False, unique=True)
       name = Column(String, nullable=False)
       # Could make this a reference to separate table
       category = Column(String, nullable=False)
       country = Column(String, nullable=False)
       # Could make this a reference to separate table
       sector = Column(String, nullable=False)
   

5. 为源维度表定义以下类(雅虎财经只需要一个条目):

   class SourceDim(Base):
       __tablename__ = 'source_dim'
       id = Column(Integer, primary_key=True)
       name = Column(String, nullable=False)
       url = Column(String)
   

填充股票价格数据库

在为股票价格数据库配方创建表格中，我们为历史股票价格数据库定义了一个模式。在这个食谱中，我们将使用雅虎财经的数据填充表格，并绘制不同时间段和业务部门的平均交易量。

股市研究人员发现了几个与季节效应有关的奇怪现象。此外，还有某些反复出现的事件，如盈利公告、股息支付和期权到期。同样，经济理论告诉我们，我们观察到的任何模式要么是幻觉，要么是所有市场参与者都已经知道的。这是真是假很难证实；然而，这种方法作为数据分析的练习是很棒的。此外，您可以使用数据库来优化您的投资组合，如优化等权重 2 资产投资组合食谱中所述。

怎么做...

代码在本书代码包的populate_database.ipynb文件中:

1. 进口情况如下:

   import database_tables as tables
   import pandas as pd
   import os
   import dautil as dl
   import ch7util
   import sqlite3
   import matplotlib.pyplot as plt
   import seaborn as sns
   from IPython.display import HTML
   

2. 定义以下函数来填充日期维度表:

   def populate_date_dim(session):
       for d in pd.date_range(start='19000101', end='20250101'):
           adate = tables.DateDim(date=d.date(), day_of_month=d.day,
                                  day_of_week=d.dayofweek, month=d.month,
                                  quarter=d.quarter, year=d.year)
           session.add(adate)
   
       session.commit()
   

3. 定义以下功能填充资产维度表:

   def populate_asset_dim(session):
       asset = tables.AssetDim(symbol='AAPL', name='Apple Inc.',
                               category='Common Stock', country='USA',
                               sector='Consumer Goods')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='INTC', name='Intel Corporation',
                               category='Common Stock', country='USA',
                               sector='Technology')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='MSFT', name='Microsoft Corporation',
                               category='Common Stock', country='USA',
                               sector='Technology')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='KO', name='The Coca-Cola Company',
                               category='Common Stock', country='USA',
                               sector='Consumer Goods')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='DIS', name='The Walt Disney Company',
                               category='Common Stock', country='USA',
                               sector='Services')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='MCD', name='McDonald\'s Corp.',
                               category='Common Stock', country='USA',
                               sector='Services')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='NKE', name='NIKE, Inc.',
                               category='Common Stock', country='USA',
                               sector='Consumer Goods')
       session.add(asset)
   
       asset = tables.AssetDim(symbol='IBM',
                               name='International Business Machines Corporation',
                               category='Common Stock', country='USA',
                               sector='Technology')
       session.add(asset)
   
       session.commit()
   

4. 定义以下函数来填充源维度表:

   def populate_source_dim(session):
       session.add(tables.SourceDim(name='Yahoo Finance',
                                    url='https://finance.yahoo.com'))
       session.commit()
   

5. 定义以下函数来填充持有股票价格的事实表:

   def populate_prices(session):
       symbols = dl.db.map_to_id(session, tables.AssetDim.symbol)
       dates = dl.db.map_to_id(session, tables.DateDim.date)
       source_id = session.query(tables.SourceDim).first().id
       ohlc = dl.data.OHLC()
       conn = sqlite3.connect(dbname)
       c = conn.cursor()
       insert = '''INSERT INTO stock_price (id, date_id,
           asset_id, source_id, open_price, high_price, low_price,
           close_price, adjusted_close, volume)  VALUES({id}, {date_id},
           {asset_id}, {source_id}, {open_price}, {high_price},
           {low_price}, {close_price}, {adj_close}, {volume})'''
       logger = dl.log_api.conf_logger(__name__)
   
       for symbol in ch7util.STOCKS:
           df = ohlc.get(symbol)
           i = 0
   
           for index, row in df.iterrows():
               date_id = dates[index.date()]
               asset_id = symbols[symbol]
               i += 1
               stmt = insert.format(id=i, date_id=date_id,
                                    asset_id=asset_id,
                                    source_id=source_id,
                                    open_price=dl.data.centify(row['Open']),
                                    high_price=dl.data.centify(row['High']),
                                    low_price=dl.data.centify(row['Low']),
                                    close_price=dl.data.centify(row['Close']),
                                    adj_close=dl.data.centify(row['Adj Close']),
                                    volume=int(row['Volume']))
               c.execute(stmt)
   
               if i % 1000 == 0:
                   logger.info("Progress %s %s", symbol, i)
   
               conn.commit()
   
           conn.commit()
   
       c.close()
       conn.close()
   

6. 定义以下函数来填充所有表:

   def populate(session):
       if session.query(tables.SourceDim).count() == 0:
           populate_source_dim(session)
           populate_asset_dim(session)
           populate_date_dim(session)
           populate_prices(session)
   

7. 定义以下函数来绘制平均体积:

   def plot_volume(col, ax):
       df = pd.read_sql(sql.format(col=col), conn)
       sns.barplot(x=col, y='AVG(P.Volume/1000)', data=df,
                   hue='sector', ax=ax)
   
       ax.legend(loc='best')
   
   dbname = os.path.join(dl.data.get_data_dir(), 'stock_prices.db')
   session = dl.db.create_session(dbname, tables.Base)
   populate(session)
   sql = '''
       SELECT
           A.sector,
           D.{col},
           AVG(P.Volume/1000)
       FROM stock_price P
       INNER JOIN date_dim D  ON (P.Date_Id = D.Id)
       INNER JOIN asset_dim A ON (P.asset_id = a.Id)
       GROUP BY
           A.sector,
           D.{col}
         '''
   

8. 用以下代码绘制平均体积:

   conn = sqlite3.connect(dbname)
   
   sp = dl.plotting.Subplotter(2, 2, context)
   plot_volume('day_of_week', sp.ax)
   sp.ax.set_xticklabels(['Mon', 'Tue', 'Wed', 'Thu', 'Fri'])
   
   plot_volume('month', sp.next_ax())
   sp.ax.set_xticklabels(dl.ts.short_months())
   
   plot_volume('day_of_month', sp.next_ax())
   plot_volume('quarter', sp.next_ax())
   HTML(sp.exit())
   

有关最终结果，请参考以下屏幕截图:

优化等权重双资产组合

买卖股票有点像购物。购物是超市和网上书店都很熟悉的事情。这些类型的业务通常应用技术，如篮子分析和推荐引擎。例如，如果你是一个写历史不准确小说的作家的粉丝，推荐引擎可能会推荐同一作家的另一部小说或其他历史不准确的小说。

股票推荐引擎不能这样工作。例如，如果你的投资组合中只有石油生产商的股票，而油价对你不利，那么整个投资组合将失去价值。所以，我们应该尝试拥有来自不同行业、行业或地理区域的股票。我们可以用回报的相关性来衡量相似性。

类似于夏普比率(参考用夏普比率和流动性排名股票食谱)，我们希望最大化我们的投资组合的平均回报，最小化投资组合回报的方差。这些想法也存在于现代投资组合理论 ( MPT )中，其发明者被授予诺贝尔奖。对于双资产投资组合，我们有以下等式:

权重 wA 和 wB 为组合权重，合计为 1。权重可以是负的——因为投资者可以卖空(不持有就卖出，这会产生借贷成本)一种证券。我们可以用线性代数方法或一般优化算法解决投资组合优化问题。然而，对于权重相等、只有少数股票的双资产投资组合来说，蛮力方法已经足够好了。

怎么做...

以下是本书代码包中portfolio_optimization.ipynb文件的分解:

1. 进口情况如下:

   import dautil as dl
   import ch7util
   import pandas as pd
   import numpy as np
   import seaborn as sns
   import matplotlib.pyplot as plt
   

2. 定义以下函数来计算预期收益(7.16):

   def expected_return(stocka, stockb, means):
       return 0.5 * (means[stocka] + means[stockb])
   

3. 定义以下函数计算投资组合回报的方差(7.17):

   def variance_return(stocka, stockb, stds):
       ohlc = dl.data.OHLC()
       dfa = ohlc.get(stocka)
       dfb = ohlc.get(stockb)
       merged = pd.merge(left=dfa, right=dfb,
                         right_index=True, left_index=True,
                         suffixes=('_A', '_B')).dropna()
       retsa = ch7util.log_rets(merged['Adj Close_A'])
       retsb = ch7util.log_rets(merged['Adj Close_B'])
       corr = np.corrcoef(retsa, retsb)[0][1]
   
       return 0.25 * (stds[stocka] ** 2 + stds[stockb] ** 2 +
                      2 * stds[stocka] * stds[stockb] * corr)
   

4. 定义以下函数计算预期收益与方差的比值:

   def calc_ratio(stocka, stockb, means, stds, ratios):
       if stocka == stockb:
           return np.nan
   
       key = stocka + '_' + stockb
       ratio = ratios.get(key, None)
   
       if ratio:
           return ratio
   
       expected = expected_return(stocka, stockb, means)
       var = variance_return(stocka, stockb, stds)
       ratio = expected/var
       ratios[key] = ratio
   
       return ratio
   

5. 计算每只股票的平均回报和标准差:

   means = {}
   stds = {}
   
   ohlc = dl.data.OHLC()
   
   for stock in ch7util.STOCKS:
       close = ohlc.get(stock)['Adj Close']
       rets = ch7util.log_rets(close)
       means[stock] = rets.mean()
       stds[stock] = rets.std()
   

6. 计算我们所有股票组合的网格比率:

   pairs = dl.collect.grid_list(ch7util.STOCKS)
   sorted_pairs = [[sorted(row[i]) for row in pairs]
                   for i in range(len(ch7util.STOCKS))]
   ratios = {}
   
   grid = [[calc_ratio(row[i][0], row[i][1], means, stds, ratios)
           for row in sorted_pairs] for i in range(len(ch7util.STOCKS))]
   

7. 在热图中绘制网格如下:

   %matplotlib inline
   plt.title('Expected Return/Return Variance for 2 Asset Portfolio')
   sns.heatmap(grid, xticklabels=ch7util.STOCKS, yticklabels=ch7util.STOCKS)
   

最终结果参见以下截图:

另见

• 关于 en.wikipedia.org/wiki/Modern… MPT 的维基百科页面(2015 年 10 月检索)

八、文本挖掘与社交网络分析

在本章中，我们将介绍以下食谱:

• 创建分类语料库
• 用句子和单词标记新闻文章
• 词干，引理，过滤和 TF-IDF 分数
• 识别命名实体
• 利用非负矩阵分解提取主题
• 实现基本术语数据库
• 计算社交网络密度
• 计算社交网络接近度中心度
• 确定中间中心性
• 估计平均聚类系数
• 计算图的分类系数
• 得到图的团数
• 创建具有余弦相似性的文档图

简介

人类通过语言交流已经有几千年了。手写文本已经存在了很长时间，古腾堡出版社当然是一个巨大的发展，但是现在我们有了计算机、互联网和社交媒体，事情肯定已经失控了。

本章将帮助你应对文本和社交媒体信息的泛滥。我们将使用的主要 Python 库是 NLTK 和 NetworkX。你必须真正理解在这些库中可以找到多少特性。用pip或conda安装 NLTK，如下所示:

$ conda/pip install nltk 

代码用 NLTK 3.0.2 进行了测试。如需下载语料库，请按照www.nltk.org/data.html给出的说明进行(2015 年 11 月检索)。

使用 pip或conda安装网络，如下所示:

$ conda/pip install networkx 

该代码用网络 1.9.1 进行了测试。

创建分类语料库

作为 Python 爱好者，我们对关于 Python 编程或相关技术的新闻感兴趣；但是，如果搜索 Python 文章，也可能会得到关于蛇的文章。这个问题的一个解决方案是训练一个识别相关文章的分类器。这需要一个训练集——一个分类的语料库，例如，“Python 编程”和“其他”类别。

NLTK 有CategorizedPlaintextCorpusReader类，用于构建分类语料库。为了让事情变得更加令人兴奋，我们将从 RSS 源中获取新闻文章的链接。我选择了英国广播公司的节目，但是你当然可以使用任何其他的节目。英国广播公司的节目已经分类了。我选择了世界新闻和科技新闻，所以这给了我们两个类别。提要不包含文章的全文，因此我们需要使用 Selenium 进行一些抓取，如第 5 章、网络挖掘、数据库和大数据所述。你可能需要对文本文件进行后处理，因为英国广播公司的网页不仅包含新闻故事的文本，还包含副刊。

做好准备

按照本章简介部分的说明安装 NLTK。安装用于处理 RSS 源的 feedparser:

$ pip/conda install feedparser

我用 feedparser 5.2.1 测试了这个配方。

怎么做...

1. 进口情况如下:

   import feedparser as fp
   import urllib
   from selenium import webdriver
   from selenium.webdriver.support.ui import WebDriverWait
   from selenium.webdriver.support import expected_conditions as EC
   from selenium.webdriver.common.by import By
   import dautil as dl
   from nltk.corpus.reader import CategorizedPlaintextCorpusReader
   import os
   

2. 创建以下变量来帮助刮擦:

   DRIVER = webdriver.PhantomJS()
   NAP_SECONDS = 10
   LOGGER = dl.log_api.conf_logger('corpus')
   

3. 定义以下功能存储文本内容:

   def store_txt(url, fname, title):
       try:
           DRIVER.get(url)
   
           elems = WebDriverWait(DRIVER, NAP_SECONDS).until(
               EC.presence_of_all_elements_located((By.XPATH, '//p'))
           )
   
           with open(fname, 'w') as txt_file:
               txt_file.write(title + '\n\n')
               lines = [e.text for e in elems]
               txt_file.write(' \n'.join(lines))
       except Exception:
           LOGGER.error("Error processing HTML", exc_info=True)
   

4. 定义以下函数来检索故事:

   def fetch_news(dir):
       base = 'http://newsrss.bbc.co.uk/rss/newsonline_uk_edition/{}/rss.xml'
   
       for category in ['world', 'technology']:
           rss = fp.parse(base.format(category))
   
           for i, entry in enumerate(rss.entries):
               fname = '{0}_bbc_{1}.txt'.format(i, category)
               fname = os.path.join(dir, fname)
   
               if not dl.conf.file_exists(fname):
                   store_txt(entry.link, fname, entry.title)
   

5. 用以下代码调用函数:

   if __name__ == "__main__":
       dir = os.path.join(dl.data.get_data_dir(), 'bbc_news_corpus')
   
       if not os.path.exists(dir):
           os.mkdir(dir)
   
       fetch_news(dir)
       reader = CategorizedPlaintextCorpusReader(dir, r'.*bbc.*\.txt',
                                                 cat_pattern=r'.*bbc_(\w+)\.txt')
       printer = dl.log_api.Printer(nelems=3)
       printer.print('Categories', reader.categories())
       printer.print('World fileids', reader.fileids(categories=['world']))
       printer.print('Technology fileids',
                     reader.fileids(categories=['technology']))
   

最终结果参见以下截图:

代码在本书代码包的corpus.py文件中。

另见

• 位于的CategorizedPlaintextCorpusReader文档

用句子和单词标记新闻文章

作为 NLTK 分布的一部分的语料库已经被标记化，所以我们可以很容易地获得单词和句子的列表。对于我们自己的语料库，我们也应该应用标记化。这个食谱演示了如何用 NLTK 实现标记化。我们将使用的文本文件在本书的代码包中。这个特殊的文本是英文的，但是 NLTK 也支持其他语言。

做好准备

按照本章简介部分的说明安装 NLTK。

怎么做...

程序在本书代码包的tokenizing.py文件中:

1. 进口情况如下:

   from nltk.tokenize import sent_tokenize
   from nltk.tokenize import word_tokenize
   import dautil as dl
   

2. 以下代码演示了标记化:

   fname = '46_bbc_world.txt'
   printer = dl.log_api.Printer(nelems=3)
   
   with open(fname, "r", encoding="utf-8") as txt_file:
       txt = txt_file.read()
       printer.print('Sentences', sent_tokenize(txt))
       printer.print('Words', word_tokenize(txt))
   

有关最终结果，请参考以下屏幕截图:

另见

• 相关文件在 www.nltk.org/api/nltk.to… = send _ tokenize # nltk . tokenize . send _ tokenize(2015 年 10 月检索)。

词干、引理、过滤和 TF-IDF 分数

单词包模型将语料库字面上表示为单词包，不考虑单词的位置，只考虑它们的数量。停止词是常见的词，如“a”、“is”和“the”，它们不增加信息价值。

TF-IDF 分数可以针对单个单词(单字)或多个连续单词(n-克)的组合进行计算。TF-IDF 大致为词频和逆文档频率的比值。我说“大致”是因为我们通常取比率的对数或应用加权方案。术语频率是文档中一个单词或 n-gram 的频率。逆文档频率是单词或 n-gram 出现的文档数量的倒数。我们可以使用 TF-IDF 分数进行聚类或者作为分类的一个特征。在用非负矩阵分解提取主题食谱中，我们将使用分数来发现主题。

NLTK 通过一个稀疏矩阵来表示分数，该矩阵为语料库中的每个文档提供一行，为每个单词或 n-gram 提供一列。即使矩阵是稀疏的，我们也应该根据我们试图解决的问题类型，尽可能多地过滤单词。过滤代码在ch8util.py中，执行以下操作:

• 将所有单词转换为小写。在英语中，句子以大写字母开头，在单词包模型中，我们不关心单词的位置。显然，如果我们想检测命名实体(如识别命名实体配方)，案例很重要。
• 忽略停止词，因为它们没有语义价值。
• 忽略仅由一个字符组成的单词，因为这些单词要么是停止单词，要么是假装单词的标点符号。
• 忽略只出现一次的单词，因为这些单词不太可能很重要。
• 只允许包含字母的单词，因此会忽略像“7”这样包含数字的单词。

我们也会用引理进行过滤。引理化类似于词干化，我也将演示这一点。这两个过程背后的想法是单词有共同的词根，例如，“分析”、“分析师”和“分析师”有共同的词根。一般来说，词干会剪切字符，所以结果不一定是一个有效的单词。相反，引理化总是产生有效的单词并执行字典查找。

本书代码包中ch8util.py文件的代码如下:

from collections import Counter
from nltk.corpus import brown
from joblib import Memory

memory = Memory(cachedir='.')

def only_letters(word):
    for c in word:
        if not c.isalpha():
            return False

    return True

@memory.cache
def filter(fid, lemmatizer, sw):
    words = [lemmatizer.lemmatize(w.lower()) for w in brown.words(fid)
             if len(w) > 1 and w.lower() not in sw]

    # Ignore words which only occur once
    counts = Counter(words)
    rare = set([w for w, c in counts.items() if c == 1])

    filtered_words = [w for w in words if w not in rare]

    return [w for w in filtered_words if only_letters(w)]

我决定将分析限制在单图，但是将分析扩展到二元模型或三元模型是非常容易的。我们将使用的 scikit-learn TfidfVectorizer类允许我们指定一个ngram_range字段，这样我们就可以同时考虑 unigrams 和 n-grams。我们将腌制这个食谱的结果，以供其他食谱重复使用。

做好准备

按照简介部分的说明安装 NLTK。

怎么做...

脚本在本书代码包的stemming_lemma.py文件中:

1. 进口情况如下:

   from nltk.corpus import brown
   from nltk.corpus import stopwords
   from nltk.stem import PorterStemmer
   from nltk.stem import WordNetLemmatizer
   import ch8util
   from sklearn.feature_extraction.text import TfidfVectorizer
   import numpy as np
   import pandas as pd
   import pickle
   import dautil as dl
   

2. 演示词干和引理如下:

   stemmer = PorterStemmer()
   lemmatizer = WordNetLemmatizer()
   
   print('stem(analyses)', stemmer.stem('analyses'))
   print('lemmatize(analyses)', lemmatizer.lemmatize('analyses'))
   

3. 过滤 NLTK 布朗语料库中的单词:

   sw = set(stopwords.words())
   texts = []
   
   fids = brown.fileids(categories='news')
   
   for fid in fids:
       texts.append(" ".join(ch8util.filter(fid, lemmatizer, sw)))
   

4. 计算 TF-IDF 分数如下:

   vectorizer = TfidfVectorizer()
   matrix = vectorizer.fit_transform(texts)
   
   with open('tfidf.pkl', 'wb') as pkl:
       pickle.dump(matrix, pkl)
   
   sums = np.array(matrix.sum(axis=0)).ravel()
   
   ranks = [(word, val) for word, val in
            zip(vectorizer.get_feature_names(), sums)]
   
   df = pd.DataFrame(ranks, columns=["term", "tfidf"])
   df.to_pickle('tfidf_df.pkl')
   df = df.sort(['tfidf'])
   dl.options.set_pd_options()
   print(df)
   

最终结果参见以下截图:

它是如何工作的

如您所见，词干不返回有效的单词。它比引理化更快；但是，如果想要重用结果，首选引理化是有意义的。TF-IDF 分数在最终 PandasDataFrame对象中按升序排序。更高的 TF-IDF 分数表示更重要的单词。

另见

• 维基百科页面关于en.wikipedia.org/wiki/Bag-of…的单词包模型(2015 年 10 月检索)
• 维基百科关于 en.wikipedia.org/wiki/Lemmat… 词条化的页面(2015 年 10 月检索)
• en.wikipedia.org/wiki/Tf%E2%… 的 TF-IDF 的维基百科页面(2015 年 10 月检索)
• 维基百科关于 en.wikipedia.org/wiki/Stop_w… 的停词页面(2015 年 10 月检索)
• 位于的TfidfVectorizer类的文档(2015 年 10 月检索)
• 位于的WordNetLemmatizer类文档

识别命名实体

命名实体识别 ( NER )试图检测文本中的人名、组织名、地点名和其他名称。有些 NER 系统几乎和人类一样好，但这不是一件容易的事情。命名实体通常以大写字母开头，例如 Ivan。因此，在应用 NER 时，我们不应该改变词语的大小写。

NLTK 支持斯坦福 NER 应用编程接口。这是一个 Java API，所以你的系统上需要有 Java。我用 Java 1.8.0_65 测试了代码。本食谱中的代码下载了截至 2015 年 10 月的最新斯坦福 NER 档案(stanford-ner-2015-04-20.zip/3.5.2)。如果你想要另一个版本，看看nlp.stanford.edu/software/CR…(检索 2015 年 10 月)。

做好准备

按照简介部分的说明安装 NLTK。您可能还需要安装 Java。

怎么做...

脚本在本书代码包的named_entity.py文件中:

1. 进口如下:

   from nltk.tag.stanford import NERTagger
   import dautil as dl
   import os
   from zipfile import ZipFile
   from nltk.corpus import brown
   

2. 定义以下功能下载 NER 档案:

   def download_ner():
       url = 'http://nlp.stanford.edu/software/stanford-ner-2015-04-20.zip'
       dir = os.path.join(dl.data.get_data_dir(), 'ner')
   
       if not os.path.exists(dir):
           os.mkdir(dir)
   
       fname = 'stanford-ner-2015-04-20.zip'
       out = os.path.join(dir, fname)
   
       if not dl.conf.file_exists(out):
           dl.data.download(url, out)
   
           with ZipFile(out) as nerzip:
               nerzip.extractall(path=dir)
   
       return os.path.join(dir, fname.replace('.zip', ''))
   

3. 将 NER 应用于布朗语料库中的一个文件:

   dir = download_ner()
   st = NERTagger(os.path.join(dir, 'classifiers',
                               'english.all.3class.distsim.crf.ser.gz'),
                  os.path.join(dir, 'stanford-ner.jar'))
   fid = brown.fileids(categories='news')[0]
   printer = dl.log_api.Printer(nelems=9)
   
   tagged = [pair for pair in dl.collect.flatten(st.tag(brown.words(fid)))
             if pair[1] != 'O']
   printer.print(tagged)
   

有关最终结果，请参考以下屏幕截图:

它是如何工作的

我们通过指定一个预先训练好的分类器和 NER JAR (Java 档案)创建了一个NerTagger对象。分类器在我们的语料库中将单词标记为组织、位置、人或其他。分类是区分大小写的，这意味着如果你把所有的单词都小写，你会得到不同的结果。

另见

• en.wikipedia.org/wiki/Named-… 关于 NER 的维基百科页面(2015 年 10 月检索)

利用非负矩阵分解提取主题

自然语言处理中的主题并不完全符合字典的定义，而是更符合一个模糊的统计概念。我们谈到话题模型和与话题相关的单词的概率分布，正如我们所知。当我们阅读文本时，我们期望出现在标题或正文中的某些单词能够捕捉文档的语义上下文。一篇关于 Python 编程的文章会有“类”和“函数”这样的词，而一个关于蛇的故事会有“蛋”和“害怕”这样的词文本通常有多个主题；例如，这个食谱是关于主题模型和非负矩阵分解，我们将很快讨论。因此，我们可以通过给主题分配不同的权重来为主题定义一个加法模型。

主题建模算法之一是非负矩阵分解 ( NMF )。该算法将一个矩阵分解为两个矩阵的乘积，使得两个矩阵没有负值。通常，我们只能在数值上逼近因式分解的解，时间复杂度是多项式的。scikit-learn NMF类实现了这个算法。NMF 还可以应用于文档聚类和信号处理。

怎么做...

我们将重用来自词干、引理、过滤和 TF-IDF 评分方法的结果:

1. 进口情况如下:

   from sklearn.decomposition import NMF
   import ch8util
   

2. 从泡菜加载 TF-IDF 矩阵和单词:

   terms = ch8util.load_terms()
   tfidf = ch8util.load_tfidf()
   

3. 将主题可视化为高级单词列表:

   nmf = NMF(n_components=44, random_state=51).fit(tfidf)
   
   for topic_idx, topic in enumerate(nmf.components_):
       label = '{}: '.format(topic_idx)
       print(label, " ".join([terms[i] for i in topic.argsort()[:-9:-1]]))
   

有关最终结果，请参考以下屏幕截图:

代码在本书代码包的topic_extraction.py文件中。

它是如何工作的

NMF类有一个components_属性，保存数据的非负成分。我们选择了components_属性中最高值对应的单词。如你所见，话题多种多样，虽然有点过时。

另见

• 位于的 NMF 课程的文档
• 关于en.wikipedia.org/wiki/Topic_…话题模型的维基百科页面(2015 年 10 月检索)
• 维基百科关于 NMF 的页面

实现基本术语数据库

众所周知，自然语言处理有许多应用:

• 商业和开源搜索引擎实现的全文搜索
• 文档聚类
• 分类，例如在产品评论的上下文中确定文本类型或情感

为了执行这些任务，我们需要计算特征，如 TF-IDF 分数(参考词干、引理、过滤和 TF-IDF 分数)。特别是在大数据集的情况下，存储特征以便于处理是有意义的。搜索引擎使用倒排索引，将单词映射到网页。这类似于关联表模式(参见实现关联表)。

我们将用三个表实现关联表模式。一个表包含单词，另一个表将实现三个表的关联表模式。一个表包含单词，另一个表保存关于文档的信息，第三个表链接其他两个表，如下图所示:

怎么做...

程序在本书代码包的terms_database.py文件中:

1. 进口为如下:

   from sqlalchemy.ext.declarative import declarative_base
   from sqlalchemy import Column
   from sqlalchemy import ForeignKey
   from sqlalchemy import Float
   from sqlalchemy import Integer
   from sqlalchemy import String
   from sqlalchemy.orm import backref
   from sqlalchemy.orm import relationship
   import os
   import dautil as dl
   from nltk.corpus import brown
   from sqlalchemy import func
   import ch8util
   
   Base = declarative_base()
   

2. 为文本文档定义以下类别:

   class Text(Base):
       __tablename__ = 'texts'
       id = Column(Integer, primary_key=True)
       file = Column(String, nullable=False, unique=True)
       terms = relationship('Term', secondary='text_terms')
   
       def __repr__(self):
           return "Id=%d file=%s" % (self.id, self.file)
   

3. 为文章中的单词定义以下类别:

   class Term(Base):
       __tablename__ = 'terms'
       id = Column(Integer, primary_key=True)
       word = Column(String, nullable=False, unique=True)
   
       def __repr__(self):
           return "Id=%d word=%s" % (self.id, self.word)
   

4. 为文档和单词的关联定义以下类:

   class TextTerm(Base):
       __tablename__ = 'text_terms'
       text_id = Column(Integer, ForeignKey('texts.id'), primary_key=True)
       term_id = Column(Integer, ForeignKey('terms.id'), primary_key=True)
       tf_idf = Column(Float)
       text = relationship('Text', backref=backref('term_assoc'))
       term = relationship('Term', backref=backref('text_assoc'))
   
       def __repr__(self):
           return "text_id=%s term_id=%s" % (self.text_id, self.term_id)
   

5. 定义以下功能在文本表中插入条目:

   def populate_texts(session):
       if dl.db.not_empty(session, Text):
           # Cowardly refusing to continue
           return
   
       fids = brown.fileids(categories='news')
   
       for fid in fids:
           session.add(Text(file=fid))
   
       session.commit()
   

6. 定义以下函数在术语表中插入条目:

   def populate_terms(session):
       if dl.db.not_empty(session, Term):
           # Cowardly refusing to continue
           return
   
       terms = ch8util.load_terms()
   
       for term in terms:
           session.add(Term(word=term))
   
       session.commit()
   

7. 定义以下函数在关联表中插入条目:

   def populate_text_terms(session):
       if dl.db.not_empty(session, TextTerm):
           # Cowardly refusing to continue
           return
   
       text_ids = dl.collect.flatten(session.query(Text.id).all())
       term_ids = dl.collect.flatten(session.query(Term.id).all())
   
       tfidf = ch8util.load_tfidf()
       logger = dl.log_api.conf_logger(__name__)
   
       for text_id, row, in zip(text_ids, tfidf):
           logger.info('Processing {}'.format(text_id))
           arr = row.toarray()[0]
           session.get_bind().execute(
               TextTerm.__table__.insert(),
               [{'text_id': text_id, 'term_id': term_id,
                 'tf_idf': arr[i]}
                for i, term_id in enumerate(term_ids)
                if arr[i] > 0]
           )
   
       session.commit()
   

8. 定义以下功能，使用关键词进行搜索:

   def search(session, keywords):
       terms = keywords.split()
       fsum = func.sum(TextTerm.tf_idf)
   
       return session.query(TextTerm.text_id, fsum).\
           join(Term, TextTerm).\
           filter(Term.word.in_(terms)).\
           group_by(TextTerm.text_id).\
           order_by(fsum.desc()).all()
   

9. 用下面的代码调用我们定义的函数:

   if __name__ == "__main__":
       dbname = os.path.join(dl.data.get_data_dir(), 'news_terms.db')
       session = dl.db.create_session(dbname, Base)
       populate_texts(session)
       populate_terms(session)
       populate_text_terms(session)
       printer = dl.log_api.Printer()
       printer.print('id, tf_idf', search(session, 'baseball game'))
   

我们搜索了“棒球比赛”有关最终结果(文件标识和 TF-IDF 总和)，请参考以下屏幕截图:

它是如何工作的

我们使用关联表数据库模式存储 TF-IDF 分数。作为使用数据库的一个例子，我们查询了“棒球比赛”查询在术语表中查找两个单词的标识，然后在关联表中汇总相关的 TF-IDF 分数。总和作为相关性得分。然后，我们给出了相应的文件标识与相关性得分的降序排列。如果您向最终用户显示结果，您必须至少再做一次查询，才能用文件名替换文件标识。碰巧，我们正在分析的文件被命名为ca01到ca44，所以查询不是严格必要的。

因为我已经有了 TF-IDF 分数，所以我发现直接存储它们很方便。但是，您也可以决定存储术语频率和反向文档频率，并从中得出 TF-IDF 分数。您只需要确定每个新文档的术语频率和文档中的单词。添加或删除文档时，需要更新所有反向文档频率。然而，通过关联表建立链接已经足以计算术语频率、反向文档频率和 TF-IDF 分数。

另见

• 维基百科中关于en.wikipedia.org/wiki/Search…搜索引擎索引的页面(2015 年 10 月检索)

计算社交网络密度

人类是群居动物，因此，社会关系非常重要。我们可以将这些联系和相关人员视为一个网络。我们将网络或子集表示为图形。图由通过边或线连接的节点或点组成。图形可以是有向的，也可以是无向的，线条可以是箭头。

我们将使用脸书 SPAN 数据，我们也在可视化网络图与蜂巢图配方中使用了该数据。脸书在 2004 年起步时规模很小，但截至 2015 年，其用户已超过 10 亿。数据没有包括所有的用户，但是对于一个像样的分析来说还是足够的。以下方程描述了无向(8.1)和有向(8.2)图的密度:

在这些等式中， n 是节点的数量， m 是边的数量。

做好准备

按照简介部分的说明安装网络。

怎么做...

代码在本书代码包的net_density.ipynb文件中:

1. 进口情况如下:

   import networkx as nx
   import dautil as dl
   

2. 创建一个网络图，如下所示:

   fb_file = dl.data.SPANFB().load()
   G = nx.read_edgelist(fb_file,
                        create_using=nx.Graph(),
                        nodetype=int)
   

3. 调用density()功能如下:

   print('Density', nx.density(G))
   

我们得到以下密度:

Density 0.010819963503439287

另见

• https://networkx . github . io/documentation/latest/reference/generated/networkx . class . function . density . html中记录的density()函数(2015 年 10 月检索)
• 关于图形的维基百科页面位于https://en . Wikipedia . org/wiki/Graph _ % 28 abstract _ data _ type % 29(2015 年 10 月检索)

计算社交网络亲密度中心度

在像脸书 SPAN 数据这样的社交网络中，我们会有有影响力的人。在图形术语中，这些是有影响力的节点。中心性找到重要节点的特征。贴近度中心度使用节点间最短路径作为特征，如下式所示:

在(8.3)中， d(u，v) 是 u 、 v 之间的最短路径， n 是节点数。一个有影响力的节点离其他节点很近，因此最短路径的总和很低。我们可以分别计算每个节点的贴近度中心度，对于一个大图来说，这可能是一个冗长的计算。NetworkX 允许我们指定我们感兴趣的节点，所以我们将只计算几个节点的接近度中心度。

做好准备

按照简介部分的说明安装网络。

怎么做...

请看一下本书代码包中的close_centrality.ipynb文件:

1. 进口情况如下:

   import networkx as nx
   import dautil as dl
   

2. 根据脸书 SPAN 数据创建一个网络图，如下所示:

   fb_file = dl.data.SPANFB().load()
   G = nx.read_edgelist(fb_file,
                        create_using=nx.Graph(),
                        nodetype=int)
   

3. 计算节点1和节点4037的贴近度中心度:

   print('Closeness Centrality Node 1',
         nx.closeness_centrality(G, 1))
   print('Closeness Centrality Node 4037',
         nx.closeness_centrality(G, 4037))
   

对于脸书跨度数据，我们得到以下结果:

Closeness Centrality Node 1 0.2613761408505405
Closeness Centrality Node 4037 0.18400546821599453

另见

• 维基百科关于亲密度中心性的页面在https://en . Wikipedia . org/wiki/Centrality #亲密度 _ 中心性(2015 年 10 月检索)

确定中间中心性

介数中心性是一种类似于接近中心性的中心性(参考计算社交网络接近中心性食谱)。该指标由以下等式给出:

它是通过一个节点的所有可能的最短路径对的分数的总和。

做好准备

按照简介部分的说明安装网络。

怎么做...

脚本在本书代码包的between_centrality.ipynb文件中:

1. 进口情况如下:

   import networkx as nx
   import dautil as dl
   import pandas as pd
   

2. 将脸书 SPAN 数据载入网络图:

   fb_file = dl.data.SPANFB().load()
   G = nx.read_edgelist(fb_file,
                        create_using=nx.Graph(),
                        nodetype=int)
   

3. 用k = 256(要使用的节点数)计算中间中心度，并将结果存储在 PandasDataFrame对象中:

   key_values = nx.betweenness_centrality(G, k=256)
   df = pd.DataFrame.from_dict(key_values, orient='index')
   
   dl.options.set_pd_options()
   print('Betweenness Centrality', df)
   

有关最终结果，请参考以下屏幕截图:

另见

• 维基百科页面，关于en.wikipedia.org/wiki/Betwee…的中间性中心性(检索于 2015 年 10 月
• 位于https://networkx . github . io/documentation/latest/reference/generated/networkx . algorithms . centrality . interness _ centrality . html的betweenness_centrality()功能的文档(2015 年 10 月检索)

估计平均聚类系数

从幼儿园开始，我们就有了朋友，亲密的朋友，永远最好的朋友，社交媒体的朋友，以及其他的朋友。因此，社交网络图应该有团块，不像你在高中派对上看到的那样。自然产生的问题是，如果我们只是邀请一群随机的陌生人参加聚会，或者在网上重新创建这个设置，会发生什么？我们希望陌生人联系的概率比朋友低。在图论中，这个概率是由 聚类系数来衡量的。

平均聚类系数是聚类系数的局部(单节点)版本。这个度量的定义考虑了由节点组成的三角形。有了三个节点，我们可以形成一个三角形，例如，三个火枪手。如果我们在混合中加入达塔格南，更多的三角形是可能的，但不是所有的三角形都必须实现。达塔尼昂可能会和三个火枪手打起来。在(8.5)中，我们将聚类系数定义为已实现三角形和可能三角形的比率以及平均聚类系数(8.6):

做好准备

按照简介部分的说明安装网络。

怎么做...

脚本在本书代码包的avg_clustering.ipynb文件中:

1. 进口情况如下:

   import networkx as nx
   import dautil as dl
   

2. 将脸书跨度数据加载到网络图中:

   fb_file = dl.data.SPANFB().load()
   G = nx.read_edgelist(fb_file,
                        create_using=nx.Graph(),
                        nodetype=int)
   

3. 计算平均聚类系数如下:

   print('Average Clustering',
         nx.average_clustering(G))
   

对于脸书跨度数据，我们得到以下结果:

Average Clustering 0.6055467186200871

另见

• 关于 en.wikipedia.org/wiki/Cluste… 聚类系数的维基百科页面(2015 年 10 月检索)
• 位于的average_clustering()功能的文档。

计算图的分类系数

在图论中，相似度是通过度分布来衡量的。度是一个节点到其他节点的连接数。在有向图中，我们有传入和传出的连接以及相应的度数和度数。朋友往往有共同点。在图论中，这个趋势由配合度系数来衡量。该系数是一对节点之间的皮尔逊相关系数，如下式所示:

qk (剩余度分布)是离开节点 k 的连接数。 ejk 是节点对剩余度数的联合概率分布。

做好准备

按照简介部分的说明安装网络。

怎么做...

代码在本书代码包的assortativity.ipynb文件中:

1. 进口情况如下:

   import networkx as nx
   import dautil as dl
   

2. 将脸书跨度数据加载到网络图中:

   fb_file = dl.data.SPANFB().load()
   G = nx.read_edgelist(fb_file,
                        create_using=nx.Graph(),
                        nodetype=int)
   

3. 计算分类系数如下:

   print('Degree Assortativity Coefficient',
         nx.degree_assortativity_coefficient(G))
   

对于脸书跨度数据，我们得到以下结果:

Degree Assortativity Coefficient 0.0635772291856

另见

• 关于【en.wikipedia.org/wiki/Degree… 年 10 月检索)](en.wikipedia.org/wiki/Degree…)
• en.wikipedia.org/wiki/Assort… 的分类的维基百科页面(2015 年 10 月检索)
• 位于的degree_assortativity_coefficient()功能的文档

求图的团数

一个完全图是一个在中的图，其中每对节点通过唯一的连接连接。一个小团体 是一个完整的子图。这相当于拉帮结派的一般概念，每个人都认识所有其他人。最大团是节点最多的团。团数是最大团中的节点数。不幸的是，找到团数需要很长时间，所以我们不会使用完整的脸书 SPAN 数据。

做好准备

按照简介部分的说明安装网络。

怎么做...

代码在本书代码包的clique_number.py文件中:

1. 进口情况如下:

   import networkx as nx
   import dautil as dl
   

2. 将脸书跨度数据加载到网络图中:

   fb_file = dl.data.SPANFB().load()
   G = nx.read_edgelist(fb_file,
                        create_using=nx.Graph(),
                        nodetype=int)
   

3. 确定子图的团数:

   print('Graph Clique Number',
         nx.graph_clique_number(G.subgraph(list(range(2048)))))
   

对于部分脸书跨度数据，我们得到以下结果:

Graph Clique Number 38

另见

• 维基百科关于 en.wikipedia.org/wiki/Comple… T2 完整图表的页面(2015 年 10 月检索)
• 维基百科关于 en.wikipedia.org/wiki/Clique… 派系的页面(2015 年 10 月检索)
• https://networkx . github . io/documentation/latest/reference/generated/networkx . algorithms . click . graph _ click _ number . html上的功能文档(2015 年 10 月检索)

创建具有余弦相似性的文档图

互联网是一个相互链接的大型文档网络。我们可以将其视为文档图，其中每个节点对应一个文档。您将期望文档链接到类似的文档；但是，网页有时会链接到其他不相关的网页。这可能是错误的或故意的，例如在广告或试图提高搜索引擎排名的背景下。像维基百科这样更值得信赖的来源可能会产生更好的图表。然而，一些维基百科页面是非常基本的存根，所以我们可能会错过高质量的链接。

余弦相似度是一种常用的距离度量，用来衡量两个文档的相似度。对于这个度量，我们需要计算两个特征向量的内积。向量的余弦相似性对应于向量之间角度的余弦，因此得名。余弦相似度由以下等式给出:

该配方中的特征向量是对应于文档的 TF-IDF 分数。文档与其自身的余弦相似度等于 1(零角度)；因此，对于相似的文档，余弦相似度应该尽可能接近 1。

我们将执行以下步骤来创建布朗语料库中新闻文章的文档图:

1. 使用我们从词干、引理、过滤和配方的 TF-IDF 分数代码中存储的 TF-IDF 分数计算余弦相似度。结果类似于相关矩阵。
2. 对于每个文档，在图中为每个文档添加一个连接，这是非常相似的。我用相似度的第 90 个百分位数作为阈值；但是，如果您愿意，可以使用另一个值。
3. 对于每个文档，使用 TF-IDF 分数选择前三个单词。我用这些词来注释文档节点。
4. 如本章所述，使用网络计算图形指标。

怎么做...

创建余弦相似性文档图的代码在本书代码包的cos_similarity.ipynb文件中:

1. 进口情况如下:

   from sklearn.metrics.pairwise import cosine_similarity
   import networkx as nx
   import matplotlib.pyplot as plt
   import numpy as np
   import dautil as dl
   import ch8util
   

2. 定义以下功能，将节点添加到每个文档前三个最重要的单词注释的网络图中:

   def add_nodes(G, nodes, start, terms):
       for n in nodes:
           words = top_3_words(tfidf, n, terms)
           G.add_node(n, words='{0}: {1}'.
                      format(n, " ".join(words.tolist())))
           G.add_edge(start, n)
   

3. 定义以下函数来查找文档的前三个单词:

   def top_3_words(tfidf, row, terms):
       indices = np.argsort(tfidf[row].toarray().ravel())[-3:]
   
       return terms[indices]
   

4. 加载必要的数据，计算余弦相似度，并创建一个网络图:

   tfidf = ch8util.load_tfidf()
   terms = ch8util.load_terms()
   
   sims = cosine_similarity(tfidf, tfidf)
   G = nx.Graph()
   

5. 迭代余弦相似度，并向图中添加节点:

   for i, row in enumerate(sims):
       over_limit = np.where(row > np.percentile(row, 90))[0]
       nodes = set(over_limit.tolist())
       nodes.remove(i)
       add_nodes(G, nodes, i, terms)
   

6. 绘制图表并使用网络打印一些指标:

   labels = nx.get_node_attributes(G, 'words')
   nx.draw_networkx(G, pos=nx.spring_layout(G), labels=labels)
   plt.axis('off')
   plt.title('Graph of News Articles in the Brown Corpus')
   print('Density', nx.density(G))
   print('Average Clustering',
         nx.average_clustering(G))
   print('Degree Assortativity Coefficient',
         nx.degree_assortativity_coefficient(G))
   print('Graph Clique Number', nx.graph_clique_number(G))
   

最终结果参见以下截图:

另见

• 计算社交网络密度的方法
• 估算平均聚类系数配方
• 计算图表的配合系数配方
• 获取图表配方的团数
• 维基百科关于en.wikipedia.org/wiki/Cosine…余弦相似度的页面(2015 年 10 月检索)
• 关于功能的文档(2015 年 10 月检索)

九、集成学习和降维

在本章中，我们将介绍以下食谱:

• 递归消除特征
• 应用主成分分析进行降维
• 应用线性判别分析进行降维
• 多个模型的堆叠和多数投票
• 随机森林学习
• 用 RANSAC 算法拟合噪声数据
• 打包以提高结果
• 促进更好的学习
• 嵌套交叉验证
• 使用 joblib 重用模型
• 分层聚类数据
• 进行一次大西洋之旅

简介

在 1983 年的战争游戏电影中，一台电脑做出了可能导致第三次世界大战的生死决定。据我所知，当时的技术还不能完成这样的壮举。然而，在 1997 年，深蓝超级计算机确实击败了一位世界象棋冠军。2005 年，一辆斯坦福自动驾驶汽车在沙漠中独自行驶了 130 多公里。2007 年，另一个车队的车在正常交通中行驶了 50 多公里。2011 年，沃森电脑在与人类对手的竞赛中获胜。如果我们假设计算机硬件是限制因素，那么我们可以尝试推断未来。雷·库兹韦尔做到了这一点，根据他的说法，我们可以预计人类水平的智力将在 2029 年左右。

在这一章中，我们将集中讨论更简单的预测第二天天气的问题。我们将假设今天的天气取决于昨天的天气。理论上，如果一只蝴蝶在一个位置扇动翅膀，这可能会引发一连串的事件，在几千公里外的地方引发一场暴风雪(蝴蝶效应)。这并非不可能，但非常不可能。然而，如果我们有很多这样的事件，类似的情况会比你想象的更频繁地发生。

不可能考虑所有可能的因素。事实上，我们会试图通过忽略一些我们现有的数据来让我们的生活变得更容易。我们将应用分类和回归算法，以及层次聚类。让我们将结果评估推迟到第 10 章、评估分类器、回归器和聚类。如果你对分类食谱中提到的混淆矩阵感到好奇，请直接跳到用混淆矩阵食谱进行分类。

事实上，如今大多数人工智能系统并不那么智能。法庭上的法官可能会做出错误的决定，因为他或她有偏见或者过得不好。一组多位评委应该表现更好。这相当于一个机器学习项目，我们担心过度拟合和欠拟合。集成学习就是这个难题的解决方案，它基本上意味着以一种巧妙的方式组合多个学习者。

本章的主要部分是关于 超参数优化——这些是分类器和回归器的参数。为了检查过拟合或欠拟合，我们可以使用学习曲线，该曲线显示训练和不同训练集大小的测试分数。我们也可以通过验证曲线改变单个超参数的值。

递归消除特征

如果我们有许多特性(解释变量)，那么将它们都包含在我们的模型中是很有诱惑力的。然而，我们会冒过度拟合的风险——得到一个对训练数据非常有效而对看不见的数据非常无效的模型。不仅如此，模型必然会相对较慢，需要大量内存。我们必须权衡准确性(或其他指标)与速度和内存需求。

我们可以尝试忽略特征或创建新的更好的复合特征。例如，在在线广告中，通常使用比率，例如与广告相关的浏览量和点击量的比率。常识或领域知识可以帮助我们选择特征。在最坏的情况下，我们可能不得不依赖相关性或其他统计方法。scikit-learn 库提供了RFE类(递归特征消除)，可以自动选择特征。我们将在这个食谱中使用这个课程。我们还需要一个外部估计器。RFE类是相对较新的，不幸的是，不能保证所有的评估人员都能与RFE类一起工作。

怎么做...

1. 进口情况如下:

   from sklearn.feature_selection import RFE
   from sklearn.svm import SVC
   from sklearn.svm import SVR
   from sklearn.preprocessing import MinMaxScaler
   import dautil as dl
   import warnings
   import numpy as np
   

2. 创建一个 SVC 分类器和一个 RFE 对象，如下所示:

   warnings.filterwarnings("ignore", category=DeprecationWarning)
   clf = SVC(random_state=42, kernel='linear')
   selector = RFE(clf)
   

3. 加载数据，使用MinMaxScaler功能进行缩放，并添加一年中的某一天作为特征:

   df = dl.data.Weather.load().dropna()
   df['RAIN'] = df['RAIN'] == 0
   df['DOY'] = [float(d.dayofyear) for d in df.index]
   scaler = MinMaxScaler()
   
   for c in df.columns:
       if c != 'RAIN':
           df[c] = scaler.fit_transform(df[c])
   

4. 打印数据的第一行作为健全性检查:

   dl.options.set_pd_options()
   print(df.head(1))
   X = df[:-1].values
   np.set_printoptions(formatter={'all': '{:.3f}'.format})
   print(X[0])
   np.set_printoptions()
   

5. 使用有雨或无雨作为类别来确定对功能的支持和排名(在分类的上下文中):

   y = df['RAIN'][1:].values
   selector = selector.fit(X, y)
   print('Rain support', df.columns[selector.support_])
   print('Rain rankings', selector.ranking_)
   

6. 使用温度作为特征来确定特征的支持度和等级:

   reg = SVR(kernel='linear')
   selector = RFE(reg)
   y = df['TEMP'][1:].values
   selector = selector.fit(X, y)
   print('Temperature support', df.columns[selector.support_])
   print('Temperature ranking', selector.ranking_)
   

最终结果参见以下截图:

本食谱的代码在本书代码包的feature_elimination.py文件中。

它是如何工作的

RFE类默认选择一半的特征。算法如下:

1. 根据数据训练外部估计器，并为特征分配权重。
2. 权重最小的要素将被移除。
3. 重复该过程，直到我们拥有必要数量的功能。

另见

• 位于的 RFE课程的文档。RFE.html(2015 年 11 月检索)

应用主成分分析进行降维

主成分分析 ( PCA )由卡尔·皮尔逊于 1901 年发明，是一种将数据转换为不相关的正交特征的算法，称为主成分。主分量是协方差矩阵的特征向量。

有时，我们通过在应用主成分分析之前缩放数据来获得更好的结果，尽管这不是严格必要的。我们可以将主成分分析解释为将数据投影到更低维度的空间。一些主成分贡献的信息相对较少(低方差)；因此，我们可以省略它们。我们有以下转变:

结果是矩阵 T _L ，行数与原始矩阵相同，但列数较低。

当然，降维对于可视化和建模以及减少过度拟合的机会是有用的。事实上，有一种技术叫做主成分回归 ( 聚合酶链反应)，它利用了这个原理。简而言之，聚合酶链反应执行以下步骤:

1. 使用主成分分析将数据转换到低维空间。
2. 在新空间中执行线性回归。
3. 将结果转换回原始坐标系。

怎么做...

1. 进口情况如下:

   import dautil as dl
   from sklearn.decomposition import PCA
   import matplotlib.pyplot as plt
   from sklearn.preprocessing import scale
   

2. 按如下方式加载数据，并按一年中的某一天分组:

   df = dl.data.Weather.load().dropna()
   df = dl.ts.groupby_yday(df).mean()
   X = df.values
   

3. 应用主成分分析将数据投影到二维空间:

   pca = PCA(n_components=2)
   X_r = pca.fit_transform(scale(X)).T
   

4. 绘制转换结果:

   plt.scatter(X_r[0], X_r[1])
   plt.xlabel('x')
   plt.ylabel('y')
   plt.title('Dimension Reducion with PCA')
   

有关最终结果，请参考以下屏幕截图:

代码在本书代码包的applying_pca.ipynb文件中。

另见

• en.wikipedia.org/wiki/Princi… 关于主成分分析的维基百科页面(2015 年 11 月检索)
• 关于聚合酶链反应的维基百科页面，网址为(2015 年 11 月检索)
• 位于的主成分分析课程的文档

应用线性判别分析进行降维

线性判别分析 ( LDA )是一种寻找特征线性组合以区分类别的算法。通过投影到低维子空间，它可以用于分类或降维。LDA 需要一个目标属性来进行分类和降维。

如果我们将类密度表示为多元高斯分布，那么线性判别分析假设类具有相同的协方差矩阵。我们可以使用训练数据来估计类分布的参数。

在 scikit-learn 中，lda.LDA在 0.17 中已被弃用，并重新命名为discriminant_analysis.LinearDiscriminantAnalysis。这个类的默认解算器使用奇异值分解，不需要计算协方差矩阵，因此速度很快。

怎么做...

代码在本书代码包的applying_lda.ipynb文件中:

1. 进口情况如下:

   import dautil as dl
   from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
   import matplotlib.pyplot as plt
   

2. 加载数据如下:

   df = dl.data.Weather.load().dropna()
   X = df.values
   y = df['WIND_DIR'].values
   

3. 应用 LDA 将数据投影到二维空间:

   lda = LinearDiscriminantAnalysis(n_components=2)
   X_r = lda.fit(X, y).transform(X).T
   

4. 绘制转换结果:

   plt.scatter(X_r[0], X_r[1])
   plt.xlabel('x')
   plt.ylabel('y')
   plt.title('Dimension Reduction with LDA')
   

有关最终结果，请参考以下屏幕截图:

另见

• en.wikipedia.org/wiki/Linear… 关于 LDA 的维基百科页面(检索到【2015 年 11 月)
• 相关 scikit-learn 文档为http://sci kit-learn . org/stable/modules/generated/sklearn . discriminal _ analysis。LinearDiscriminantAnalysis.html(2015 年 11 月检索)

多模型叠加和多数投票

一般认为，两个人单独认识的不止一个人。民主应该比独裁更有效。在机器学习中，我们没有人类做决定，而是算法。当我们有多个分类器或回归器一起工作时，我们说的是 集成学习。

有很多集成学习方案。最简单的设置是对分类进行多数投票，对回归进行平均。在 scikit-learn 0.17 中，可以使用VotingClassifier类进行多数投票。该分类器允许您使用权重来强调或抑制分类器。

堆叠获取机器学习估计器的输出，然后将这些输出用作另一个算法的输入。当然，您可以将高级算法的输出馈送给另一个预测器。可以使用任意拓扑，但出于实际原因，您应该首先尝试简单的设置。

怎么做...

1. 进口情况如下:

   import dautil as dl
   from sklearn.tree import DecisionTreeClassifier
   import numpy as np
   import ch9util
   from sklearn.ensemble import VotingClassifier
   from sklearn.grid_search import GridSearchCV
   from IPython.display import HTML
   

2. 加载数据并创建三个决策树分类器:

   X_train, X_test, y_train, y_test = ch9util.rain_split()
   default = DecisionTreeClassifier(random_state=53, min_samples_leaf=3,
                                    max_depth=4)
   entropy = DecisionTreeClassifier(criterion='entropy',
                                    min_samples_leaf=3, max_depth=4,
                                    random_state=57)
   random = DecisionTreeClassifier(splitter='random', min_samples_leaf=3,
                                   max_depth=4, random_state=5)
   

3. 使用分类器进行投票:

   clf = VotingClassifier([('default', default), 
                           ('entropy', entropy), ('random', random)])
   params = {'voting': ['soft', 'hard'],
            'weights': [None, (2, 1, 1), (1, 2, 1), (1, 1, 2)]}
   gscv = GridSearchCV(clf, param_grid=params, n_jobs=-1, cv=5)
   gscv.fit(X_train, y_train)
   votes = gscv.predict(X_test)
   
   preds = []
   
   for clf in [default, entropy, random]:
       clf.fit(X_train, y_train)
       preds.append(clf.predict(X_test))
   
   preds = np.array(preds)
   

4. 绘制基于投票的预测的混淆矩阵:

   %matplotlib inline
   context = dl.nb.Context('stacking_multiple')
   dl.nb.RcWidget(context)
   
   sp = dl.plotting.Subplotter(2, 2, context)
   html = ch9util.report_rain(votes, y_test, gscv.best_params_, sp.ax)
   sp.ax.set_title(sp.ax.get_title() + ' | Voting')
   

5. 绘制基于堆叠的预测的混淆矩阵:

   default.fit(preds_train.T, y_train)
   stacked_preds = default.predict(preds.T)
   html += ch9util.report_rain(stacked_preds, 
                               y_test, default.get_params(), sp.next_ax())
   sp.ax.set_title(sp.ax.get_title() + ' | Stacking')
   ch9util.report_rain(default.predict(preds.T), y_test)
   

6. 绘制投票和堆叠分类器的学习曲线【T1:

   ch9util.plot_learn_curve(sp.next_ax(), gscv.best_estimator_, X_train,
                            y_train, title='Voting')
   
   ch9util.plot_learn_curve(sp.next_ax(), default, X_train,
                            y_train, title='Stacking')
   

有关最终结果，请参考以下屏幕截图:

代码是本书代码包中stacking_multiple.ipynb文件中的。

另见

• 位于的 VotingClassifier课程的文档
• 维基百科关于en.wikipedia.org/wiki/Ensemb…堆积的部分(2015 年 11 月检索)

随机森林学习

if a: else b语句是 Python 编程中最常见的语句之一。通过嵌套和组合这样的语句，我们可以构建一个所谓的决策树。这类似于老式的流程图，尽管流程图也允许循环。决策树在机器学习中的应用叫做 决策树学习。决策树学习中的树的末端节点，也称为叶子，包含一个分类问题的类标签。每个非叶节点都与一个包含特征值的布尔条件相关联。

决策树可以用来推导相对简单的规则。能够产生这样的结果，当然是一个巨大的优势。然而，你不得不怀疑这些规则有多好。如果我们添加新数据，我们会得到相同的规则吗？

如果一棵决策树是好的，那么整个森林应该会更好。多棵树应该可以减少过度拟合的机会。然而，在真实的森林中，我们不希望只有一种类型的树。显然，我们将不得不平均或通过多数投票决定什么是适当的结果。

在这个食谱中，我们将应用利奥·布雷曼和阿黛尔·卡特勒发明的 随机森林算法。名称中的“随机”是指从数据中随机选择特征。我们使用所有数据，但不在同一个决策树中。

随机森林也应用 装袋 ( 引导聚集)，我们将在装袋中讨论以提高结果配方。决策树的打包包括以下步骤:

1. 带有替换的示例训练示例，并将其分配到树中。
2. 根据指定的数据训练树。

我们可以通过交叉验证或绘制测试和训练误差相对于树的数量来确定正确的树的数量。

怎么做...

代码在本书代码包的random_forest.ipynb文件中:

1. 进口情况如下:

   import dautil as dl
   from sklearn.grid_search import GridSearchCV
   from sklearn.ensemble import RandomForestClassifier
   import ch9util
   import numpy as np
   from IPython.display import HTML
   

2. 加载数据，做如下预测:

   X_train, X_test, y_train, y_test = ch9util.rain_split()
   clf = RandomForestClassifier(random_state=44)
   params = {
       'max_depth': [2,  4],
       'min_samples_leaf': [1, 3],
       'criterion': ['gini', 'entropy'],
       'n_estimators': [100, 200]
   }
   
   rfc = GridSearchCV(estimator=RandomForestClassifier(),
                      param_grid=params, cv=5, n_jobs=-1)
   rfc.fit(X_train, y_train)
   preds = rfc.predict(X_test)
   

3. 将降雨预报混淆矩阵绘制如下:

   sp = dl.plotting.Subplotter(2, 2, context)
   html = ch9util.report_rain(preds, y_test, rfc.best_params_, sp.ax)
   

4. 绘制一系列森林大小的验证曲线:

   ntrees = 2 ** np.arange(9)
   ch9util.plot_validation(sp.next_ax(), rfc.best_estimator_, 
                           X_train, y_train, 'n_estimators', ntrees)
   

5. 绘制深度范围的验证曲线:

   depths = np.arange(2, 9)
   ch9util.plot_validation(sp.next_ax(), rfc.best_estimator_, 
                           X_train, y_train, 'max_depth', depths)
   

6. 绘制最佳估计量的学习曲线:

   ch9util.plot_learn_curve(sp.next_ax(), 
                            rfc.best_estimator_, X_train,y_train)
   HTML(html + sp.exit())
   

有关最终结果，请参考以下屏幕截图:

还有更多…

随机森林分类被认为是一种通用的算法，我们几乎可以将其用于任何分类任务。遗传 算法和遗传编程一般可以做一个网格搜索或者优化。

我们可以把程序看作是一个生成结果的操作符和操作数的序列。当然，这是一个非常简化的编程模型。然而，在这样的模型中，有可能使用模仿生物学理论的自然选择来进化程序。遗传程序是自我修改的，具有巨大的适应性，但我们得到的决定论水平较低。

TPOT 项目试图发展机器学习管道(目前使用少量分类器，包括随机森林)。我在 GitHub 上分叉了 TPOT 0.1.3 并做了一些修改。TPOT 使用deap作为遗传编程部件，可以安装如下:

$ pip install deap

我用 deap 1.0.2 测试了代码。在标签r1下安装我的更改，如下所示:

$ git clone git@github.com:ivanidris/tpot.git
$ cd tpot
$ git checkout r1
$ python setup.py install

您也可以从github.com/ivanidris/t…获取代码。本书代码包中rain_pot.py文件的以下代码演示了如何用 TPOT 拟合和评分降雨预测:

import ch9util
from tpot import TPOT

X_train, X_test, y_train, y_test = ch9util.rain_split()
tpot = TPOT(generations=7, population_size=110, verbosity=2)
tpot.fit(X_train, y_train)
print(tpot.score(X_train, y_train, X_test, y_test))

另见

• 维基百科关于 en.wikipedia.org/wiki/Random… 随机森林的页面(2015 年 11 月检索)
• 位于的RandomForestClassifier类的文档(2015 年 11 月检索)

用 RANSAC 算法拟合噪声数据

我们在本书的其他地方讨论了回归背景下的异常值问题(参见本食谱末尾的也参见部分)。问题很明显——异常值使得我们很难恰当地拟合我们的模型。随机样本一致性算法 ( RANSAC )尽最大努力以迭代方式拟合我们的数据。RANSAC 是由 Fishler 和 Bolles 在 1981 年推出的。

我们通常对我们的数据有一些了解，例如数据可能遵循正态分布。或者，数据可以是由具有不同特征的多个过程产生的混合。由于数据转换中的故障或错误，我们也可能有异常数据。在这种情况下，识别异常值并对其进行适当处理应该很容易。RANSAC 算法不知道你的数据，但它也假设有内联和外联。

算法经过固定次数的迭代。目标是找到一组指定大小的内联集(共识集)。

RANSAC 执行以下步骤:

1. 随机选择尽可能小的数据子集并拟合模型。
2. 检查每个数据点是否与上一步拟合的模型一致。使用残差阈值将不一致的点标记为异常值。
3. 如果找到足够多的内联，则接受该模型。
4. 用完全一致集重新估计参数。

scikit-learn RANSACRegressor类可以使用合适的估计器进行拟合。我们将使用默认的LinearRegression估计器。我们还可以指定拟合的最小样本数、残差阈值、异常值的决策函数、决定模型是否有效的函数、最大迭代次数以及共识集中所需的内联数。

怎么做...

代码在本书代码包的fit_ransac.ipynb文件中:

1. 进口情况如下:

   import ch9util
   from sklearn import linear_model
   from sklearn.grid_search import GridSearchCV
   import numpy as np
   import dautil as dl
   from IPython.display import HTML
   

2. 加载数据并进行温度预测，如下所示:

   X_train, X_test, y_train, y_test = ch9util.temp_split()
   ransac = linear_model.RANSACRegressor(random_state=27)
   params = {
       'max_trials': [50, 100, 200],
       'stop_probability': [0.98, 0.99]
   }
   
   gscv = GridSearchCV(estimator=ransac, param_grid=params, cv=5)
   gscv.fit(X_train, y_train)
   preds = gscv.predict(X_test)
   

3. 将预测值与实际值进行散点图:

   sp = dl.plotting.Subplotter(2, 2, context)
   html = ch9util.scatter_predictions(preds, y_test, gscv.best_params_,
                                      gscv.best_score_, sp.ax)
   

4. 绘制一系列试验号的验证曲线:

   trials = 10 * np.arange(5, 20)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'max_trials', trials)
   

5. 绘制一系列停止概率的验证曲线:

   probs = 0.01 * np.arange(90, 99)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'stop_probability', probs)
   

6. 绘制一系列一致集合大小的验证曲线:

   ninliers = 2 ** np.arange(4, 14)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'stop_n_inliers', ninliers)
   HTML(html + sp.exit())
   

有关最终结果，请参考以下屏幕截图:

另见

• en.wikipedia.org/wiki/RANSAC 关于 RANSAC 算法的维基百科页面(2015 年 11 月检索)
• 相关 scikit-learn 文档位于http://sci kit-learn . org/stable/modules/generated/sklearn . linear _ model。RANSACRegressor.html(2015 年 11 月检索)
• 拟合稳健线性模型的方法
• 用加权最小二乘法考虑方差的配方

套袋提高效果

Bootstrap 聚合或 bagging 是 Leo Breiman 在 1994 年引入的一种算法，它将 Bootstrap 应用于机器学习问题。随机森林学习食谱中也提到了套袋。

该算法旨在通过以下步骤减少过拟合的机会:

1. 我们通过替换采样从输入训练数据生成新的训练集。
2. 使模型适合每个生成的训练集。
3. 通过平均或多数投票来组合模型的结果。

scikit-learn BaggingClassifier类允许我们引导训练示例，我们还可以引导随机森林算法中的特征。当我们执行网格搜索时，我们引用前缀为 base_estimator__的基本估计量的超参数。我们将使用决策树作为基础估计，这样我们就可以重用来自随机森林学习方法的一些超参数配置。

怎么做...

代码在本书代码包的bagging.ipynb文件中:

1. 进口情况如下:

   import ch9util
   from sklearn.ensemble import BaggingClassifier
   from sklearn.grid_search import GridSearchCV
   from sklearn.tree import DecisionTreeClassifier
   import numpy as np
   import dautil as dl
   from IPython.display import HTML
   

2. 加载数据并创建BaggingClassifier :

   X_train, X_test, y_train, y_test = ch9util.rain_split()
   clf = BaggingClassifier(base_estimator=DecisionTreeClassifier(
       min_samples_leaf=3, max_depth=4), random_state=43)
   

3. 网格搜索、拟合和预测如下:

   params = {
       'n_estimators': [320, 640],
       'bootstrap_features': [True, False],
       'base_estimator__criterion': ['gini', 'entropy']
   }
   
   gscv = GridSearchCV(estimator=clf, param_grid=params,
                       cv=5, n_jobs=-1)
   
   gscv.fit(X_train, y_train)
   preds = gscv.predict(X_test)
   

4. 将降雨预报混淆矩阵绘制如下:

   sp = dl.plotting.Subplotter(2, 2, context)
   html = ch9util.report_rain(preds, y_test, gscv.best_params_, sp.ax)
   

5. 绘制一系列集合大小的验证曲线:

   ntrees = 2 ** np.arange(4, 11)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'n_estimators', ntrees)
   

6. 绘制max_samples参数的验证曲线:

   nsamples = 2 ** np.arange(4, 14)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'max_samples', nsamples)
   

7. 绘制学习曲线如下:

   ch9util.plot_learn_curve(sp.next_ax(), gscv.best_estimator_, 
                            X_train, y_train)
   HTML(html + sp.exit())
   

有关最终结果，请参考以下屏幕截图:

另见

• en.wikipedia.org/wiki/Bootst… 装袋的维基百科页面(2015 年 11 月检索)
• 位于的 BaggingClassifier文档

促进更好的学习

数量上的优势是大国比小国更成功的原因。这并不意味着一个大国的人生活得更好。但是从全局来看，个人并没有那么重要，就像在决策树的集合中，如果我们有足够多的树，单棵树的结果可以被忽略。

在分类的背景下，我们将 **【弱学习者】**定义为只比基线好一点的学习者，例如随机分配班级。虽然学习能力弱的人个人也很弱，就像蚂蚁一样，但是在一起，他们可以像蚂蚁一样做惊人的事情。

使用权重考虑每个学习者的力量是有意义的。这个总的思路叫做 助推。有很多助推算法，其中我们将在本食谱中使用 AdaBoost 。增强算法的不同之处主要在于它们的加权方案。

AdaBoost 使用加权求和来产生最终结果。这是一种自适应算法，试图提高单个训练示例的结果。如果你已经为考试而学习，你可能已经应用了类似的技术，识别你有问题的问题类型，并专注于困难的问题。在 AdaBoost 的情况下，提升是通过调整弱学习者来完成的。

怎么做...

程序在本书代码包的boosting.ipynb文件中:

1. 进口情况如下:

   import ch9util
   from sklearn.grid_search import GridSearchCV
   from sklearn.ensemble import AdaBoostRegressor
   from sklearn.tree import DecisionTreeRegressor
   import numpy as np
   import dautil as dl
   from IPython.display import HTML
   

2. 加载数据并创建AdaBoostRegressor类:

   X_train, X_test, y_train, y_test = ch9util.temp_split()
   params = {
       'loss': ['linear', 'square', 'exponential'],
       'base_estimator__min_samples_leaf': [1, 2]
   }
   reg = AdaBoostRegressor(base_estimator=DecisionTreeRegressor(random_state=28),
                           random_state=17)
   

3. 网格搜索、拟合和预测如下:

   gscv = GridSearchCV(estimator=reg,
                       param_grid=params, cv=5, n_jobs=-1)
   gscv.fit(X_train, y_train)
   preds = gscv.predict(X_test)
   

4. 将预测值与实际值进行散点图:

   sp = dl.plotting.Subplotter(2, 2, context)
   html = ch9util.scatter_predictions(preds, y_test, gscv.best_params_,  
                                      gscv.best_score_, sp.ax)
   

5. 绘制一系列集合大小的验证曲线:

   nestimators = 2 ** np.arange(3, 9)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'n_estimators', nestimators)
   

6. 绘制一系列学习率的验证曲线:

   learn_rate = np.linspace(0.1, 1, 9)
   ch9util.plot_validation(sp.next_ax(), gscv.best_estimator_, 
                           X_train, y_train, 'learning_rate', learn_rate)
   

7. 绘制学习曲线如下:

   ch9util.plot_learn_curve(sp.next_ax(), gscv.best_estimator_, 
                            X_train, y_train)
   HTML(html + sp.exit())
   

有关最终结果，请参考以下屏幕截图:

另见

• 维基百科关于提升的页面在https://en . Wikipedia . org/wiki/Boosting _ % 28 machine _ learning % 29(2015 年 11 月检索)
• en.wikipedia.org/wiki/AdaBoo… 关于 AdaBoost 的维基百科页面(2015 年 11 月检索)
• 位于的AdaBoostRegressor课程的文档

嵌套交叉验证

如果我们将数据拟合到一条直线上，数学模型的参数将是直线的斜率和截距。当我们确定模型的参数时，我们在数据的子集(训练集)上拟合模型，并在其余数据(测试集)上评估模型的性能。这叫做 验证还有更精细的方案。例如，scikit-learn GridSearchCV类使用 k 倍交叉验证。

分类器和回归器通常需要额外的参数(超参数)，例如集成的组件数量，这通常与第一句中提到的线性模型无关。谈论模型有点令人困惑，因为我们有带有普通参数的模型和带有超参数的更大的模型。

让我们称较大的模型为 2 级模型，尽管据我所知这不是标准术语。如果我们使用GridSearchCV来获得 2 级模型的超参数，我们还有另一组参数(不是超参数或 1 级参数)需要担心——折叠的数量和用于比较的度量。评估指标还没有通过评审(参考第 10 章、评估分类器、回归器和聚类，但是比我们在本章中使用的指标更多。此外，我们可能会担心我们是否在用于评估结果的相同数据上确定了超参数。一种解决方案是应用 嵌套交叉验证。

嵌套交叉验证由以下交叉验证组成:

• 内部交叉验证进行超参数优化，例如使用网格搜索
• 外部交叉验证用于评估绩效和进行统计分析

在这个配方中，我们将查看以下分布:

• 所有分数的分布
• GridSearchCV报告的每个外部交叉验证迭代的最佳分数分布
• 每个折叠的平均分数分布
• GridSearchCV迭代中分数标准差的分布

怎么做...

代码在本书代码包的nested_cv.ipynb文件中:

1. 进口情况如下:

   from sklearn.grid_search import GridSearchCV
   from sklearn.cross_validation import ShuffleSplit
   from sklearn.cross_validation import cross_val_score
   import dautil as dl
   from sklearn.ensemble import ExtraTreesRegressor
   from joblib import Memory
   import numpy as np
   from IPython.display import HTML
   
   memory = Memory(cachedir='.')
   

2. 按照上一节

   @memory.cache
   def get_scores():
       df = dl.data.Weather.load()[['WIND_SPEED', 'TEMP', 'PRESSURE']].dropna()
       X = df.values[:-1]
       y = df['TEMP'][1:]
   
       params = { 'min_samples_split': [1, 3],
               'min_samples_leaf': [3, 4]}
   
       gscv = GridSearchCV(ExtraTreesRegressor(bootstrap=True,
                                               random_state=37),
                           param_grid=params, n_jobs=-1, cv=5)
       cv_outer = ShuffleSplit(len(X), n_iter=500,
                               test_size=0.3, random_state=55)
       r2 = []
       best = []
       means = []
       stds = []
   
       for train_indices, test_indices in cv_outer:
           train_i = X[train_indices], y[train_indices]
           gscv.fit(*train_i)
           test_i = X[test_indices]
           gscv.predict(test_i)
           grid_scores = dl.collect.flatten([g.cv_validation_scores
               for g in gscv.grid_scores_])
           r2.extend(grid_scores)
           means.extend(dl.collect.flatten([g.mean_validation_score
               for g in gscv.grid_scores_]))
           stds.append(np.std(grid_scores))
           best.append(gscv.best_score_)
   
       return {'r2': r2, 'best': best, 'mean': means, 'std': stds}
   

   的描述，得到 R 平方的分数

3. 获取分数并将其载入 NumPy 数组:

   scores = get_scores()
   r2 = np.array(scores['r2'])
   avgs = np.array(scores['mean'])
   stds = np.array(scores['std'])
   best = np.array(scores['best'])
   

4. 将分布绘制如下:

   sp = dl.plotting.Subplotter(2, 2, context)
   dl.plotting.hist_norm_pdf(sp.ax, r2)
   sp.label()
   
   dl.plotting.hist_norm_pdf(sp.next_ax(), best)
   sp.label()
   
   dl.plotting.hist_norm_pdf(sp.next_ax(), avgs)
   sp.label()
   
   dl.plotting.hist_norm_pdf(sp.next_ax(), stds)
   sp.label()
   HTML(sp.exit())
   

有关最终结果(交叉验证结果的分布)，请参考以下屏幕截图:

另见

• 关于交叉验证的维基百科页面
• 关于en.wikipedia.org/wiki/Hyperp…超参数优化的维基百科页面(2015 年 11 月检索)

使用 joblib 重用模型

joblib Memory类是一个实用程序类，便于将函数或方法结果缓存到磁盘。我们通过指定一个缓存目录来创建一个Memory对象。然后，我们可以修饰函数以进行缓存，或者在类构造函数中指定要缓存的方法。如果愿意，可以指定要忽略的参数。Memory类的默认行为是在函数修改或输入值改变时移除缓存。显然，您也可以通过移动或删除缓存目录和文件来手动删除缓存。

在这个食谱中，我描述了如何重用 scikit-learn 回归器或分类器。天真的方法是将对象存储在标准 Python 泡菜中或使用 joblib。然而，在大多数情况下，最好存储估计量的超参数。

我们将使用ExtraTreesRegressor类作为估计量。额外树 ( 极随机树)是随机森林算法的变种，在随机森林学习食谱中有所涉及。

怎么做...

1. 进口情况如下:

   from sklearn.grid_search import GridSearchCV
   from sklearn.ensemble import ExtraTreesRegressor
   import ch9util
   from tempfile import mkdtemp
   import os
   import joblib
   

2. 加载数据并定义超参数网格搜索字典:

   X_train, X_test, y_train, y_test = ch9util.temp_split()
   params = {'min_samples_split': [1, 3],
             'bootstrap': [True, False],
             'min_samples_leaf': [3, 4]}
   

3. 按照以下步骤进行网格搜索:

   gscv = GridSearchCV(ExtraTreesRegressor(random_state=41),
                       param_grid=params, cv=5)
   

4. 拟合和预测如下:

   gscv.fit(X_train, y_train)
   preds = gscv.predict(X_test)
   

5. 存储网格搜索找到的最佳参数:

   dir = mkdtemp()
   pkl = os.path.join(dir, 'params.pkl')
   joblib.dump(gscv.best_params_, pkl)
   params = joblib.load(pkl)
   print('Best params', gscv.best_params_)
   print('From pkl', params)
   

6. 创建一个新的估计值，并比较预测值:

   est = ExtraTreesRegressor(random_state=41)
   est.set_params(**params)
   est.fit(X_train, y_train)
   preds2 = est.predict(X_test)
   print('Max diff', (preds - preds2).max())
   

最终结果参见以下截图:

代码在本书代码包的reusing_models.py文件中。

另见

• pythonhosted.org/joblib/memo…的Memory班的文件(2015 年 11 月检索)
• 维基百科关于 en.wikipedia.org/wiki/Random… T2 随机森林的页面(2015 年 11 月检索)

分层聚类数据

在 Python 数据分析中，您学习了聚类——在不提供任何提示的情况下将数据分成聚类——这是无监督学习的一种形式。有时，我们需要猜测集群的数量，就像我们在中使用 Spark 方法对流式数据进行集群一样。

不限制簇包含其他簇。在这种情况下，我们谈到层次聚类。我们需要一个距离度量来分隔数据点。看看下面的等式:

在本食谱中，我们将使用 SciPy pdist()函数提供的欧几里德距离(9.2)。点集合之间的距离由链接标准给出。在本食谱中，我们将使用 SciPy linkage()功能提供的单连锁标准(9.3)。

怎么做...

脚本在本书代码包的clustering_hierarchy.ipynb文件中:

1. 进口情况如下:

   from scipy.spatial.distance import pdist
   from scipy.cluster.hierarchy import linkage
   from scipy.cluster.hierarchy import dendrogram
   import dautil as dl
   import matplotlib.pyplot as plt
   

2. 加载数据，重采样至年值，计算距离:

   df = dl.data.Weather.load().resample('A').dropna()
   dist = pdist(df)
   

3. 将分层聚类绘制如下:

   dendrogram(linkage(dist), labels=[d.year for d in df.index],
              orientation='right')
   plt.tick_params(labelsize=8)
   plt.xlabel('Cluster')
   plt.ylabel('Year')
   

有关最终结果，请参考以下屏幕截图:

另见

• en.wikipedia.org/wiki/Hierar… 的等级聚类的维基百科页面(2015 年 11 月检索)
• 位于的 pdist()功能文档
• 位于的 linkage()功能文档

进行一次茶道旅行

antao 是蒙特利尔的一个机器学习小组创建的 Python 库，通常与深度学习相关联，尽管这不一定是它的核心目的。antio 与 NumPy 紧密集成，可以在 CPU 或 GPU 上运行代码。如果您对图形处理器选项感兴趣，请参考部分列出的文档，另请参见部分。also 还支持通过符号变量进行符号微分。

根据它的文档，antao 是 NumPy 和 SymPy 的杂交。用 Anano 实现机器学习算法是可能的，但它没有使用 scikit-learn 那么容易和方便。但是，您可能会获得更高的并行性和数值稳定性的潜在优势。

在本食谱中，我们将使用梯度下降对温度数据进行线性回归。梯度下降是一种优化算法，我们可以在回归环境中使用它来最小化拟合残差。梯度衡量一个函数有多陡。为了找到一个局部极小值，该算法需要许多与梯度有多陡成正比的步骤。我们正试图走下坡路，但我们不知道在哪个方向可以找到局部最小值。所以，大幅度下跌平均会让我们下跌得更快，但不能保证。在某些情况下，它可能有助于平滑功能(更平滑的山丘)，因此我们不会花很多时间振荡。

做好准备

使用以下命令安装天线:

$ pip install --no-deps git+git://github.com/Theano/Theano.git 

截至 2015 年 11 月，我用出血边缘版本测试了代码。

怎么做...

代码在本书代码包的theano_tour.ipynb文件中:

1. 进口情况如下:

   import theano
   import numpy as np
   import theano.tensor as T
   import ch9util
   from sklearn.cross_validation import train_test_split
   from sklearn.metrics import r2_score
   import dautil as dl
   from IPython.display import HTML
   

2. 加载温度数据并定义无符号变量:

   temp = dl.data.Weather.load()['TEMP'].dropna()
   X = temp.values[:-1]
   y = temp.values[1:]
   X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=16)
   w = theano.shared(0., name ='w')
   c = theano.shared(0., name ='c')
   
   x = T.vector('x')
   y = T.vector('y')
   

3. 定义预测和成本(损失)函数以最小化:

   prediction = T.dot(x, w) + c
   cost = T.sum(T.pow(prediction - y, 2))/(2 * X_train.shape[0])
   Define gradient functions as follows:
   gw = T.grad(cost, w)
   gc = T.grad(cost, c)
   
   learning_rate = 0.01
   training_steps = 10000
   

4. 将训练函数定义如下:

   train = theano.function([x, y], cost, updates =
                           [(w, w - learning_rate * gw),
                            (c, c - learning_rate * gc)])
   predict = theano.function([x], prediction)
   

5. 按照以下步骤训练估算器:

   for i in range(training_steps):
       train(X_train.astype(np.float), y_train)
   

6. 预测和可视化预测如下:

   preds = predict(X_test)
   r2 = r2_score(preds, y_test)
   HTML(ch9util.scatter_predictions(preds, y_test, '', r2))
   

有关最终结果，请参考以下屏幕截图:

另见

• deeplearning.net/software/th… 年 11 月检索)
• 维基百科关于en.wikipedia.org/wiki/Gradie…梯度下降的页面(2015 年 11 月检索)