Python-数据分析秘籍-二-Python 数据分析秘籍（二）四、处理数据和数字问题本章中的配方如下: 裁剪和过滤

Python 数据分析秘籍（二）

原文：Mastering Python Data Analysis

协议：CC BY-NC-SA 4.0

四、处理数据和数字问题

本章中的配方如下:

裁剪和过滤异常值
正在获取数据
测量噪声数据的中心趋势
用 Box-Cox 变换归一化
用权力阶梯转换数据
用对数转换数据
重新绑定数据
应用logit()变换比例
拟合稳健的线性模型
用加权最小二乘法考虑方差
使用任意精度进行优化
线性代数使用任意精度

简介

在现实世界中，数据很少与教科书的定义和例子相匹配。我们必须处理诸如硬件故障、不合作的客户和不满的同事等问题。很难预测你会遇到什么样的问题，但可以肯定的是，这些问题会非常丰富和具有挑战性。在本章中，我将概述一些处理有噪声数据的常见方法，这些方法更多地基于经验法则，而不是严格的科学。幸运的是，数据分析的试错部分是有限的。

本章的大部分内容是关于异常值管理的。离群值是我们认为不正常的值。当然，这不是你唯一会遇到的问题，但这是一个偷偷摸摸的问题。一个常见的问题是缺少或无效的值，所以我将简要地提到屏蔽数组和 pandas 特性，例如我在本书中使用的dropna()函数。

我还写了两个关于使用 mpmath 进行任意精度计算的食谱。我不建议使用 mpmath，除非你真的必须这样做，因为你必须付出性能代价。通常我们可以解决数值问题，所以很少需要任意精度的库。

裁剪和过滤异常值

异常值是数据分析中常见的问题。虽然不存在异常值的精确定义，但我们知道异常值会影响均值和回归结果。异常值是异常的值。通常，异常值是由测量误差引起的，但异常值有时是真实的。在第二种情况下，我们可能要处理两种或两种以上与不同现象相关的数据。

该配方的数据在中描述。它由某一星团中 47 颗恒星的对数有效温度和对数光强组成。任何一个天文学家读了这一段都会知道 T4【赫茨普鲁-罗素图。从数据分析的角度来看，该图是一个散点图，但对于天文学家来说，它当然不止于此。赫茨普龙罗素图是在 1910 年左右定义的，其特征是一条对角线(不完全是直线)，称为主序列。我们数据集中的大多数恒星应该在主序列上，左上角有四个异常值。这些异常值被归类为巨人。

我们有许多策略来处理异常值。在这个食谱中，我们将使用两个最简单的策略:用 NumPy clip()函数裁剪和完全去除异常值。对于本例，我将异常值定义为从第一个和第三个四分位数定义的框中移除的 1.5 个四分位数区间的值。

怎么做...

以下步骤显示了如何裁剪和过滤异常值:

进口情况如下:

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

定义以下函数过滤异常值:

def filter_outliers(a):
    b = a.copy()
    bmin, bmax = dl.stats.outliers(b)
    b[bmin > b] = np.nan
    b[b > bmax] = np.nan

    return b

如下所示加载和裁剪异常值:

starsCYG = sm.datasets.get_rdataset("starsCYG", "robustbase", cache=True).data

clipped = starsCYG.apply(dl.stats.clip_outliers)

过滤异常值如下:

filtered = starsCYG.copy()
filtered['log.Te'] = filter_outliers(filtered['log.Te'].values)
filtered['log.light'] = filter_outliers(filtered['log.light'].values)
filtered.dropna()

用以下代码绘制结果:

sp = dl.plotting.Subplotter(3, 1, context)
sp.label()
sns.regplot(x='log.Te', y='log.light', data=starsCYG, ax=sp.ax)
sp.label(advance=True)
sns.regplot(x='log.Te', y='log.light', data=clipped, ax=sp.ax)
sp.label(advance=True)
sns.regplot(x='log.Te', y='log.light', data=filtered, ax=sp.ax)
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见下面的截图(参见本书代码包中的outliers.ipynb文件):

How to do it...

另见

https://docs . scipy . org/doc/NumPy/reference/generated/NumPy . clip . html # NumPy . clip记录的 NumPy clip()功能(2015 年 8 月检索)
您可以在https://en . Wikipedia . org/wiki/赫茨普龙%E2%80%93Russell_diagram 上阅读更多关于 Hertzsprung-Russell 图的信息(检索于 2015 年 8 月)

winsoring 数据

winsoring是处理异常值的另一种技术，以 Charles Winsor 命名。实际上，Winsorization 以对称的方式将异常值裁剪到给定的百分位数。例如，我们可以剪辑到第 5 和第 95 百分位。SciPy 有一个winsorize()功能，执行这个程序。该配方的数据与裁剪和过滤异常值配方的数据相同。

怎么做...

使用以下程序读取数据:

进口情况如下:

rom scipy.stats.mstats import winsorize
import statsmodels.api as sm
import seaborn as sns
import matplotlib.pyplot as plt
import dautil as dl
from IPython.display import HTML

加载并保存有效温度的数据(限制设置为 15%):

starsCYG = sm.datasets.get_rdataset("starsCYG", "robustbase", cache=True).data
limit = 0.15
winsorized_x = starsCYG.copy()
winsorized_x['log.Te'] = winsorize(starsCYG['log.Te'], limits=limit)

按照以下方式调节光强:

winsorized_y = starsCYG.copy()
winsorized_y['log.light'] = winsorize(starsCYG['log.light'], limits=limit)
winsorized_xy = starsCYG.apply(winsorize, limits=[limit, limit])

用回归线绘制赫茨普鲁格-罗素图(不是通常天文图的一部分):

sp = dl.plotting.Subplotter(2, 2, context)
sp.label()
sns.regplot(x='log.Te', y='log.light', data=starsCYG, ax=sp.ax)

sp.label(advance=True)
sns.regplot(x='log.Te', y='log.light', data=winsorized_x, ax=sp.ax)

sp.label(advance=True)
sns.regplot(x='log.Te', y='log.light', data=winsorized_y, ax=sp.ax)

sp.label(advance=True)
sns.regplot(x='log.Te', y='log.light', data=winsorized_xy, ax=sp.ax)
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果见下图截图(参见本书代码包中的winsorising_data.ipynb文件):

How to do it...

另见

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

测量噪声数据的中心趋势

我们可以用平均值和中位数来测量中心趋势。这些措施使用所有可用的数据。通过丢弃数据集中较高端和较低端的数据来消除异常值是一个普遍接受的想法。截断表示 或 修剪表示，其派生词如 四分位数间平均值 ( IQM )和曲美、也采用这种思路。看看下面的等式:

Measuring central tendency of noisy data

截断平均值在给定的百分位数丢弃数据，例如，从最低值到第 5 个百分位数，从第 95 个百分位数到最高值。曲美安 (4.1)是中位数、第一四分位数和第三四分位数的加权平均值。对于 IQM (4.2)，我们丢弃数据的最低和最高四分位数，因此这是截断平均值的特殊情况。我们将使用 SciPy tmean()和trima()函数计算这些度量。

怎么做...

我们将通过以下步骤研究不同截断水平的主要趋势:

进口情况如下:

import matplotlib.pyplot as plt
from scipy.stats import tmean
from scipy.stats.mstats import trima
import numpy as np
import dautil as dl
import seaborn as sns
from IPython.display import HTML

context = dl.nb.Context('central_tendency')

定义以下函数计算四分位数平均值:

def iqm(a):
    return truncated_mean(a, 25)

定义以下函数来绘制分布:

def plotdists(var, ax):
    displot_label = 'From {0} to {1} percentiles'
    cyc = dl.plotting.Cycler()

    for i in range(1, 9, 3):
        limits = dl.stats.outliers(var, method='percentiles', 
                               percentiles=(i, 100 - i))
        truncated = trima(var, limits=limits).compressed()
        sns.distplot(truncated, ax=ax, color=cyc.color(),
                     hist_kws={'histtype': 'stepfilled', 'alpha': 1/i,
                               'linewidth': cyc.lw()},
                     label=displot_label.format(i, 100 - i))

定义以下函数计算截断平均值:

def truncated_mean(a, percentile):
    limits = dl.stats.outliers(a, method='percentiles', 
                               percentiles=(percentile, 100 - percentile))

    return tmean(a, limits=limits)

加载数据，计算方式如下:

df = dl.data.Weather.load().resample('M').dropna()
x = range(9)
temp_means = [truncated_mean(df['TEMP'], i) for i in x]
ws_means = [truncated_mean(df['WIND_SPEED'], i) for i in x]

用以下代码绘制平均值和分布:

sp = dl.plotting.Subplotter(2, 2, context)
cp = dl.plotting.CyclePlotter(sp.ax)
cp.plot(x, temp_means, label='Truncated mean')
cp.plot(x, dl.stats.trimean(df['TEMP']) * np.ones_like(x), label='Trimean')
cp.plot(x, iqm(df['TEMP']) * np.ones_like(x), label='IQM')
sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(x, ws_means, label='Truncated mean')
cp.plot(x, dl.stats.trimean(df['WIND_SPEED']) * np.ones_like(x),
        label='Trimean')
cp.plot(x, iqm(df['WIND_SPEED']) * np.ones_like(x), label='IQM')
sp.label(ylabel_params=dl.data.Weather.get_header('WIND_SPEED'))

plotdists(df['TEMP'], sp.next_ax())
sp.label(xlabel_params=dl.data.Weather.get_header('TEMP'))

plotdists(df['WIND_SPEED'], sp.next_ax())
sp.label(xlabel_params=dl.data.Weather.get_header('WIND_SPEED'))
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见以下截图(参见本书代码包中的central_tendency.ipynb文件):

How to do it...

另见

位于的trima()的 SciPy 文档(检索日期:2015 年 8 月)
位于的tmean()的 SciPy 文档(检索日期:2015 年 8 月)

用 Box-Cox 变换归一化

不遵循已知分布(如正态分布)的数据通常难以管理。控制数据的一个流行策略是应用 Box-Cox 变换。它由以下等式给出:

Normalizing with the Box-Cox transformation

scipy.stats.boxcox()函数可以对正数据应用转换。我们将使用与裁剪和过滤异常值配方中相同的数据。通过 Q-Q 图，我们将表明 Box-Cox 变换确实使数据看起来更正常。

怎么做...

以下步骤显示了如何使用 Box-Cox 转换对数据进行规范化:

进口情况如下:

import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy.stats import boxcox
import seaborn as sns
import dautil as dl
from IPython.display import HTML

加载数据并转换如下:

context = dl.nb.Context('normalizing_boxcox')

starsCYG = sm.datasets.get_rdataset("starsCYG", "robustbase", cache=True).data

var = 'log.Te'

# Data must be positive
transformed, _ = boxcox(starsCYG[var])

显示 Q-Q 图和分布图如下:

sp = dl.plotting.Subplotter(2, 2, context)
sp.label()
sm.qqplot(starsCYG[var], fit=True, line='s', ax=sp.ax)

sp.label(advance=True)
sm.qqplot(transformed, fit=True, line='s', ax=sp.ax)

sp.label(advance=True)
sns.distplot(starsCYG[var], ax=sp.ax)

sp.label(advance=True)
sns.distplot(transformed, ax=sp.ax)                                       
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见下面的截图(参见本书代码包中的normalizing_boxcox.ipynb文件):

How to do it...

它是如何工作的

在之前的截图中，Q-Q 图将正态分布的理论分位数与实际数据的分位数进行了对比。为了帮助评估是否符合正态分布，我显示了一条线，它应该与完全正常的数据一致。数据越符合直线，就越正常。如您所见，转换后的数据更符合直线，因此更正常。分布图应该有助于你确认这一点。

另见

相关维基百科页面在en.wikipedia.org/wiki/Power_…(2015 年 8 月检索)
《转换分析》，英国皇家统计学会杂志，第 26 期，第 211-252 页(1964 年)。

用权力阶梯转换数据

线性关系在科学和数据分析中很常见。显然，线性模型比非线性模型更容易理解。所以历史上，线性模型的工具最先被开发出来。在某些情况下，线性化(使数据线性)以简化分析是有好处的。一个简单的策略有时会奏效，那就是将一个或多个变量平方或立方。同样，我们可以通过取平方根或立方根，将数据转换成一个假想的幂阶梯。

在本食谱中，我们将使用中描述的邓肯数据集的数据。这些数据是在 1961 年前后收集的，大约有 45 种职业，分为四栏——类型、收入、教育和声望。我们来看看收入和声望。这些变量似乎通过三次多项式联系在一起，所以我们可以取收入的立方根或声望的立方根。为了检查结果，我们将把回归的残差可视化。期望是残差是随机分布的，这意味着我们不期望它们遵循可识别的模式。

怎么做...

在以下步骤中，我将演示基本的数据转换:

进口情况如下:

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

加载并转换数据如下:

df = sm.datasets.get_rdataset("Duncan", "car", cache=True).data
transformed = df.copy()
transformed['income'] = np.power(transformed['income'], 1.0/3)

用 Seaborn 回归图(三次多项式)绘制原始数据，如下所示:

sp = dl.plotting.Subplotter(2, 2, context)
sp.label()
sns.regplot(x='income', y='prestige', data=df, order=3, ax=sp.ax)

用以下线条绘制转换后的数据:

sp.label(advance=True)
sns.regplot(x='income', y='prestige', data=transformed, ax=sp.ax)

绘制三次多项式的残差图:

sp.label(advance=True)
sns.residplot(x='income', y='prestige', data=df, order=3, ax=sp.ax)

绘制残差图绘制转换数据如下:

sp.label(advance=True)
sp.ax.set_xlim([1, 5])
sns.residplot(x='income', y='prestige', data=transformed, ax=sp.ax)
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见以下截图(参见本书代码包中的transforming_up.ipynb文件):

How to do it...

用对数转换数据

当数据按数量级变化时，用对数转换数据是一个显而易见的策略。在我的经验中，用指数函数做相反的变换不太常见。通常在探索时，我们将成对变量的对数或半对数散点图可视化。

为了展示这一转变，我们将使用世界银行的数据，计算每 1000 例活产婴儿死亡率和可用国家的人均国内生产总值(T2)(T4)。如果我们对两个变量都应用以 10 为底的对数，我们通过拟合数据得到的直线的斜率有一个有用的性质。一个变量增加 1%对应于另一个变量斜率给出的百分比变化。

怎么做...

按照以下步骤使用对数转换数据:

进口情况如下:

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

下载 2010 年数据，代码如下:

wb = dl.data.Worldbank()
countries = wb.get_countries()[['name', 'iso2c']]
inf_mort = wb.get_name('inf_mort')
gdp_pcap = wb.get_name('gdp_pcap')
df = wb.download(country=countries['iso2c'],
                 indicator=[inf_mort, gdp_pcap],
                 start=2010, end=2010).dropna()

使用以下代码片段应用日志转换:

loglog = df.applymap(np.log10)
x = loglog[gdp_pcap]
y = loglog[inf_mort]

绘制变换前后的数据:

sp = dl.plotting.Subplotter(2, 1, context)
xvar = 'GDP per capita'
sp.label(xlabel_params=xvar)
sp.ax.set_ylim([0, 200])
sp.ax.scatter(df[gdp_pcap], df[inf_mort])

sp.next_ax()
sp.ax.scatter(x, y, label='Transformed')
dl.plotting.plot_polyfit(sp.ax, x, y)
sp.label(xlabel_params=xvar)
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见下面的截图(参见本书代码包中的transforming_down.ipynb文件):

How to do it...

重新绑定数据

通常，我们拥有的数据并不是按照我们想要的方式构建的。我们可以使用的一种结构化技术叫做(统计)数据宁滨或逆势。该策略用一个代表值替换一个区间(一个箱)内的值。在这个过程中，我们可能会丢失信息；但是，我们可以更好地控制数据和效率。

在天气数据集中，我们有以度为单位的风向和以米/秒为单位的风速，它们可以用不同的方式表示。在这个食谱中，我选择用基本方向(北、南等)来表示风向。对于风速，我使用了博福特标度(访问en.wikipedia.org/wiki/Beaufo…)。

怎么做...

按照以下说明重新绑定数据:

进口情况如下:

import dautil as dl
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from IPython.display import HTML

按照以下方式加载和重新绑定数据(风向为 0-360 度；我们重新绑定到基本方向，如北、西南等):

df = dl.data.Weather.load()[['WIND_SPEED', 'WIND_DIR']].dropna()
categorized = df.copy()
categorized['WIND_DIR'] = dl.data.Weather.categorize_wind_dir(df)
categorized['WIND_SPEED'] = dl.data.Weather.beaufort_scale(df)

用以下代码显示分布和计数图:

sp = dl.plotting.Subplotter(2, 2, context)
sns.distplot(df['WIND_SPEED'], ax=sp.ax)
sp.label(xlabel_params=dl.data.Weather.get_header('WIND_SPEED'))

sns.distplot(df['WIND_DIR'], ax=sp.next_ax())
sp.label(xlabel_params=dl.data.Weather.get_header('WIND_DIR'))

sns.countplot(x='WIND_SPEED', data=categorized, ax=sp.next_ax())
sp.label()

sns.countplot(x='WIND_DIR', data=categorized, ax=sp.next_ax())
sp.label()
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见下面的截图(参见本书代码包中的rebinning_data.ipynb文件):

How to do it...

应用 logit()变换比例

我们可以使用 SciPy logit()函数转换比例或比率。结果应该是更高斯的分布。该函数是由以下等式定义的:

Applying logit() to transform proportions

正如你在等式(4.4)中看到的，logit 是赔率的对数。我们希望通过这种转换获得更对称的分布——接近于零的偏斜。当比例接近 0 和 1 时，logit 渐近地接近负无穷大和无穷大，所以在这些情况下我们必须小心。

作为一个比例的例子，我们将采取每月的雨天比例。我们通过将降雨量转化为二进制变量，然后对每个月进行平均，得到这些比例。

怎么做...

按照本指南通过转换比率:

进口情况如下:

import dautil as dl
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import math
import statsmodels.api as sm
from scipy.special import logit
from IPython.display import HTML

加载数据并用以下代码转换:

rain = dl.data.Weather.load()['RAIN'].dropna()
rain = rain > 0
rain = rain.resample('M').dropna()
transformed = rain.apply(logit)
transformed = dl.data.dropinf(transformed.values)

用分布图和 Q-Q 图绘制转换的结果:

sp = dl.plotting.Subplotter(2, 2, context)
sns.distplot(rain, ax=sp.ax)
sp.label()

sp.label(advance=True)
sns.distplot(transformed, ax=sp.ax)

sp.label(advance=True)
sm.qqplot(rain, line='s', ax=sp.ax)

sp.label(advance=True)
sm.qqplot(transformed, line='s', ax=sp.ax)
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果见下图截图(参见本书代码包中的transforming_ratios.ipynb文件):

How to do it...

拟合稳健线性模型

稳健回归旨在比普通回归更好地处理数据中的异常值。这种类型的回归使用特殊的稳健估计量，这也受到状态模型的支持。显然，没有最佳估计量，所以估计量的选择取决于数据和模型。

在这个配方中，我们将拟合 statsmodels 中可用的年太阳黑子数数据。我们将定义一个简单的模型，其中当前计数线性依赖于先前的值。为了演示异常值的影响，我添加了一个相当大的值，我们将比较稳健回归模型和普通最小二乘模型。

怎么做...

以下步骤描述了如何应用稳健线性模型:

进口情况如下:

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

定义以下功能来设置地块的标签:

def set_labels(ax):
    ax.set_xlabel('Year')
    ax.set_ylabel('Sunactivity')

定义以下功能来绘制模型拟合:

def plot_fit(df, ax, results):
    x = df['YEAR']
    cp = dl.plotting.CyclePlotter(ax)
    cp.plot(x[1:], df['SUNACTIVITY'][1:], label='Data')
    cp.plot(x[2:], results.predict()[1:], label='Fit')
    ax.legend(loc='best')

加载数据并添加一个异常值用于演示目的:

df = sm.datasets.sunspots.load_pandas().data
vals = df['SUNACTIVITY'].values

# Outlier added by malicious person, because noone
# laughs at his jokes.
vals[0] = 100

如下所示拟合稳健模型:

rlm_model = sm.RLM(vals[1:], sm.add_constant(vals[:-1]),
                   M=sm.robust.norms.TrimmedMean())

rlm_results = rlm_model.fit()
hb = dl.report.HTMLBuilder()
hb.h1('Fitting a robust linear model')
hb.h2('Robust Linear Model')
hb.add(rlm_results.summary().tables[1].as_html())

拟合一个普通的最小二乘模型:

hb.h2('Ordinary Linear Model')
ols_model = sm.OLS(vals[1:], sm.add_constant(vals[:-1]))
ols_results = ols_model.fit()
hb.add(ols_results.summary().tables[1].as_html())

用以下代码绘制数据和模型结果:

fig, [ax, ax2] = plt.subplots(2, 1)

plot_fit(df, ax, rlm_results)
ax.set_title('Robust Linear Model')
set_labels(ax)

ax2.set_title('Ordinary Least Squares')
plot_fit(df, ax2, ols_results)
set_labels(ax2)
plt.tight_layout()
HTML(hb.html)

最终结果参见下面的截图(参见本书代码包中的rlm_demo.ipynb文件):

How to do it...

另见

相关的 statsmodels 文档位于http://stats models . source forge . net/0 . 6 . 0/generated/stats models . robust . robust _ linear _ model。RLM.html(2015 年 8 月检索)

用加权最小二乘法考虑方差

statsmodels 库允许我们为回归定义每个数据点的任意权重。用简单的拇指规则有时很容易发现异常值。其中一个经验法则是基于四分位数范围，即第一个四分位数和第三个四分位数数据之间的差异。通过四分位数范围，我们可以定义加权最小二乘回归的权重。

我们将使用来自的数据和模型来拟合鲁棒的线性模式，但是具有任意的权重。我们怀疑是异常值的点将获得较低的权重，这是刚才提到的四分位数范围值的倒数。

怎么做...

使用以下方法用加权最小二乘法拟合数据:

进口情况如下:

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

加载数据并添加一个异常值:

temp = dl.data.Weather.load()['TEMP'].dropna()
temp = dl.ts.groupby_yday(temp).mean()

# Outlier added by malicious person, because noone
# laughs at his jokes.
temp.values[0] = 100

使用普通最小二乘模型拟合:

ntemp = len(temp)
x = np.arange(1, ntemp + 1)
factor = 2 * np.pi/365.25
cos_x = sm.add_constant(np.cos(-factor * x - factor * 337))
ols_model = sm.OLS(temp, cos_x)
ols_results = ols_model.fit()
hb = dl.report.HTMLBuilder()
hb.h1('Taking variance into account with weighted least squares')
hb.h2('Ordinary least squares')
hb.add(ols_results.summary().tables[1].as_html())
ols_preds = ols_results.predict()

使用四分位数范围计算权重，并拟合加权最小二乘模型:

box = dl.stats.Box(temp)
iqrs = box.iqr_from_box()
# Adding 1 to avoid div by 0
weights = 1./(iqrs + 1)
wls_model = sm.WLS(temp, cos_x, weights=weights)
wls_results = wls_model.fit()

hb.h2('Weighted least squares')
hb.add(wls_results.summary().tables[1].as_html())

绘制模型结果和权重:

sp = dl.plotting.Subplotter(2, 2, context)

sp.ax.plot(x[1:], temp[1:], 'o', label='Data')
sp.ax.plot(x[1:], ols_preds[1:], label='Fit')
sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

sp.label(advance=True)
sp.ax.plot(x, iqrs, 'o')

sp.next_ax().plot(x[1:], temp[1:], 'o', label='Data')
sp.ax.plot(x[1:], wls_results.predict()[1:], label='Fit')
sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

sp.label(advance=True)
sp.ax.plot(x, weights, 'o')
plt.tight_layout()
HTML(hb.html)

最终结果参见下面的截图(参见本书代码包中的weighted_ls.ipynb文件):

How to do it...

另见

相关 statsmodels 文档位于http://stats models . SourceForge . net/0 . 5 . 0/examples/generated/example _ WLS . html(2015 年 8 月检索)。

使用任意精度进行优化

这本书的目标读者应该知道浮点数的问题。我会提醒你，我们不能准确地表示浮点数。偶数整数表示是有限的。对于某些应用，例如金融计算或涉及已知解析表达式的工作，我们需要比 NumPy 等数值软件更高的精度。Python 标准库提供了Decimal类，我们可以用它来实现任意精度。然而，专门的mpmath库更适合更高级的用途。

温度遵循季节性模式，因此包含余弦的模型似乎是自然的。我们将应用这样的模型。使用任意精度的好处是，您可以轻松地进行分析、微分、求根和逼近多项式。

做好准备

使用以下任一命令安装专用mpmath库:

$ conda install mpmath
$ pip install mpmath

我通过 Anaconda 用 mpmath 0.19 测试了代码。

怎么做...

以下说明描述了如何使用任意精度进行优化:

进口情况如下:

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

为模型和一阶导数定义以下函数:

def model(t):
    mu, C, w, phi = (9.6848106, -7.59870042, -0.01766333, -5.83349705)

    return mu + C * mpmath.cos(w * t + phi)

def diff_model(t):
    return mpmath.diff(model, t)

加载数据，求一阶导数的根:

vals = dl.data.Weather.load()['TEMP'].dropna()
vals = dl.ts.groupby_yday(vals).mean()
diff_root = mpmath.findroot(diff_model, (1, 366), solver='anderson')

获得模型的多项式近似，如下所示:

days = range(1, 367)
poly = mpmath.chebyfit(model, (1, 366), 3)
poly = np.array([float(c) for c in poly])

用以下代码绘制数据、建模结果和近似值:

sp = dl.plotting.Subplotter(2, 1, context)
cp = dl.plotting.CyclePlotter(sp.ax)
cp.plot(days, [model(i) for i in days], label='Model')
cp.plot(days, vals, label='Data')
sp.ax.annotate(s='Root of derivative', xy=(diff_root, vals.max() - 1),
            xytext=(diff_root, vals.max() - 8),
            arrowprops=dict(arrowstyle='->'))
yvar = dl.data.Weather.get_header('TEMP')
sp.label(ylabel_params=yvar)

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(days, vals, label='Data')
cp.plot(days, np.polyval(poly, days), label='Approximation')
sp.label(ylabel_params=yvar)
plt.tight_layout()
HTML(dl.report.HTMLBuilder().watermark())

最终结果参见下面的截图(参见本书代码包中的mpmath_fit.ipynb文件):

How to do it...

另见

位于的chebyfit()函数的文档(检索于 2015 年 8 月)
位于的 findroot()函数的文档

对线性代数使用任意精度

很多模型可以归结为线性方程组，其中是线性代数的定义域。中提到的使用任意精度进行优化的 mpmath 库配方也可以做任意精度的线性代数。

理论上，我们可以将任何可微函数近似为多项式级数。为了找到多项式的系数，我们可以定义一个线性方程组，基本上是取一个数据向量(向量作为数学项)的幂，用一个 1 的向量来表示多项式中的常数。我们将用 mpmath lu_solve()函数来求解这样一个系统。作为示例数据，我们将使用按一年中的某一天分组的风速数据。

做好准备

有关说明，请参考使用任意精度进行优化配方。

怎么做...

按照以下步骤对线性代数使用任意精度:

进口情况如下:

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

定义以下函数，用 mpmath 计算算术平均值:

def mpmean(arr):
    mpfs = [mpmath.mpf(a) for a in arr]

    return sum(mpfs)/len(arr)

加载数据，用lu_solve() :

vals = dl.data.Weather.load()['WIND_SPEED'].dropna()
vals = dl.ts.groupby_yday(vals).apply(mpmean)

days = np.arange(1, 367, dtype=int)
A = [[], [], []]
A[0] = np.ones_like(days, dtype=int).tolist()
A[1] = days.tolist()
A[2] = (days ** 2).tolist()
A = mpmath.matrix(A).transpose()

params = mpmath.lu_solve(A, vals)

result = dl.report.HTMLBuilder()
result.h1('Arbitrary Precision Linear Algebra')
result.h2('Polynomial fit')
dfb = dl.report.DFBuilder(['Coefficient 0', 'Coefficient 1', 'Coefficient 2'])
dfb.row(params)
result.add_df(dfb.build())

求解系统

定义以下函数来计算我们得到的多项式:

def poly(x):
    return mpmath.polyval(params[::-1], x)

使用fourier()函数获得三角近似:

cs = mpmath.fourier(poly, days.tolist(), 1)
result.h2('Cosine and sine terms')
dfb = dl.report.DFBuilder(['Coefficient 1', 'Coefficient 2'])
dfb.row(cs[0])
dfb.row(cs[1])
result.add_df(dfb.build(index=['Cosine', 'Sine']), index=True)

将数据、模型结果和近似值绘制如下:

sp = dl.plotting.Subplotter(2, 1, context)

cp = dl.plotting.CyclePlotter(sp.ax)
cp.plot(days, vals, label='Data')
cp.plot(days, poly(days), label='Fit')
yvar = dl.data.Weather.get_header('WIND_SPEED')
sp.label(ylabel_params=yvar)

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(days, vals, label='Data')
cp.plot(days, [mpmath.fourierval(cs, days, d) for d in days], label='Approximation')
sp.label(ylabel_params=yvar)
plt.tight_layout()
HTML(result.html)

最终结果参见下面的截图(参见本书代码包中的mpmath_linalg.ipynb文件):

How to do it...

另见

在https://mpmath . readed docs . org/en/latest/calculation/近似值. html 中记录的fourier()函数？高光=傅立叶#mpmath .傅立叶(检索于 2015 年 8 月)
lu_solve()功能在mpmath.readthedocs.org/en/latest/m…(2015 年 8 月检索)

五、网络挖掘、数据库和大数据

本章菜单上有以下食谱:

模拟网络浏览
刮网
处理非 ASCII 文本和 HTML 实体
实现关联表
设置数据库迁移脚本
向现有表中添加表列
创建表后添加索引
设置测试 web 服务器
用事实表和维度表实现星型模式
利用 HDFS
设置火花
使用 Spark 对数据进行聚类

简介

这一章侧重于数学，但更侧重于技术主题。科技可以为数据分析师提供很多东西。数据库已经存在了一段时间，但是大多数人熟悉的关系数据库可以追溯到 20 世纪 70 年代。埃德加·科德提出了许多想法，这些想法后来导致了关系模型和 SQL 的创建。从那以后，关系数据库一直是一项主导技术。在 20 世纪 80 年代，面向对象的编程语言引起了范式的转变，并不幸与关系数据库不匹配。

面向对象编程语言支持继承等概念，而关系数据库和 SQL 不支持这些概念(当然也有一些例外)。Python 生态系统有几个试图解决这种不匹配问题的对象关系映射 ( ORM )框架。这是不可能的，也没有必要全部覆盖，所以我选择了 SQLAlchemy 作为这里的食谱。我们还将把数据库模式迁移视为一个常见的热门话题，尤其是对于生产系统。

大数据是你可能听说过的流行语之一。Hadoop 和 Spark 可能听起来也很熟悉。我们将在本章中研究这些框架。如果您使用我的 Docker 映像，您将很不幸地在其中找不到 Selenium、Hadoop 和 Spark，因为为了节省空间，我决定不包含它们。

另一个重要的技术发展是万维网，也称为互联网。互联网是最终的数据来源；然而，以易于分析的形式获得这些数据有时是一个相当大的挑战。作为最后的资源，我们可能不得不抓取网页。成功没有保证，因为网站所有者可以在不警告我们的情况下更改内容。这是由你来保持网页抓取食谱的最新代码。

模拟网页浏览

公司网站通常由团队或部门使用专门的工具和模板制作。很多内容是动态生成的，由很大一部分 JavaScript 和 CSS 组成。这意味着即使我们下载了内容，我们仍然至少要评估 JavaScript 代码。我们可以通过 Python 程序做到这一点的一种方法是使用 **【硒】**应用编程接口。Selenium 的主要目的实际上是测试网站，但没有什么能阻止我们用它来刮网站。

我们将刮除本书代码包中的test_widget.ipynb文件，而不是刮除网站。为了模拟浏览这个网页，我们在test_simulating_browsing.py中提供了一个单元测试类。如果你想知道，这不是测试 IPython 笔记本的推荐方法。

出于历史原因，我更喜欢使用 XPath 来查找 HTML 元素。XPath 是一种查询语言，也适用于 HTML。这不是唯一的方法，您还可以使用 CSS 选择器、标签名或 id。为了找到正确的 XPath 表达式，你可以为你最喜欢的浏览器安装一个相关的插件，或者在谷歌浏览器中，你可以检查一个元素的 XPath。

做好准备

使用以下命令安装硒:

$ pip install selenium

我用 Selenium 2.47.1 测试了代码。

怎么做…

以下步骤向您展示了如何使用我制作的 IPython 小部件模拟网页浏览。该配方的代码在本书代码包的test_simulating_browsing.py文件中:

第一步是运行以下内容:
```
$ ipython notebook
```

进口情况如下:

from selenium import webdriver
import time
import unittest
import dautil as dl

NAP_SECS = 10

定义以下函数，创建一个火狐浏览器实例:

class SeleniumTest(unittest.TestCase):
    def setUp(self):
        self.logger = dl.log_api.conf_logger(__name__)
        self.browser = webdriver.Firefox()

定义测试完成后要清理的以下功能:

    def tearDown(self):
        self.browser.quit()

以下功能点击小部件标签(我们必须等待用户界面响应):

def wait_and_click(self, toggle, text):
        xpath = "//a[@data-toggle='{0}' and contains(text(), '{1}')]"
        xpath = xpath.format(toggle, text)
        elem = dl.web.wait_browser(self.browser, xpath)
        elem.click()

定义以下函数，该函数执行测试，包括评估笔记本单元格并单击 IPython 小部件中的几个选项卡(我们使用端口 8888):

    def test_widget(self):
        self.browser.implicitly_wait(NAP_SECS)
        self.browser.get('http://localhost:8888/notebooks/test_widget.ipynb')

        try:
            # Cell menu
            xpath = '//*[@id="menus"]/div/div/ul/li[5]/a'
            link = dl.web.wait_browser(self.browser, xpath)
            link.click()
            time.sleep(1)

            # Run all
            xpath = '//*[@id="run_all_cells"]/a'
            link = dl.web.wait_browser(self.browser, xpath)
            link.click()
            time.sleep(1)

            self.wait_and_click('tab', 'Figure')
            self.wait_and_click('collapse', 'figure.figsize')
        except Exception:
            self.logger.warning('Error while waiting to click', exc_info=True)
            self.browser.quit()

        time.sleep(NAP_SECS)
        self.browser.save_screenshot('widgets_screenshot.png')

if __name__ == "__main__":
    unittest.main()

以下截图由代码创建:

How to do it…

另见

Selenium 文档位于https://Selenium-python . readed docs . org/en/latest/installation . html(2015 年 9 月检索)
维基百科中关于 XPath 的页面位于en.wikipedia.org/wiki/XPath(检索于 2015 年 9 月)

刮网

我们知道搜索引擎发出被称为机器人的自主程序，在互联网上寻找信息。通常，这导致创建类似于电话簿或字典的巨大索引。Python 3 用户目前的情况(2015 年 9 月)在刮网方面并不理想。大多数框架只支持 Python 2。然而，圭多·范罗森斯，一生的仁慈独裁者 ( BDFL )刚刚在 GitHub 上贡献了一个使用 AsyncIO API 的爬虫。向 BDFL 致敬！

为了保存抓取的网址，我分叉了存储库并做了一些小的修改。我也让爬虫提前退出了。这些变化不是很优雅，但这是我在有限的时间内所能做的。无论如何，我不能指望比 BDFL 本人做得更好。

一旦我们有了网页链接列表，我们将从 Selenium 加载这些网页(参考模拟网页浏览食谱)。我选择了 PhantomJS，一款无头浏览器，应该比 Firefox 占用空间更小。虽然这不是绝对必要的，但我认为有时下载你正在抓取的网页是有意义的，因为你可以在本地测试抓取。您还可以更改下载的 HTML 中的链接，以指向本地文件。这与设置测试网络服务器配方有关。抓取的一个常见用例是创建一个用于语言分析的文本语料库。这是我们在这个食谱中的目标。

做好准备

按照模拟网页浏览食谱所述安装硒。我在这个食谱中使用了 PhantomJS，但这并不是硬性要求。您可以使用 Selenium 支持的任何其他浏览器。我的修改在github.com/ivanidris/5…的 0.0.1 标签下(2015 年 9 月检索)。下载其中一个源档案并将其解压缩。导航至crawler目录及其code子目录。

使用以下命令启动(可选步骤)爬虫(我以 CNN 为例):

$ python crawl.py edition.cnn.com

怎么做…

您可以在本书的代码包中使用带有链接的 CSV 文件，也可以按照我在上一节中解释的那样自己制作。以下过程描述了如何创建新闻文章的文本语料库(参见本书代码包中的download_html.py文件):

进口情况如下:

import dautil as dl
import csv
import os
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 urllib.parse as urlparse
import urllib.request as urlrequest

定义以下全局常数:

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

定义以下函数从一个 HTML 页面中提取文本并保存:

def write_text(fname):
    elems = []

    try:
        DRIVER.get(dl.web.path2url(fname))

        elems = WebDriverWait(DRIVER, NAP_SECONDS).until(
            EC.presence_of_all_elements_located((By.XPATH, '//p'))
        )

        LOGGER.info('Elems', elems)

        with open(fname.replace('.html', '_phantomjs.html'), 'w') as pjs_file:
            LOGGER.warning('Writing to %s', pjs_file.name)
            pjs_file.write(DRIVER.page_source)

    except Exception:
        LOGGER.error("Error processing HTML", exc_info=True)

    new_name = fname.replace('html', 'txt')

    if not os.path.exists(new_name):
        with open(new_name, 'w') as txt_file:
            LOGGER.warning('Writing to %s', txt_file.name)

            lines = [e.text for e in elems]
            LOGGER.info('lines', lines)
            txt_file.write(' \n'.join(lines))

定义如下main()函数，读取带有链接的 CSV 文件，调用前面步骤中的函数:

def main():
    filedir = os.path.join(dl.data.get_data_dir(), 'edition.cnn.com')

    with open('saved_urls.csv') as csvfile:
        reader = csv.reader(csvfile)

        for line in reader:
            timestamp, count, basename, url = line
            fname = '_'.join([count, basename])
            fname = os.path.join(filedir, fname)

            if not os.path.exists(fname):
                dl.data.download(url, fname)

            write_text(fname)

if __name__ == '__main__':
    DRIVER.implicitly_wait(NAP_SECONDS)
    main()
    DRIVER.quit()

处理非 ASCII 文本和 HTML 实体

HTML 不像来自数据库查询或 Pandas 的数据那样结构化。你可能会被用正则表达式或字符串函数操纵 HTML。然而，这种方法只在有限的情况下有效。您最好使用专门的 Python 库来处理 HTML。在这个食谱中，我们将使用 lxml 库的clean_html()功能。这个函数从一个网页中剥离所有的 JavaScript 和 CSS。

美国信息交换标准码 ( ASCII )在 2007 年底之前一直是互联网上的主导编码标准，UTF-8 (8 位 Unicode)占据首位。ASCII 仅限于英文字母，不支持不同语言的字母。Unicode 对字母有更广泛的支持。但是，我们有时需要将自己限制在 ASCII 中，所以这个食谱为您提供了一个如何忽略非 ASCII 字符的示例。

做好准备

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

$ pip install lxml
$ conda install lxml

我用 Anaconda 的 lxml 3.4.2 测试了代码。

怎么做…

代码在本书的代码包中的processing_html.py文件中，分为以下步骤:

进口情况如下:

from lxml.html.clean import clean_html
from difflib import Differ
import unicodedata
import dautil as dl

PRINT = dl.log_api.Printer()

定义以下函数来区分两个文件:

def diff_files(text, cleaned):
    d = Differ()
    diff = list(d.compare(text.splitlines(keepends=True),
                          cleaned.splitlines(keepends=True)))
    PRINT.print(diff)

以下代码块打开一个 HTML 文件，对其进行清理，并将清理后的文件与原始文件进行比较:

with open('460_cc_phantomjs.html') as html_file:
    text = html_file.read()
    cleaned = clean_html(text)
    diff_files(text, cleaned)
    PRINT.print(dl.web.find_hrefs(cleaned))

下面的代码片段演示了对非 ASCII 文本的处理:

bulgarian = 'Питон is Bulgarian for Python'
PRINT.print('Bulgarian', bulgarian)
PRINT.print('Bulgarian ignored', unicodedata.normalize('NFKD', bulgarian).encode('ascii', 'ignore'))

最终结果请参考下面的截图(为了简洁起见，我省略了一些输出):

How to do it…

另见

lxml 文档是位于lxml.de/index.html的(2015 年 9 月检索)

实现关联表

关联表充当数据库表之间的桥梁，数据库表具有多对多的关系。该表包含与它所连接的表的主键相链接的外键。

在这个食谱中，我们将把网页和网页内的链接联系起来。一个页面有很多链接，链接可以在很多页面中。我们将只关注其他网站的链接，但这不是必需的。如果您试图在本地机器上复制一个网站进行测试或分析，那么您还需要存储图像和 JavaScript 链接。看看下面的关系模式图:

Implementing association tables

做好准备

我用 Anaconda 安装了 SQLAlchemy 0.9.9，如下所示:

$ conda install sqlalchemy

如果愿意，也可以使用以下命令安装 SQLAlchemy:

$ pip install sqlalchemy

怎么做…

本书代码包中impl_association.py文件的以下代码实现了关联表模式:

进口情况如下:

from sqlalchemy import create_engine
from sqlalchemy import Column
from sqlalchemy import ForeignKey
from sqlalchemy import Integer
from sqlalchemy import String
from sqlalchemy import Table
from sqlalchemy.orm import backref
from sqlalchemy.orm import relationship
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy.orm import sessionmaker
from sqlalchemy.exc import IntegrityError
import dautil as dl
import os

Base = declarative_base()

定义以下类来表示网页:

class Page(Base):
    __tablename__ = 'pages'
    id = Column(Integer, primary_key=True)
    filename = Column(String, nullable=False, unique=True)
    links = relationship('Link', secondary='page_links')

    def __repr__(self):
        return "Id=%d filename=%s" %(self.id, self.filename)

定义以下类来表示网络链接:

class Link(Base):
    __tablename__ = 'links'
    id = Column(Integer, primary_key=True)
    url = Column(String, nullable=False, unique=True)

    def __repr__(self):
        return "Id=%d url=%s" %(self.id, self.url)

定义以下类来表示页面和链接之间的关联:

class PageLink(Base):
    __tablename__ = 'page_links'
    page_id = Column(Integer, ForeignKey('pages.id'), primary_key=True)
    link_id = Column(Integer, ForeignKey('links.id'), primary_key=True)
    page = relationship('Page', backref=backref('link_assoc'))
    link = relationship('Link', backref=backref('page_assoc'))

    def __repr__(self):
        return "page_id=%s link_id=%s" %(self.page_id, self.link_id)

定义以下功能来浏览 HTML 文件并更新相关表格:

def process_file(fname, session):
    with open(fname) as html_file:
        text = html_file.read()

        if dl.db.count_where(session, Page.filename, fname):
            # Cowardly refusing to continue
            return

        page = Page(filename=fname)
        hrefs = dl.web.find_hrefs(text)

        for href in set(hrefs):
            # Only saving http links
            if href.startswith('http'):
                if dl.db.count_where(session, Link.url, href):
                    continue

                link = Link(url=href)
                session.add(PageLink(page=page, link=link))

        session.commit()

定义以下函数来填充数据库:

def populate():
    dir = dl.data.get_data_dir()
    path = os.path.join(dir, 'crawled_pages.db')
    engine = create_engine('sqlite:///' + path)
    DBSession = sessionmaker(bind=engine)
    Base.metadata.create_all(engine)
    session = DBSession()

    files  = ['460_cc_phantomjs.html', '468_live_phantomjs.html']

    for file in files:
        process_file(file, session)

    return session

下面的代码片段使用了我们定义的函数和类:

if __name__ == "__main__":
    session = populate()
    printer = dl.log_api.Printer(nelems=3)
    pages = session.query(Page).all()
    printer.print('Pages', pages)

    links = session.query(Link).all()
    printer.print('Links', links)

    page_links = session.query(PageLink).all()
    printer.print('PageLinks', page_links)

最终结果参见下面的截图:

How to do it…

设置数据库迁移脚本

你在编程课上学到的第一件事是，没有人能第一时间把复杂的程序做对。软件随着时间的推移而发展，我们希望最好的。自动化测试中的自动化有助于确保我们的程序随着时间的推移而改进。然而，当谈到进化数据库模式时，自动化似乎并不那么明显。尤其是在大型企业中，数据库模式是数据库管理员和专家的领域。当然，还有与更改模式相关的安全和操作问题，在生产数据库中更是如此。在任何情况下，您都可以在本地测试环境中实现数据库迁移，并为生产团队记录建议的更改。

我们将使用 Alembic 来演示如何着手设置迁移脚本。在我看来，Alembic 是这项工作的合适工具，尽管它截至 2015 年 9 月仍处于测试阶段。

做好准备

使用以下命令安装 Alembic:

$ pip install alembic

Alembic 依赖于 SQLAlchemy，如果需要，它会自动安装 SQLAlchemy。你现在应该有alembic命令在你的路径上。这一章我用了 Alembic 0.8.2。

怎么做…

以下步骤描述了如何设置一个迁移步骤。当我们在目录中运行 Alembic 初始化脚本时，它会创建一个alembic目录和一个名为alembic.ini的配置文件:

导航到适当的目录并初始化迁移项目，如下所示:
```
$ alembic init alembic
```
根据需要编辑alembic.ini文件。例如，更改sqlalchemy.url属性以指向正确的数据库。

另见

相关的 Alembic 文档位于alembic.readthedocs.org/en/latest/t…(2015 年 9 月检索)

向现有表格添加表格列

如果我们使用对象关系映射器(ORM)，比如 SQLAlchemy，我们将类映射到表，将类属性映射到表列。通常，由于新的业务需求，我们需要添加一个表列和相应的类属性。我们可能需要在添加后立即填充该列。

如果我们处理生产数据库，那么您可能没有直接访问权限。幸运的是，我们可以用 Alembic 生成 SQL，数据库管理员可以查看它。

做好准备

参考设置数据库迁移脚本配方。

怎么做…

Alembic 有自己的版本控制系统，需要额外的表。它还用生成的 Python 代码文件在alembic目录下创建一个versions目录。我们需要在这些文件中指定迁移所需的更改类型:

创建新版本，如下所示:
```
$ alembic revision -m "Add a column"
```

打开生成的 Python 文件(例如27218d73000_add_a_column.py)。用下面的代码替换其中的两个函数，增加link_type字符串列:

def upgrade():
    # MODIFIED Ivan Idris
    op.add_column('links', sa.Column('link_type', sa.String(20)))

def downgrade():
    # MODIFIED Ivan Idris
    op.drop_column('links', 'link_type')

生成 SQL，如下所示:
```
$ alembic upgrade head --sql
```

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

How to do it…

创建表后添加索引

指数是计算中的一个通用概念。这本书还有一个索引，可以更快地查找，将概念与页码匹配起来。索引占用空间；就这本书而言，几页。数据库索引还有一个额外的缺点，由于更新索引的额外开销，它们会使插入和更新变得更慢。通常，主键和外键会自动获得索引，但这取决于数据库实现。

添加索引不能掉以轻心，最好在咨询数据库管理员后进行。Alembic 的索引添加功能类似于我们在中看到的功能，将一个表列添加到现有的表配方中。

做好准备

参考设置数据库迁移脚本配方。

怎么做…

这个配方和在现有的表格配方上增加一个表格列有一些重叠，所以我就不重复所有的细节了:

创建新版本，如下所示:
```
$ alembic revision -m "Add indices"
```

打开生成的 Python 文件(例如21579ecccd8_add_indices.py)并修改代码，使其具有以下功能，这些功能负责添加索引:

def upgrade():
    # MODIFIED Ivan Idris
    op.create_index('idx_links_url', 'links', ['url'])
    op.create_index('idx_pages_filename', 'pages', ['filename'])

def downgrade():
    # MODIFIED Ivan Idris
    op.drop_index('idx_links_url')
    op.drop_index('idx_pages_filename')

生成 SQL，如下所示:
```
$ alembic upgrade head --sql
```

最终结果参见下面的截图:

How to do it…

它是如何工作的…

create_index()函数添加给定索引名、表和表列列表的索引。drop_index()函数则相反，移除给定索引名的索引。

另见

关于数据库索引的维基百科页面位于en.wikipedia.org/wiki/Databa…(检索于 2015 年 9 月)

设置测试网络服务器

在第 1 章、为可再现数据分析奠定基础中，我们讨论了为什么单元测试是一个好主意。纯粹主义者会告诉你，你只需要单元测试。然而，普遍的共识是，更高级别的测试也是有用的。

显然，这本书是关于数据分析的，而不是关于 web 开发的。不过，通过网站或网络服务分享你的结果或数据是一个常见的要求。当你挖掘网络或者做其他与网络相关的事情时，经常需要重现某些用例，比如登录表单。正如您对成熟语言的期望，Python 有许多优秀的网络框架。我选择了 Flask，这是一个简单的 Pythonic 网络框架，因为它看起来很容易设置，但是你应该使用你自己的判断，因为我不知道你的要求是什么。

做好准备

我用 Anaconda 的 Flask 0.10.1 测试了代码。用conda或pip安装烧瓶，如下所示:

$ conda install flask
$ pip install flask

怎么做…

在这个食谱中，我们将设置一个带有登录表单的安全页面，您可以使用它进行测试。代码由app.py Python 文件和templates目录下的 HTML 文件组成(我就不详细讨论 HTML 了):

进口情况如下:

from flask import Flask
from flask import render_template
from flask import request
from flask import redirect
from flask import url_for

app = Flask(__name__)

定义以下功能来处理主页请求:

@app.route('/')
def home():
    return "Test Site"

定义以下功能来处理登录尝试:

@app.route('/secure', methods=['GET', 'POST'])
def login():
    error = None
    if request.method == 'POST':
        if request.form['username'] != 'admin' or\
                request.form['password'] != 'admin':
            error = 'Invalid password or user name.'
        else:
            return redirect(url_for('home'))
    return render_template('admin.html', error=error)

以下块运行服务器(不要将debug=True用于面向公众的网站):
```
if __name__ == '__main__':
    app.run(debug=True)
```
运行$ python app.py，在http://127.0.0.1:5000/和http://127.0.0.1:5000/secure打开网页浏览器。

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

How to do it…

用事实表和维度表实现星型模式

星型模式是一种便于报告的数据库模式。星型模式适用于事件的处理，如网站访问、广告点击或金融交易。事件信息(指标，如温度或购买金额)存储在事实表中，事实表链接到小得多的维度表。星型模式是非规范化的，这将完整性检查的责任放在了应用程序代码上。因此，我们应该只以受控的方式写入数据库。如果使用 SQLAlchemy 进行大容量插入，应该选择核心应用编程接口而不是 ORM 应用编程接口，或者使用直接的 SQL。你可以在docs.sqlalchemy.org/en/rel_1_0/…阅读更多原因(2015 年 9 月检索)。

时间是报告中的一个常见维度。例如，我们可以在维度表中存储每日天气测量的日期。对于数据中的每个日期，我们可以保存日期、年份、月份和一年中的某一天。我们可以在处理事件之前预先填充这个表，然后根据需要添加新的日期。如果我们假设只需要维护数据库一个世纪，我们甚至不需要向时间维度表中添加新记录。在这种情况下，我们将在时间维度表中预先填充我们想要支持的范围内的所有可能日期。如果我们正在处理二进制或分类变量，预先填充维度表也是可能的。

在本食谱中，我们将为中描述的直接营销数据实施星型模式。数据在一个直接营销活动的 CSV 文件中。为了简单起见，我们将忽略一些列。作为一个指标，我们将带采购金额的spend列。对于维度，我选择了渠道(电话、网络或多渠道)、邮政编码(农村、郊区或城市)和细分市场(男性、女性或没有电子邮件)。请参考以下实体关系图:

Implementing a star schema with fact and dimension tables

怎么做…

以下代码下载数据，加载到数据库中，然后查询数据库(参考本书代码包中的star_schema.py文件):

进口情况如下:

from sqlalchemy import Column
from sqlalchemy import create_engine
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy import ForeignKey
from sqlalchemy import Integer
from sqlalchemy import String
from sqlalchemy.orm import sessionmaker
from sqlalchemy import func
import dautil as dl
from tabulate import tabulate
import sqlite3
import os
from joblib import Memory

Base = declarative_base()
memory = Memory(cachedir='.')

定义以下类来表示邮政编码维度:

class DimZipCode(Base):
    __tablename__ = 'dim_zip_code'
    id = Column(Integer, primary_key=True)
    # Urban, Suburban, or Rural.
    zip_code = Column(String(8), nullable=False, unique=True)

定义以下类来表示线段尺寸:

class DimSegment(Base):
    __tablename__ = 'dim_segment'
    id = Column(Integer, primary_key=True)
    # Mens E-Mail, Womens E-Mail or No E-Mail
    segment = Column(String(14), nullable=False, unique=True)

定义以下类来表示通道维度:

class DimChannel(Base):
    __tablename__ = 'dim_channel'
    id = Column(Integer, primary_key=True)
    channel = Column(String)

定义以下类来表示事实表:

class FactSales(Base):
    __tablename__ = 'fact_sales'
    id = Column(Integer, primary_key=True)
    zip_code_id = Column(Integer, ForeignKey('dim_zip_code.id'),
                         primary_key=True)
    segment_id = Column(Integer, ForeignKey('dim_segment.id'),
                        primary_key=True)
    channel_id = Column(Integer, ForeignKey('dim_channel.id'),
                        primary_key=True)

    # Storing amount as cents
    spend = Column(Integer)

    def __repr__(self):
        return "zip_code_id={0} channel_id={1} segment_id={2}".format(
            self.zip_code_id, self.channel_id, self.segment_id)

定义以下功能创建 SQLAlchemy 会话:

def create_session(dbname):
    engine = create_engine('sqlite:///{}'.format(dbname))
    DBSession = sessionmaker(bind=engine)
    Base.metadata.create_all(engine)

    return DBSession()

定义以下功能来填充段维度表:

def populate_dim_segment(session):
    options = ['Mens E-Mail', 'Womens E-Mail', 'No E-Mail']

    for option in options:
        if not dl.db.count_where(session, DimSegment.segment, option):
            session.add(DimSegment(segment=option))

    session.commit()

定义以下函数来填充邮政编码维度表:

def populate_dim_zip_code(session):
    # Note the interesting spelling
    options = ['Urban', 'Surburban', 'Rural']

    for option in options:
        if not dl.db.count_where(session, DimZipCode.zip_code, option):
            session.add(DimZipCode(zip_code=option))

    session.commit()

定义以下函数来填充通道维度表:

def populate_dim_channels(session):
    options = ['Phone', 'Web', 'Multichannel']

    for option in options:
        if not dl.db.count_where(session, DimChannel.channel, option):
            session.add(DimChannel(channel=option))

    session.commit()

定义以下函数来填充事实表(出于性能原因，它使用直接的 SQL:

```py

def load(csv_rows, session, dbname):
    channels = dl.db.map_to_id(session, DimChannel.channel)
    segments = dl.db.map_to_id(session, DimSegment.segment)
    zip_codes = dl.db.map_to_id(session, DimZipCode.zip_code)
    conn = sqlite3.connect(dbname)
    c = conn.cursor()
    logger = dl.log_api.conf_logger(__name__)

    for i, row in enumerate(csv_rows):
        channel_id = channels[row['channel']]
        segment_id = segments[row['segment']]
        zip_code_id = zip_codes[row['zip_code']]
        spend = dl.data.centify(row['spend'])

        insert = "INSERT INTO fact_sales (id, segment_id,\
            zip_code_id, channel_id, spend) VALUES({id}, \
            {sid}, {zid}, {cid}, {spend})"
        c.execute(insert.format(id=i, sid=segment_id,
                                zid=zip_code_id, cid=channel_id,     spend=spend))

        if i % 1000 == 0:
            logger.info("Progress %s/64000", i)
            conn.commit()

    conn.commit()
    c.close()
    conn.close()
```

11. 定义以下函数下载并解析数据:

```py
@memory.cache
def get_and_parse():
    out = dl.data.get_direct_marketing_csv()
    return dl.data.read_csv(out)
```

12. 以下模块使用我们定义的函数和类【T3:

```py
if __name__ == "__main__":
    dbname = os.path.join(dl.data.get_data_dir(), 'marketing.db')
    session = create_session(dbname)
    populate_dim_segment(session)
    populate_dim_zip_code(session)
    populate_dim_channels(session)

    if session.query(FactSales).count() < 64000:
        load(get_and_parse(), session, dbname)

    fsum = func.sum(FactSales.spend)
    query = session.query(DimSegment.segment, DimChannel.channel,
                          DimZipCode.zip_code, fsum)
    dim_cols = (DimSegment.segment, DimChannel.channel, DimZipCode.zip_code)
    dim_entities = [dl.db.entity_from_column(col) for col in dim_cols]
    spend_totals = query.join(FactSales,
                              *dim_entities)\
                        .group_by(*dim_cols).order_by(fsum.desc()).all()
    print(tabulate(spend_totals, tablefmt='psql',
                   headers=['Segment', 'Channel', 'Zip Code', 'Spend']))
```

有关最终结果(以美分为单位的支出金额)，请参考以下屏幕截图:

How to do it…

另见

en.wikipedia.org/wiki/Star_s… 的星模式维基百科页面(2015 年 9 月检索)

使用 HDFS

Hadoop 分布式文件系统 ( HDFS )是大数据 Hadoop 框架的存储组件。HDFS 是一个分布式文件系统，它在多个系统上传播数据，其灵感来自于谷歌用于其搜索引擎的谷歌文件系统。HDFS 需要一个 Java 运行时环境 ( JRE )，它使用一个NameNode服务器来跟踪文件。系统还会复制数据，这样丢失几个节点不会导致数据丢失。HDFS 的典型用例是处理大型只读文件。Apache Spark 也包含在本章中，它也可以使用 HDFS。

做好准备

安装 Hadoop 和一个 JRE。因为这些不是 Python 框架，所以您必须检查什么程序适合您的操作系统。这个食谱我用的是 Hadoop 2.7.1 和 Java 1.7.0_60。这可能是一个复杂的过程，但网上有许多资源可以帮助您解决特定系统的故障。

怎么做…

我们可以用在你的 Hadoop 安装中找到的几个 XML 文件来配置 HDFS。本节中的一些步骤仅作为示例，您应该根据您的操作系统、环境和个人偏好来实施它们:

编辑core-site.xml文件，使其具有以下内容(注释省略):

<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet type="text/xsl" href="configuration.xsl"?>
<configuration>
    <property>
        <name>fs.default.name</name>
        <value>hdfs://localhost:8020</value>
    </property>
</configuration>

编辑hdfs-site.xml文件，使其具有以下内容(注释省略)，将每个文件的复制设置为仅 1，以在本地运行 HDFS:

<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet type="text/xsl" href="configuration.xsl"?>
<configuration>
    <property>
        <name>dfs.replication</name>
        <value>1</value>
    </property>
</configuration>

如有必要，在您的系统上启用远程登录来 SSH 到本地主机并生成密钥(Windows 用户可以使用 putty):
```
$ ssh-keygen -t dsa -f ~/.ssh/id_dsa
$ cat ~/.ssh/id_dsa.pub >> ~/.ssh/authorized_keys
```
从 Hadoop 目录的根目录格式化文件系统:
```
$ bin/hdfs namenode –format
```
启动 NameNode 服务器，如下(相反命令为$ sbin/stop-dfs.sh ):
```
$ sbin/start-dfs.sh
```
使用以下命令在 HDFS 创建目录:
```
$ hadoop fs -mkdir direct_marketing
```
或者，如果您想使用火花配方中的direct_marketing.csv文件，您需要将其复制到 HDFS，如下所示:
```
$ hadoop fs -copyFromLocal <path to file>/direct_marketing.csv direct_marketing
```

另见

位于的 HDFS 用户指南

设置火花

Apache Spark 是 Hadoop 生态系统中的一个项目(参考使用 HDFS 配方)，据称其性能优于 Hadoop 的 MapReduce。Spark 尽可能将数据加载到内存中，对机器学习有很好的支持。在用 Spark 配方聚类数据中，我们将通过 Spark 应用机器学习算法。

Spark 可以独立工作，但它被设计为使用 HDFS 与 Hadoop 一起工作。弹性分布式数据集 ( RDDs )是 Spark 的中心结构，代表分布式数据。Spark 对 Scala(一种 JVM 语言)有很好的支持，对 Python 的支持有些滞后。例如，pyspark API 中对流的支持有点滞后。Spark 也有 DataFrames 的概念，但是它不是通过 pandas 实现的，而是通过 Spark 实现的。

做好准备

从spark.apache.org/downloads.h…的下载页面下载星火(2015 年 9 月检索)。我下载了星火 1.5.0 的spark-1.5.0-bin-hadoop2.6.tgz档案。

将归档文件解压缩到适当的目录中。

怎么做…

以下步骤通过几个可选步骤说明了 Spark 的基本设置:

如果要使用不同于系统 Python 的 Python 版本，请通过操作系统的 GUI 或 CLI 设置PYSPARK_PYTHON环境变量，如下所示:
```
$ export PYSPARK_PYTHON=/path/to/anaconda/bin/python
```

设置SPARK_HOME环境变量，如下所示:

$ export SPARK_HOME=<path/to/spark/>spark-1.5.0-bin-hadoop2.6

将python目录添加到您的PYTHONPATH环境变量中，如下所示:
```
$ export PYTHONPATH=$SPARK_HOME/python:$PYTHONPATH
```

将py4j的 ZIP 添加到您的PYTHONPATH环境变量中，如下所示:

$ export PYTHONPATH=$SPARK_HOME/python/lib/py4j-0.8.2.1-src.zip:$PYTHONPATH

如果 Spark 的日志记录过于冗长，将$SPARK_HOME/conf目录中的log4j.properties.template文件复制到log4j.properties并将 INFO 级别更改为 WARN。

另见

官方星火网站在spark.apache.org/(2015 年 9 月检索)

用 Spark 聚类数据

在上一个食谱设置火花中，我们介绍了火花的基本设置。如果你按照使用 HDFS 食谱，你可以选择提供来自 Hadoop 的数据。在这种情况下，您需要以这种方式指定文件的网址，hdfs://hdfs-host:port/path/direct_marketing.csv。

我们将使用与我们在中使用的相同的数据，用事实和维度表实现星型模式。但是，这次我们将使用spend、history和recency列。第一列对应直销活动后的最近购买金额，第二列对应历史购买金额，第三列对应最近几个月的购买金额。数据描述见http://blog . minethatdata . com/2008/03/minethatdata-e-mail-analytics-and-data . html(2015 年 9 月检索)。我们将应用流行的 K-means 机器学习算法对数据进行聚类。第九章、集成学习和降维，更关注机器学习算法。K-means 算法试图为给定数量的聚类的数据集找到最佳聚类。我们应该要么知道这个数字，要么通过反复试验找到它。在本食谱中，我通过集内平方和误差 ( WSSSE )也称为集内平方和 ( WCSS )来评估聚类。该度量计算每个点和其指定簇之间距离的平方误差之和。您可以在第 10 章、评估分类器、回归器和聚类中了解更多关于评估指标的信息。

做好准备

按照设置火花配方中的说明进行操作。

怎么做…

该配方的代码在本书代码包的clustering_spark.py文件中:

进口情况如下:

from pyspark.mllib.clustering import KMeans
from pyspark import SparkContext
import dautil as dl
import csv
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.colors import Normalize

定义以下函数计算误差:

def error(point, clusters):
    center = clusters.centers[clusters.predict(point)]

    return dl.stats.wssse(point, center)

读取并解析数据，如下所示:

sc = SparkContext()
csv_file = dl.data.get_direct_marketing_csv()
lines = sc.textFile(csv_file)
header = lines.first().split(',')
cols_set = set(['recency', 'history', 'spend'])
select_cols = [i for i, col in enumerate(header) if col in cols_set]

设置以下 rdd:

header_rdd = lines.filter(lambda l: 'recency' in l)
noheader_rdd = lines.subtract(header_rdd)
temp = noheader_rdd.map(lambda v: list(csv.reader([v]))[0])\
                   .map(lambda p: (int(p[select_cols[0]]),
                        dl.data.centify(p[select_cols[1]]),
                        dl.data.centify(p[select_cols[2]])))

# spend > 0
temp = temp.filter(lambda x: x[2] > 0)

用 k-means 算法对数据进行聚类:

points = []
clusters = None

for i in range(2, 28):
    clusters = KMeans.train(temp, i, maxIterations=10,
                            runs=10, initializationMode="random")

    val = temp.map(lambda point: error(point, clusters))\
              .reduce(lambda x, y: x + y)
    points.append((i, val))

绘制集群，如下所示:

dl.options.mimic_seaborn()
fig, [ax, ax2] = plt.subplots(2, 1)
ax.set_title('k-means Clusters')
ax.set_xlabel('Number of clusters')
ax.set_ylabel('WSSSE')
dl.plotting.plot_points(ax, points)

collected = temp.collect()
recency, history, spend = zip(*collected)
indices = [clusters.predict(c) for c in collected]
ax2.set_title('Clusters for spend, history and recency')
ax2.set_xlabel('history (cents)')
ax2.set_ylabel('spend (cents)')
markers = dl.plotting.map_markers(indices)
colors = dl.plotting.sample_hex_cmap(name='hot', ncolors=len(set(recency)))

for h, s, r, m in zip(history, spend, recency, markers):
    ax2.scatter(h, s, s=20 + r, marker=m, c=colors[r-1])

cma = mpl.colors.ListedColormap(colors, name='from_list', N=None)
norm = Normalize(min(recency), max(recency))
msm = mpl.cm.ScalarMappable(cmap=cma, norm=norm)
msm.set_array([])
fig.colorbar(msm, label='Recency')

for i, center in enumerate(clusters.clusterCenters):
    recency, history, spend = center
    ax2.text(history, spend, str(i))

plt.tight_layout()
plt.show()

最终结果见以下截图(图中数字对应聚类中心):

How to do it…

它是如何工作的…

k 均值聚类将数据点分配给 k 个聚类。聚类问题不是直接可以解决的，但是我们可以应用启发式算法，获得一个可接受的结果。k 均值算法在两个步骤之间迭代，不包括 k 均值的初始化(通常是随机的):

给每个数据点分配一个 WCSS 均值最低的聚类
重新计算聚类中心作为聚类点坐标的平均值

当集群分配变得稳定时，算法停止。

还有更多…

Spark 1.5.0 增加了对流 K 均值的实验支持。由于这些新特性的实验性质，我决定不详细讨论它们。我在本书的代码包streaming_clustering.py文件中添加了以下示例代码:

import dautil as dl
from pyspark.mllib.clustering import StreamingKMeansModel
from pyspark import SparkContext

csv_file = dl.data.get_direct_marketing_csv()
csv_rows = dl.data.read_csv(csv_file)

stkm = StreamingKMeansModel(28 * [[0., 0., 0.]], 28 * [1.])
sc = SparkContext()

for row in csv_rows:
    spend = dl.data.centify(row['spend'])

    if spend > 0:
        history = dl.data.centify(row['history'])
        data = sc.parallelize([[int(row['recency']),
                               history, spend]])
        stkm = stkm.update(data, 0., 'points')

print(stkm.centers)

另见

用 Python 构建机器学习系统、威利·里歇特和路易斯·佩德罗科埃略(2013)
维基百科关于 K-means 聚类的页面在 en.wikipedia.org/wiki/K-mean… 年 9 月检索)

六、信号处理和时间序列

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

用周期图进行光谱分析
用韦尔奇方法估计功率谱密度
分析峰值
测量相位同步
指数平滑法
评估平滑度
使用隆布-斯卡格尔周期图
分析音频的频谱
用离散余弦变换分析信号
块自举时间序列数据
移动块自举时间序列数据
应用离散小波变换

简介

时间是科学和日常生活中的一个重要维度。时间序列数据非常丰富，需要特殊的技术。通常，我们对趋势和季节性或周期性感兴趣。在数学术语中，这意味着我们试图用(通常是线性的)多项式或三角函数，或两者的组合来表示数据。

当我们研究季节性时，我们通常区分时域和频域分析。在时域，我们可以使用十几个 Pandas 函数来滚动窗口。我们还可以平滑数据以去除噪声，同时希望保持足够的信号。平滑在许多方面类似于拟合，这很方便，因为我们可以重用一些我们知道的回归工具。

为了进入频域，我们应用变换，例如快速傅立叶变换和离散余弦变换。然后我们可以用周期图进一步分析信号。

用周期图进行光谱分析

我们可以认为周期性的信号是由多个频率组成的。比如声音是由多种音调组成，光是由多种颜色组成的。频率范围称为频谱。当我们分析一个信号的频谱时，我们很自然地会看到这个信号的傅里叶变换的结果。周期图对此进行了扩展，并等于傅里叶变换的平方幅度，如下所示:

Spectral analysis with periodograms

我们将查看以下变量的周期图:

KNMI De Bilt 气象数据中的雨量值
雨量值的二阶差(相当于微积分中的二阶导数)
使用 365 天窗口的雨数值的滚动总和
使用 365 天窗口的雨值的滚动平均值

怎么做...

进口情况如下:

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

加载数据如下:

fs = 365
rain = dl.data.Weather.load()['RAIN'].dropna()

定义以下函数绘制周期图:

def plot_periodogram(arr, ax):
    f, Pxx_den = signal.periodogram(arr, fs)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('PSD')
    ax.semilogy(f, Pxx_den)

用以下代码绘制周期图:

sp = dl.plotting.Subplotter(2, 2, context)
sp.label()
plot_periodogram(rain, sp.ax)
sp.label(advance=True)
plot_periodogram(np.diff(rain, 2), sp.ax)
sp.label(advance=True)
plot_periodogram(pd.rolling_sum(rain, fs).dropna(), sp.ax)
sp.label(advance=True)
plot_periodogram(pd.rolling_mean(rain, fs).dropna(), sp.ax)
HTML(sp.exit())

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

How to do it...

代码来自本书代码包演示周期图中的periodograms.ipynb文件。

另见

https://docs . scipy . org/doc/scipy-0 . 16 . 0/reference/generated/scipy . signal . periogram . html # scipy . signal . periogram上的功能文档(2015 年 9 月检索)
en.wikipedia.org/wiki/Spectr… 相关维基百科页面(2015 年 9 月检索)

用韦尔奇方法估计功率谱密度

韦尔奇方法是对周期图技术的改进(它降低了噪声)，以 P.D .韦尔奇命名。通过以下步骤降低功率谱的噪声:

我们用固定数量的重叠点分割信号。如果重叠为 0，那么我们有巴特利特方法。
在时域中，我们将窗口函数应用于步骤 1 的每个片段。
我们计算每个片段的周期图，如光谱分析周期图配方中所述。
我们对周期图进行平均，从而减少噪音。平均有效地平滑了信号。然而，我们现在处理的是频率仓(就像直方图)。

我们还将探讨 法诺因子，给出如下:

Estimating power spectral density with the Welch method

这是一个加窗的方差均值比。除以平均值基本上将这些值归一化，我们得到一个归一化的离差度量。作为输入数据，我们将使用温度数据。

怎么做...

进口情况如下:

from scipy import signal
import matplotlib.pyplot as plt
import dautil as dl
from IPython.display import HTML

加载数据并计算 Fano 因子:

fs = 365
temp = dl.data.Weather.load()['TEMP'].dropna()
fano_factor = dl.ts.fano_factor(temp, fs)

定义以下函数绘制周期图:

def plot_welch(arr, ax):
    f, Pxx_den = signal.welch(arr, fs)
    ax.semilogy(f, Pxx_den)

绘制输入数据和对应的周期图:

sp = dl.plotting.Subplotter(2, 2, context)
temp.plot(ax=sp.ax)
sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))
sp.label(advance=True)
sp.ax.plot(temp.index, fano_factor)
sp.label(advance=True)
plot_welch(temp, sp.ax)
sp.label(advance=True)
plot_welch(fano_factor.dropna(), sp.ax)
HTML(sp.exit())

最终结果参见以下截图:

How to do it...

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

另见

维基百科页面关于en.wikipedia.org/wiki/Fano_f…的法诺因子(2015 年 9 月检索)
关于 en.wikipedia.org/wiki/Welch'… T2 韦尔奇方法的维基百科页面(2015 年 9 月检索)
welch()功能记录在https://docs . scipy . org/doc/scipy/reference/generated/scipy . signal . Welch . html # scipy . signal . Welch(2015 年 9 月检索)

分析峰值

峰的分析与谷的分析相似，因为两者都是极值。SciPy 具有找到相对最大值的argrelmax()函数。当我们将此函数应用于每日温度值时，它不仅会发现夏季的热天，还会发现冬季的热天，除非我们使该函数考虑更大的时间范围。当然，我们也可以检查数值是否高于阈值，或者仅使用先验知识选择夏季数据。

当我们分析时间序列数据中的峰值时，我们可以应用两种方法。第一种方法是考虑一年、一个月或另一个固定时间间隔内的最高峰值，并用这些值构建一个序列。第二种方法是将任何高于阈值的值定义为峰值。在这个食谱中，我们将使用第 95 个百分位数作为阈值。在这种方法中，我们可以在一个序列中有多个峰值。例如，在热浪的情况下，长条纹会产生负面影响。

怎么做...

进口情况如下:

import dautil as dl
from scipy import signal
import matplotlib.pyplot as plt
import seaborn as sns
from IPython.display import HTML

加载并重新采样数据:

temp = dl.data.Weather.load()['TEMP'].dropna()
monthly = temp.resample('M')

绘制峰值，注意还考虑了冬季的热天:

sp = dl.plotting.Subplotter(2, 2, context)
max_locs = signal.argrelmax(monthly.values)
sp.ax.plot(monthly.index, monthly, label='Monthly means')
sp.ax.plot(monthly.index[max_locs], monthly.values[max_locs], 
           'o', label='Tops')
sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

绘制年度最大系列:

annual_max = dl.ts.groupby_year(temp).max()
sp.next_ax().plot(annual_max.index, annual_max, label='Annual Maximum Series')
dl.plotting.plot_polyfit(sp.ax, annual_max.index, annual_max.values)
sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

绘制超过第 95 百分位阈值的全年最长高温天数:

_, threshhold = dl.stats.outliers(temp, method='percentiles')
over_threshhold = temp > threshhold
streaks = dl.ts.groupby_year(over_threshhold).apply(
    lambda x: dl.collect.longest_streak(x, 1))
sp.next_ax().plot(streaks.index, streaks)
dl.plotting.plot_polyfit(sp.ax, streaks.index, streaks.values)
over_threshhold = dl.ts.groupby_year(over_threshhold).mean()
sp.label()

绘制年度最大系列分布:

sp.label(advance=True)
sns.distplot(annual_max, ax=sp.ax)
sp.label(xlabel_params=dl.data.Weather.get_header('TEMP'))
HTML(sp.exit())

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

How to do it...

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

另见

argrelmax()功能记录在https://docs . scipy . org/doc/scipy/reference/generated/scipy . signal . argrelmax . html(2015 年 9 月检索)

测量相位同步

两个信号可以完全同步，也可以不同步，或者介于两者之间。我们通常以弧度测量相位同步。 瞬时相位的相关量可以用 NumPy angle()功能测量。对于实值数据，我们需要获得信号的解析表示，它由希尔伯特变换给出。希尔伯特变换也有 SciPy 和 NumPy 版本。

互相关使用滑动内积测量两个信号之间的相关性。我们可以使用互相关来测量两个信号之间的时间延迟。NumPy 提供correlate()函数，计算两个数组之间的互相关。

怎么做...

进口情况如下:

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

加载数据，计算瞬时相位:

df = dl.data.Weather.load().dropna()
df = dl.ts.groupby_yday(df).mean().dropna()
ws_phase = dl.ts.instant_phase(df['WIND_SPEED'])
wd_phase = dl.ts.instant_phase(df['WIND_DIR'])

绘制风向和风速 z 分数:

sp = dl.plotting.Subplotter(2, 2, context)
cp = dl.plotting.CyclePlotter(sp.ax)
cp.plot(df.index, dl.stats.zscores(df['WIND_DIR'].values),
       label='Wind direction')
cp.plot(df.index, dl.stats.zscores(df['WIND_SPEED'].values),
       label='Wind speed')
sp.label()

将瞬时相位绘制如下:

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(df.index, ws_phase, label='Wind speed')
cp.plot(df.index, wd_phase, label='Wind direction')
sp.label()

绘制风速和风向的相关性:

sp.label(advance=True)
sp.ax.plot(np.correlate(df['WIND_SPEED'], df['WIND_DIR'], 'same'))

用快速傅里叶变换绘制相移图:

sp.label(advance=True)
sp.ax.plot(np.angle(np.fft.fft(df['WIND_SPEED'])/np.fft.fft(df['WIND_DIR'])))
HTML(sp.exit())

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

How to do it...

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

另见

en.wikipedia.org/wiki/Instan… 年 9 月检索)的瞬时阶段维基百科页面
en.wikipedia.org/wiki/Analyt… 的分析信号维基百科页面(2015 年 9 月检索)
en.wikipedia.org/wiki/Cross-… 的互相关的维基百科页面(2015 年 9 月检索)
位于的 angle()功能文档
功能的文档

指数平滑

指数平滑是一个低通滤波器，旨在去除噪声。在这个配方中，我们将应用单指数平滑和双指数平滑，如下式所示:

Exponential smoothing

单指数平滑(6.3)需要 平滑因子 α ，其中 0 < α < 1 。双指数平滑(6.4 和 6.5)试图通过 趋势平滑因子 β 处理数据中的趋势，其中 0 < β < 1 。

我们还将研究风速的滚动偏差，这类似于 z 分数，但它们适用于滚动窗口。平滑与回归相关，虽然平滑的目标是去噪。然而，与回归相关的指标，如均方误差 ( 均方误差)也适用于平滑。

怎么做...

进口情况如下:

import dautil as dl
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from IPython.display import HTML

定义以下函数，帮助可视化双指数平滑的结果:

def grid_mse(i, j, devs):
    alpha = 0.1 * i
    beta = 0.1 * j
    cell = dl.ts.double_exp_smoothing(devs.values, alpha, beta)

    return dl.stats.mse(devs, cell)

加载风速数据，计算年平均值和滚动偏差:

wind = dl.data.Weather.load()['WIND_SPEED'].dropna()
wind = dl.ts.groupby_year(wind).mean()
devs = dl.ts.rolling_deviations(wind, 12).dropna()

绘制风速数据的年平均值:

sp = dl.plotting.Subplotter(2, 2, context)
sp.label(ylabel_params=dl.data.Weather.get_header('WIND_SPEED'))
sp.ax.plot(wind.index, wind)

用 0.7 的 α 绘制滚动偏差:

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(devs.index, devs, label='Rolling Deviations')
cp.plot(devs.index, dl.ts.exp_smoothing(devs.values, 0.7), label='Smoothing')
sp.label()

绘制不同平滑因子的均方误差:

alphas = 0.01 * np.arange(1, 100)
errors = [dl.stats.mse(devs, dl.ts.exp_smoothing(devs.values, alpha)
          for alpha in alphas]
sp.label(advance=True)
sp.ax.plot(alphas, errors)

绘制 α 和 β 值网格的均方误差:

sp.label(advance=True)
rng = range(1, 10)
df = dl.report.map_grid(rng, rng, ["alpha", "beta", "mse"], grid_mse, devs) 
sns.heatmap(df, cmap='Blues', square=True, annot=True, fmt='.1f',
            ax=sp.ax)

HTML(sp.exit())

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

How to do it...

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

另见

en.wikipedia.org/wiki/Expone… 指数平滑的维基百科页面(2015 年 9 月检索)

评估平滑度

平滑的很多方面都堪比回归；因此，您也可以将第 10 章、评估分类器、回归器和聚类中的一些技术应用于平滑。在本食谱中，我们将使用Savitzky-Golay 滤波器进行平滑，该滤波器符合以下等式:

Evaluating smoothing

该滤波器将大小为 n 的滚动窗口内的点拟合为一个数量级为 m 的多项式。亚伯拉罕·萨维奇基和马塞尔·戈雷在 1964 年左右创造了这种算法，并首次将其应用于化学问题。过滤器有两个自然形成网格的参数。就像在回归问题中一样，我们将看一看差异，在这种情况下，原始信号和平滑信号之间的差异。我们假设，就像我们拟合数据一样，残差是随机的，并且遵循高斯分布。

怎么做...

以下步骤来自本书代码包中的eval_smooth.ipynb文件:

进口情况如下:

import dautil as dl
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
import pandas as pd
import numpy as np
import seaborn as sns
from IPython.display import HTML

定义以下助手函数:

def error(data, fit):
    return data - fit

def win_rng():
    return range(3, 25, 2)

def calc_mape(i, j, pres):
    return dl.stats.mape(pres, savgol_filter(pres, i, j))

加载大气压力数据如下:

pres = dl.data.Weather.load()['PRESSURE'].dropna()
pres = pres.resample('A')

绘制原始数据和具有窗口大小 11 和各种多项式阶次的过滤器:

sp = dl.plotting.Subplotter(2, 2, context)
cp = dl.plotting.CyclePlotter(sp.ax)
cp.plot(pres.index, pres, label='Pressure')
cp.plot(pres.index, savgol_filter(pres, 11, 2), label='Poly order 2')
cp.plot(pres.index, savgol_filter(pres, 11, 3), label='Poly order 3')
cp.plot(pres.index, savgol_filter(pres, 11, 4), label='Poly order 4')
sp.label(ylabel_params=dl.data.Weather.get_header('PRESSURE'))

绘制不同窗口大小的滤波器残差的标准偏差:

cp = dl.plotting.CyclePlotter(sp.next_ax())
stds = [error(pres, savgol_filter(pres, i, 2)).std()
        for i in win_rng()]
cp.plot(win_rng(), stds, label='Filtered')
stds = [error(pres, pd.rolling_mean(pres, i)).std()
        for i in win_rng()]
cp.plot(win_rng(), stds, label='Rolling mean')
sp.label()

绘制滤波残差的箱线图:

sp.label(advance=True)
sp.ax.boxplot([error(pres, savgol_filter(pres, i, 2))
            for i in win_rng()])
sp.ax.set_xticklabels(win_rng())

绘制窗口大小和多项式阶数网格的 MAPE 图:

sp.label(advance=True)
df = dl.report.map_grid(win_rng()[1:], range(1, 5),
                 ['win_size', 'poly', 'mape'], calc_mape, pres)
sns.heatmap(df, cmap='Blues', ax=sp.ax)
HTML(sp.exit())

最终结果参见以下截图:

How to do it...

另见

维基百科页面关于 en.wikipedia.org/wiki/Savitz… 的萨维奇基-戈莱过滤器(2015 年 9 月检索)
savgol_filter()功能记录在https://docs . scipy . org/doc/scipy/reference/generated/scipy . signal . sav gol _ filter . html(2015 年 9 月检索)

使用隆布-斯卡格尔周期图

隆布-斯卡格尔周期图是一种对数据进行正弦拟合的频谱估计方法，经常用于采样不均匀的数据。方法以尼古拉斯·r·隆布和杰弗里·d·斯卡格尔的名字命名。该算法发表于 1976 年左右，此后一直在改进。Scargle 引入了一个时间延迟参数，用于分离正弦和余弦波形。以下等式定义了时间延迟(6.7)和周期图(6.8)。

Using the Lomb-Scargle periodogram

怎么做...

进口情况如下:

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

加载太阳黑子数据如下:

df = sm.datasets.sunspots.load_pandas().data
sunspots = df['SUNACTIVITY'].values
size = len(sunspots)
t = np.linspace(-2 * np.pi, 2 * np.pi, size)
sine = dl.ts.sine_like(sunspots)
f = np.linspace(0.01, 2, 10 * size)

绘制如下正弦波形:

sp = dl.plotting.Subplotter(2, 2, context)
sp.ax.plot(t, sine)
sp.label()

sp.next_ax().plot(df['YEAR'], df['SUNACTIVITY'])
sp.label()

将周期图应用于正弦:

pgram = signal.lombscargle(t, sine, f)
sp.next_ax().plot(f, 2 * np.sqrt(pgram/size))
sp.label()

将周期图应用于太阳黑子数据:

pgram = signal.lombscargle(np.arange(size, dtype=float), sunspots, f)
sp.next_ax().plot(f, 2 * np.sqrt(pgram/size))
sp.label()
HTML(sp.exit())

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

How to do it...

前面的代码是本书代码包中lomb_scargle.ipynb文件的分解。

另见

相关维基百科页面位于(2015 年 9 月检索)
lombscargle()功能记录在https://docs . scipy . org/doc/scipy/reference/generated/scipy . signal . lombscargle . html(2015 年 9 月检索)

分析音频频谱

我们可以应用许多技术来分析音频，因此，我们可以详细讨论哪些技术最合适。最明显的方法据称是快速傅立叶变换。作为变体，我们可以使用 短时傅立叶变换 ( STFT )。STFT 将时域中的信号分成相等的部分，然后对每个部分应用快速傅立叶变换。我们将使用的另一种算法是倒谱，最初用于分析地震，但后来成功应用于语音分析。功率倒谱由下式给出:

Analyzing the frequency spectrum of audio

算法如下:

计算傅里叶变换。
计算变换的平方幅度。
取上一个结果的对数。
应用傅里叶逆变换。
再次计算大小的平方。

通常，当我们在频域中有大的变化时，倒谱是有用的。倒谱的一个重要用途是形成音频分类的特征向量。这需要从频率到 mel 标度的映射(参考中提到的维基百科页面，另请参见部分)。

怎么做...

进口情况如下:

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

定义以下功能，通过快速傅立叶变换计算信号幅度:
```
def amplitude(arr):
    return np.abs(np.fft.fft(arr))
```
加载数据如下:
```
rate, audio = read_wav()
```

绘制音频波形:

sp = dl.plotting.Subplotter(2, 2, context)
t = np.arange(0, len(audio)/float(rate), 1./rate)
sp.ax.plot(t, audio)
freqs = np.fft.fftfreq(audio.size, 1./rate)
indices = np.where(freqs > 0)[0]
sp.label()

绘制振幅谱:

magnitude = amplitude(audio)
sp.next_ax().semilogy(freqs[indices], magnitude[indices])
sp.label()

将倒谱绘制如下:

cepstrum = dl.ts.power(np.fft.ifft(np.log(magnitude ** 2)))
sp.next_ax().semilogy(cepstrum)
sp.label()

将 STFT 绘制成等高线图:

npieces = 200
stft_amps = []

for i, c in enumerate(dl.collect.chunk(audio[: npieces ** 2], len(audio)/npieces)):
    amps = amplitude(c)
    stft_amps.extend(amps)

stft_freqs = np.linspace(0, rate, npieces)
stft_times = np.linspace(0, len(stft_amps)/float(rate), npieces)
sp.next_ax().contour(stft_freqs/rate, stft_freqs,
           np.log(stft_amps).reshape(npieces, npieces))
sp.label()

HTML(sp.exit())

最终结果参见以下截图:

How to do it...

示例代码在本书代码包的analyzing_audio.ipynb文件中。

另见

维基百科关于 en.wikipedia.org/wiki/Short-… STFT 的页面(2015 年 9 月检索)
维基百科关于 en.wikipedia.org/wiki/Cepstr… T2 倒谱的页面(2015 年 9 月检索)
en.wikipedia.org/wiki/Mel_sc… 的 mel 量表的维基百科页面(2015 年 9 月检索)

用离散余弦变换分析信号

离散余弦变换 ( 离散余弦变换)是一种类似于傅里叶变换的变换，但它试图仅通过余弦项的和来表示信号(参见等式 6.11)。离散余弦变换用于信号压缩和计算 mel 频率频谱，这是我在分析音频频谱配方中提到的。我们可以通过以下等式将正常频率转换为 mel 频率(更适合分析语音和音乐的频率):

Analyzing signals with the discrete cosine transform

创建 mel 频谱的步骤并不复杂，但有很多步骤。相关的维基百科页面可在en.wikipedia.org/wiki/Mel-fr…获得(2015 年 9 月检索)。如果你做一个快速的网络搜索，你可以找到几个实现该算法的 Python 库。我在这个配方中实现了一个非常简单的计算版本。

怎么做...

进口情况如下:

import dautil as dl
from scipy.fftpack import dct
import matplotlib.pyplot as plt
import ch6util
import seaborn as sns
import numpy as np
from IPython.display import HTML

加载数据并转换如下:

rate, audio = ch6util.read_wav()
transformed = dct(audio)

使用离散余弦变换绘制振幅频谱:

sp = dl.plotting.Subplotter(2, 2, context)
freqs = np.fft.fftfreq(audio.size, 1./rate)
indices = np.where(freqs > 0)[0]
sp.ax.semilogy(np.abs(transformed)[indices])
sp.label()

绘制振幅的分布:

sns.distplot(np.log(np.abs(transformed)), ax=sp.next_ax())
sp.label()

绘制相位分布图:

sns.distplot(np.angle(transformed), ax=sp.next_ax())
sp.label()

将 mel 振幅谱绘制如下:

magnitude = ch6util.amplitude(audio)
cepstrum = dl.ts.power(np.fft.ifft(np.log(magnitude ** 2)))
mel = 1127 * np.log(1 + freqs[indices]/700)
sp.next_ax().plot(mel, ch6util.amplitude(dct(np.log(magnitude[indices] ** 2))))
sp.label()
HTML(sp.exit())

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

How to do it...

这个食谱的代码在这个书的代码包中的analyzing_dct.ipynb文件中。

另见

en.wikipedia.org/wiki/Mel_sc… 的 mel 量表的维基百科页面(2015 年 9 月检索)
en.wikipedia.org/wiki/Discre… 的离散余弦变换的维基百科页面(2015 年 9 月检索)

块自举时间序列数据

通常的自举方法不能保持时间序列数据的顺序，因此不适合趋势估计。在块自举方法中，我们将数据分割成大小相等的非重叠块，并使用这些块生成新样本。在这个食谱中，我们将应用一个非常简单且易于实现的线性模型，其中包含年温度数据。该配方的程序如下:

将数据分割成块并生成新的数据样本。
将数据拟合到一行，或者计算新数据的第一个差异。
重复上一步，建立第一个差异的斜率或中间值列表。

怎么做...

进口情况如下:

import dautil as dl
import random
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns
import ch6util
from IPython.display import HTML

定义以下函数引导数据:

def shuffle(temp, blocks):
    random.shuffle(blocks)
    df = pd.DataFrame({'TEMP': dl.collect.flatten(blocks)},
                      index=temp.index)
    df = df.resample('A')

    return df

加载数据并从中创建块:

temp = dl.data.Weather.load()['TEMP'].resample('M').dropna()
blocks = list(dl.collect.chunk(temp.values, 100))
random.seed(12033)

绘制一些随机实现作为健全性检查:

sp = dl.plotting.Subplotter(2, 2, context)
cp = dl.plotting.CyclePlotter(sp.ax)
medians = []
slopes = []

for i in range(240):
    df = shuffle(temp, blocks)
    slopes.append(ch6util.fit(df))
    medians.append(ch6util.diff_median(df))

    if i < 5:
        cp.plot(df.index, df.values)

sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

使用自举数据绘制第一个差值中位数的分布:

sns.distplot(medians, ax=sp.next_ax(), norm_hist=True)
sp.label()

使用自举数据

sns.distplot(slopes, ax=sp.next_ax(), norm_hist=True)
sp.label()

绘制线性回归斜率的分布

绘制不同数量引导的置信区间:

mins = []
tops = []
xrng = range(30, len(medians))

for i in xrng:
    min, max = dl.stats.outliers(medians[:i])
    mins.append(min)
    tops.append(max)

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(xrng, mins, label='5 %')
cp.plot(xrng, tops, label='95 %')
sp.label()
HTML(sp.exit())

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

How to do it...

以下代码来自本书代码包中的block_boot.ipynb文件。

另见

相关维基百科页面位于https://en . Wikipedia . org/wiki/Bootstrapping _ % 28 statistics % 29 # Block _ Bootstrapping(2015 年 9 月检索)

移动块自举时间序列数据

如果你按照块自举时间序列数据的方法，你现在知道了一个简单的时间序列数据自举方案。移动块自举算法有点复杂。在该方案中，我们通过移动固定大小的窗口来生成重叠块，类似于移动平均。然后，我们组装这些块来创建新的数据样本。

在本食谱中，我们将把移动块自举应用于年温度数据，以生成第二差值中位数和 AR(1)模型斜率的列表。这是一个滞后为 1 的自回归模型。此外，我们将尝试使用中值滤波器来中和异常值和噪声。

怎么做...

以下代码片段来自本书代码包中的moving_boot.ipynb文件:

进口情况如下:

import dautil as dl
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import ch6util
from scipy.signal import medfilt
from IPython.display import HTML

定义以下函数引导数据:

def shuffle(temp):
    indices = np.random.choice(start, n/12)
    sample = dl.collect.flatten([temp.values[i: i + 12] for i in indices])
    sample = medfilt(sample)
    df = pd.DataFrame({'TEMP': sample}, index=temp.index[:len(sample)])
    df = df.resample('A', how=np.median)

    return df

加载数据如下:

temp = dl.data.Weather.load()['TEMP'].resample('M', how=np.median).dropna()
n = len(temp)
start = np.arange(n - 11)
np.random.seed(2609787)

绘制一些随机实现作为健全性检查:

sp = dl.plotting.Subplotter(2, 2, context)
cp = dl.plotting.CyclePlotter(sp.ax)
medians = []
slopes = []

for i in range(240):
    df = shuffle(temp)
    slopes.append(dl.ts.ar1(df.values.flatten())['slope'])
    medians.append(ch6util.diff_median(df, 2))

    if i < 5:
        cp.plot(df.index, df.values)

sp.label(ylabel_params=dl.data.Weather.get_header('TEMP'))

使用自举数据绘制第二差值中位数的分布:
```
sns.distplot(medians, ax=sp.next_ax())
sp.label()
```
使用自举数据绘制 AR(1)模型斜率的分布:
```
sns.distplot(slopes, ax=sp.next_ax())
sp.label()
```

绘制不同引导次数的置信区间:

mins = []
tops = []
xrng = range(30, len(medians))

for i in xrng:
    min, max = dl.stats.outliers(medians[:i])
    mins.append(min)
    tops.append(max)

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(xrng, mins, label='5 %')
cp.plot(xrng, tops, label='95 %')
sp.label()
HTML(sp.exit())

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

How to do it...

另见

相关维基百科页面位于https://en . Wikipedia . org/wiki/Bootstrapping _ % 28 statistics % 29 # Block _ Bootstrapping(2015 年 9 月检索)
位于的medfilt()文件 https://docs . scipy . org/doc/scipy/reference/generated/scipy . signal . med filt . html(检索【2015 年 9 月)

应用离散小波变换

离散小波变换 ( 离散小波变换)捕获时域和频域的信息。数学家阿尔弗雷德·哈尔创造了第一个小波。我们在这个食谱中也会用到这个 哈尔小波。变换返回近似和细节系数，我们需要一起使用才能得到原始信号。逼近系数是低通滤波器的结果。高通滤波器产生细节系数。哈尔小波算法的阶数为 0(n)，类似于 STFT 算法(参考分析音频频谱配方)，结合了频率和时间信息。

与傅里叶变换的不同之处在于，我们将信号表示为正弦和余弦项的和，而小波由单个波(小波函数)表示。就像在 STFT 一样，我们在时域中分割信号，然后将小波函数应用于每个片段。DWT 在这个配方中可以有多个级别，我们不会超过第一个级别。为了获得下一级，我们将小波应用于前一级的近似系数。这意味着我们可以有多个层次的细节系数。

作为数据集，我们将看一看著名的尼罗河流量，这甚至是希腊历史学家希罗多德写的。最近，在上个世纪，水文学家赫斯特发现了一个幂律，适用于当年尼罗河 T2 流量的重新标度范围。更多信息请参考参见部分。重新缩放的范围并不难计算，但有许多步骤，如下式所述:

Applying the discrete wavelet transform

来自权力定律的赫斯特指数是趋势的指示器。我们还可以用更有效的方法从小波系数中得到赫斯特指数。

开始

安装pywavelets，如下:

$ pip install pywavelets

这个食谱我用的是pywavelets 0.3.0。

怎么做...

进口情况如下:

from statsmodels import datasets
import matplotlib.pyplot as plt
import pywt
import pandas as pd
import dautil as dl
import numpy as np
import seaborn as sns
import warnings
from IPython.display import HTML

过滤警告如下(可选步骤):

warnings.filterwarnings(action='ignore',
                        message='.*Mean of empty slice.*')
warnings.filterwarnings(action='ignore',
                        message='.*Degrees of freedom <= 0 for slice.*')

定义以下函数来计算重新缩放的范围:

def calc_rescaled_range(X):
    N = len(X)

    # 1\. Mean
    mean = X.mean()

    # 2\. Y mean adjusted
    Y = X - mean

    # 3\. Z cumulative deviates
    Z = np.array([Y[:i].sum() for i in range(N)])

    # 4\. Range R
    R = np.array([0] + [np.ptp(Z[:i]) for i in range(1, N)])

    # 5\. Standard deviation S
    S = np.array([X[:i].std() for i in range(N)])

    # 6\. Average partial R/S
    return [np.nanmean(R[:i]/S[:i]) for i in range(N)]

加载数据并用哈尔小波变换:

data = datasets.get_rdataset('Nile', cache=True).data
cA, cD = pywt.dwt(data['Nile'].values, 'haar')
coeff = pd.DataFrame({'cA': cA, 'cD': cD})

绘制尼罗河流量如下图:

sp = dl.plotting.Subplotter(2, 2, context)
sp.ax.plot(data['time'], data['Nile'])
sp.label()

绘制转换数据的近似值和细节系数:

cp = dl.plotting.CyclePlotter(sp.next_ax())
cp.plot(range(len(cA)), cA, label='Approximation coefficients')
cp.plot(range(len(cD)), cD, label='Detail coefficients')
sp.label()

绘制系数的重新缩放的范围如下:

sp.next_ax().loglog(range(len(cA)), calc_rescaled_range(cA), 
                    label='Approximation coefficients')
sp.ax.loglog(range(len(cD)), calc_rescaled_range(cD), 
             label='Detail coefficients')
sp.label()

拟合绘制尼罗河流量数据的重新标度范围:

range_df = pd.DataFrame(data={'Year': data.index,
                              'Rescaled range':
                              calc_rescaled_range(data['Nile'])})
sp.next_ax().set(xscale="log", yscale="log")
sns.regplot('Year', 'Rescaled range', range_df, ax=sp.ax, order=1,
            scatter_kws={"s": 100})
sp.label()
HTML(sp.exit())

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

How to do it...

相关代码在本书的代码包discrete_wavelet.ipynb文件中。

另见

en.wikipedia.org/wiki/Discre… 离散小波变换的维基百科页面(2015 年 9 月检索)
维基百科关于 en.wikipedia.org/wiki/Rescal… 的重新标度范围的页面(2015 年 9 月检索)
维基百科关于 en.wikipedia.org/wiki/Hurst_… 赫斯特指数的页面(2015 年 9 月检索)