Clojure 数据科学（五）
原文：
annas-archive.org/md5/bd3fb3fcef82cc10763b55da27a85129译者：飞龙
第九章 时间序列
| “再次经过的时间。” | ||
|---|---|---|
| --卡罗琳·基恩，《老钟的秘密》 |
在之前的几章中，我们看到如何应用迭代算法来解决复杂方程。我们首先在梯度下降法中遇到这一点——包括批量和随机梯度下降——但最近我们在使用图并行计算模型进行图社区检测时也看到了这一点。
本章讨论的是时间序列数据。时间序列是指按测量时间排列的某一数量的定期观测数据系列。为了使本章中的许多技术能够工作，我们要求连续观测之间的间隔是相等的。测量间隔可以是销售数据的月度、降雨量或股市波动的日度，或者高流量网站的访问量以分钟为单位。
为了能够预测时间序列的未来值，我们需要未来的值在某种程度上依赖于之前的值。因此，本章也涉及到递归：我们如何构建一个序列，其中每个新值都是前一个值的函数。通过将实际时间序列建模为一个新值以这种方式生成的过程，我们希望能够将序列向前推演，并做出预测。
不过，在我们开始讨论递归之前，我们将学习如何利用已经接触过的一种技术——线性回归——来拟合时间序列数据的曲线。
关于数据
本章将使用两个 Incanter 预装的数据集：Longley 数据集，包含了 1947 年到 1962 年间美国七个经济变量的数据，以及Airline 数据集，包含了 1949 年 1 月到 1960 年 12 月的每月航空乘客总数数据。
注意
你可以从github.com/clojuredatascience/ch9-time-series下载本章的源代码。
Airline 数据集是我们在本章中将花费大部分时间的部分，但首先让我们看看 Longley 数据集。它包含的列包括国内生产总值（GDP）、就业与失业人口数量、人口和军队规模。由于许多预测变量本身是相关的，这使得它成为分析多重共线性的经典数据集。幸运的是，这不会影响我们的分析，因为我们一次只使用一个预测变量。
加载 Longley 数据
由于 Incanter 将 Longley 数据集作为其示例数据集库的一部分，因此加载数据只需调用incanter.datasets/get-dataset并将:longley作为唯一参数。一旦加载完毕，我们可以通过incanter.core/view来查看数据集：
(defn ex-9-1 []
(-> (d/get-dataset :longley)
(i/view)))
数据应大致如下所示：
这些数据最初由国家标准与技术研究院（NIST）发布，作为统计参考数据集，列的描述可以在他们的网站上找到：www.itl.nist.gov/div898/strd/lls/data/LINKS/i-Longley.shtml。我们将考虑最后三列：x4：武装部队的规模，x5：14 岁及以上的“非机构”人口，以及x6：年份。
首先，让我们看看人口随时间的变化：
(defn ex-9-2 []
(let [data (d/get-dataset :longley)]
(-> (c/scatter-plot (i/$ :x6 data)
(i/$ :x5 data)
:x-label "Year"
:y-label "Population")
(i/view))))
前面的代码生成了以下图表：
年份与人口的图表显示出一种非常明显的近似线性关系。轻微的曲线表明，随着人口增加，人口增长率也在上升。
注意
回想一下第三章中的吉布拉特法则，相关性，公司增长率与其规模成正比。在分析适用吉布拉特法则的人口时，常常会看到类似前面图表的增长曲线：增长率会随着时间的推移而增加。
我们已经看到如何使用 Incanter 的线性模型拟合直线。或许令人惊讶的是，linear-model函数同样可以用来拟合曲线。
使用线性模型拟合曲线
首先，让我们回顾一下如何使用 Incanter 的linear-model函数拟合直线。我们要从数据集中提取x5和x6列，并将它们（按顺序：x6，年份，是我们的预测变量）应用于incanter.stats/linear-model函数。
(defn ex-9-3 []
(let [data (d/get-dataset :longley)
model (s/linear-model (i/$ :x5 data)
(i/$ :x6 data))]
(println "R-square" (:r-square model))
(-> (c/scatter-plot (i/$ :x6 data)
(i/$ :x5 data)
:x-label "Year"
:y-label "Population")
(c/add-lines (i/$ :x6 data)
(:fitted model))
(i/view))))
;; R-square 0.9879
前面的代码生成了以下图表：
尽管直线拟合数据的效果相当好，R²值超过 0.98，但它并未捕捉到线的曲线。特别是，我们可以看到图表两端以及中间的点偏离了直线。我们简单的模型有很高的偏差，且根据年份的不同，系统性地低估和高估了人口。残差的图表清楚地显示，误差并不是均匀分布的。
linear-model函数之所以被称为线性模型，是因为它生成的模型与其参数之间存在线性关系。然而，或许出乎意料的是，只要我们提供非线性特征，它也能够生成非线性预测。例如，我们可以在年份之外，添加年份的平方作为一个参数。在以下代码中，我们使用 Incanter 的bind-columns函数将这两个特征结合起来，形成一个矩阵：
(defn ex-9-4 []
(let [data (d/get-dataset :longley)
x (i/$ :x6 data)
xs (i/bind-columns x (i/sq x))
model (s/linear-model (i/$ :x5 data) xs)]
(println "R-square" (:r-square model))
(-> (c/scatter-plot (i/$ :x6 data)
(i/$ :x5 data)
:x-label "Year"
:y-label "Population")
(c/add-lines (i/$ :x6 data)
(:fitted model))
(i/view))))
;; 0.9983
我们的R²值有所提高，得到以下图表：
这显然更适合数据。我们可以通过创建一个 forecast 函数来使用我们的模型进行预测，该函数接受模型的系数并返回一个关于 x（年份）的函数，将系数乘以我们定义的特征：
(defn forecast [coefs]
(fn [x]
(first
(i/mmult (i/trans coefs)
(i/matrix [1.0 x (i/sq x)])))))
系数包括了一个偏置项参数，因此我们将系数乘以 1.0、x 和 x²。
(defn ex-9-5 []
(let [data (d/get-dataset :longley)
x (i/$ :x6 data)
xs (i/bind-columns x (i/sq x))
model (s/linear-model (i/$ :x5 data) xs)]
(-> (c/scatter-plot (i/$ :x6 data)
(i/$ :x5 data)
:x-label "Year"
:y-label "Population")
(c/add-function (forecast (:coefs model))
1947 1970)
(i/view))))
接下来，我们将函数图扩展到 1970 年，以更清晰地看到拟合模型的曲线，如下所示：
当然，我们正在进行超出数据范围的外推。如同在第三章《相关性》中讨论的那样，外推太远通常是不明智的。为了更清楚地说明这一点，让我们关注一下 Longley 数据集中的另一个列，即武装力量的规模：x6。
我们可以像之前一样绘制它：
(defn ex-9-6 []
(let [data (d/get-dataset :longley)]
(-> (c/scatter-plot (i/$ :x6 data)
(i/$ :x4 data)
:x-label "Year"
:y-label "Size of Armed Forces")
(i/view))))
这生成了以下图表：
这显然是一个更复杂的序列。我们可以看到，在 1950 年到 1952 年间，武装力量的规模急剧增加，随后出现了缓慢下降。在 1950 年 6 月 27 日，杜鲁门总统命令空军和海军支援韩国，这成为了后来所称的朝鲜战争。
要拟合这些数据的曲线，我们需要生成高阶多项式。首先，让我们构建一个 polynomial-forecast 函数，它将基于单一的 x 和要创建的最高次多项式自动生成更高阶的特征：
(defn polynomial-forecast [coefs degree]
(fn [x]
(first
(i/mmult (i/trans coefs)
(for [i (range (inc degree))]
(i/pow x i))))))
例如，我们可以使用以下代码训练一个模型，直到 x11：
(defn ex-9-7 []
(let [data (d/get-dataset :longley)
degree 11
x (s/sweep (i/$ :x6 data))
xs (reduce i/bind-columns
(for [i (range (inc degree))]
(i/pow x i)))
model (s/linear-model (i/$ :x4 data) xs
:intercept false)]
(println "R-square" (:r-square model))
(-> (c/scatter-plot (i/$ 1 xs) (i/$ :x4 data)
:x-label "Distance from Mean (Years)"
:y-label "Size of Armed Forces")
(c/add-function (polynomial-forecast (:coefs model)
degree)
-7.5 7.5)
(i/view))))
;; R-square 0.9755
上述代码生成了以下图表：
曲线拟合数据相当好，R²超过 0.97。然而，现在你应该不感到惊讶，我们实际上是在过拟合数据。我们构建的模型可能没有太多的预测能力。事实上，如果我们像在 ex-9-8 中那样将图表的范围向右延伸，展示未来的预测，我们将得到以下结果：
在最后一次测量数据点后的两年半，我们的模型预测军事人数将增长超过 500%，达到超过 175,000 人。
时间序列分解
我们在建模军事时间序列时面临的一个问题是，数据量不足以生成一个能够代表产生该系列过程的通用模型。建模时间序列的常见方法是将序列分解为若干个独立的组件：
-
趋势：序列是否随时间总体上升或下降？这种趋势是否像我们在观察人口时看到的那样是指数曲线？
-
季节性：序列是否表现出周期性的涨跌，且周期数目固定？对于月度数据，通常可以观察到 12 个月的周期性循环。
-
周期：数据集中是否存在跨越多个季节的长期周期？例如，在金融数据中，我们可能会观察到与扩张和衰退周期相对应的多年度周期。
另一个指定军事数据问题的方式是，由于没有足够的信息，我们无法确定是否存在趋势，或者观察到的峰值是否属于季节性或周期性模式的一部分。尽管数据似乎有上升的趋势，但也有可能我们正在密切观察一个最终会回落的周期。
本章我们将研究的一个数据集是经典的时间序列数据，观察的是 1949 到 1960 年间的月度航空公司乘客人数。该数据集较大，并且明显展示了趋势和季节性成分。
检查航空公司数据
与 Longley 数据集类似，航空公司数据集也是 Incanter 数据集库的一部分。我们加载incanter.datasets库为d，并将incanter.code库加载为i。
(defn ex-9-9 []
(-> (d/get-dataset :airline-passengers)
(i/view)))
前几行数据应如下所示：
在分析时间序列时，确保数据按时间顺序排列是非常重要的。该数据按年份和月份排序。所有一月的数据都排在所有二月数据之前，以此类推。为了进一步处理，我们需要将年份和月份列转换为一个可以排序的单一列。为此，我们将再次使用clj-time库（github.com/clj-time/clj-time）。
可视化航空公司数据
在第三章相关性中解析时间时，我们能够利用时间的字符串表示是 clj-time 默认理解的格式。自然，clj-time 并不能自动推断所有时间表示的格式。特别是，美国格式的mm/dd/yyyy与世界大多数地方偏好的dd/mm/yyyy格式之间的差异尤为问题。在clj-time.format命名空间中提供了一个parse函数，允许我们传递一个格式字符串，指示库如何解释该字符串。在以下代码中，我们将format命名空间作为tf引入，并指定我们的时间将采用"MMM YYYY"格式。
注意
一个由 clj-time 使用的格式化字符串列表可以在www.joda.org/joda-time/key_format.html找到。
换句话说，是三个字符表示“月”，后面跟着四个字符表示“年”。
(def time-format
(tf/formatter "MMM YYYY"))
(defn to-time [month year]
(tf/parse time-format (str month " " year)))
有了之前的函数，我们可以将年份和月份列解析成一个单一的时间，按顺序排列它们，并提取乘客人数：
(defn airline-passengers []
(->> (d/get-dataset :airline-passengers)
(i/add-derived-column :time [:month :year] to-time)
(i/$order :time :asc)
(i/$ :passengers)))
结果是一组按时间顺序排列的乘客人数数据。现在让我们将其作为折线图来可视化：
(defn timeseries-plot [series]
(-> (c/xy-plot (range (count series)) series
:x-label "Time"
:y-label "Value")
(i/view)))
(defn ex-9-10 []
(-> (airline-passengers)
(timeseries-plot)))
上述代码生成了以下图表：
你可以看到数据呈现出明显的季节性模式（每 12 个月重复一次）、上升趋势和缓慢增长曲线。
图表右侧的方差大于左侧的方差，因此我们说数据表现出一定的异方差性。我们需要去除数据集中的方差增加和上升趋势。这样会得到一个平稳的时间序列。
平稳性
一个平稳的时间序列是指其统计特性随时间保持不变的序列。大多数统计预测方法假设序列已经被转换为平稳序列。对于平稳的时间序列，预测变得更加容易：我们假设该序列的统计特性在未来与过去相同。为了去除航空公司数据中的方差增长和增长曲线，我们可以简单地对乘客数量取对数：
(defn ex-9-11 []
(-> (airline-passengers)
(i/log)
(timeseries-plot)))
这生成了以下图表：
取对数的效果是双重的。首先，初始图表中显现的异方差性被去除了。其次，指数增长曲线被转化为一条直线。
这使得数据变得更容易处理，但我们仍然在序列中看到趋势，也称为漂移。为了获得真正平稳的时间序列，我们还需要使均值稳定。实现这一点有几种方法。
去趋势和差分
第一种方法是去趋势处理序列。在取对数后，航空公司数据集包含了一个非常强的线性趋势。我们可以为这些数据拟合一个线性模型，然后绘制残差图：
(defn ex-9-12 []
(let [data (i/log (airline-passengers))
xs (range (count data))
model (s/linear-model data xs)]
(-> (c/xy-plot xs (:residuals model)
:x-label "Time"
:y-label "Residual")
(i/view))))
这生成了以下图表：
残差图显示了一个均值比原始序列更稳定的序列，并且上升趋势已被完全去除。然而，不幸的是，残差似乎并未围绕新均值呈现正态分布。特别是图表中间似乎出现了一个“驼峰”。这表明我们的线性模型在航空公司数据上表现不佳。我们可以像本章开始时那样拟合一条曲线，但让我们改为查看另一种使时间序列平稳的方法。
第二种方法是差分。如果我们从时间序列中的每个点中减去其直接前一个点的值，我们将得到一个新的时间序列（少一个数据点），其中只包含相邻点之间的差异。
(defn difference [series]
(map - (drop 1 series) series))
(defn ex-9-13 []
(-> (airline-passengers)
(i/log)
(difference)
(timeseries-plot)))
我们可以在以下图表中看到这一效果：
注意，上升趋势已经被替换为围绕常数均值波动的序列。均值略高于零，这对应于差异更可能为正，并导致我们观察到的上升趋势。
这两种技术的目标都是使序列的均值保持恒定。在某些情况下，可能需要对序列进行多次差分，或者在去趋势之后应用差分，以获得一个真正均值稳定的序列。例如，在去趋势之后，序列中仍然可以看到一些漂移，因此在本章的其余部分，我们将使用差分后的数据。
在继续讨论如何为预测建模这样的时间序列之前，让我们绕个弯，思考一下什么是时间序列，以及我们如何将时间序列建模为一个递归过程。
离散时间模型
离散时间模型，如我们迄今为止所看的那些，将时间划分为定期的时间片。为了使我们能够预测未来时间片的值，我们假设它们依赖于过去的时间片。
注意
时间序列也可以根据频率而不是时间进行分析。我们在本章中不讨论频域分析，但书籍的维基页面wiki.clojuredatascience.com包含了更多资源的链接。
在以下内容中，令y[t]表示时间t时刻观察值的值。最简单的时间序列是每个时间片的值与前一个时间片的值相同。此类序列的预测器为：
这就是说，在时间t + 1的预测值给定时间t时等于时间t的观察值。请注意，这个定义是递归的：时间t的值依赖于时间t - 1的值。时间t - 1的值依赖于时间t - 2的值，依此类推。
我们可以将这个“常数”时间序列建模为 Clojure 中的惰性序列，其中序列中的每个值都是常数值：
(defn constant-series [y]
(cons y (lazy-seq (constant-series y))))
(defn ex-9-14 []
(take 5 (constant-series 42)))
;; (42 42 42 42 42)
注意constant-series的定义中包含了对其自身的引用。这是一个递归函数定义，它创建了一个无限的惰性序列，我们可以从中获取值。
下一时刻，即时间t + 1，实际观察到的值为y[t+1]。如果此值与我们的预测值不同，则我们可以将这个差异计算为预测的误差：
通过结合前两个方程，我们得到了时间序列的随机模型。
换句话说，当前时间片的值等于前一个时间片的值加上一些误差。
随机漫步
最简单的随机过程之一是随机漫步。让我们将constant-series扩展为一个random-walk过程。我们希望我们的误差服从均值为零且方差恒定的正态分布。我们将通过调用 Incanter 的stats/sample-normal函数来模拟随机噪声。
(defn random-walk [y]
(let [e (s/sample-normal 1)
y (+ y e)]
(cons y (lazy-seq (random-walk y)))))
(defn ex-9-15 []
(->> (random-walk 0)
(take 50)
(timeseries-plot)))
当然，您会得到不同的结果，但它应该类似于以下图表：
随机游走模型在金融学和计量经济学中非常常见。
注意
随机游走这一术语最早由卡尔·皮尔逊在 1905 年提出。许多过程——从波动的股价到分子在气体中运动时所描绘的路径——都可以被建模为简单的随机游走。1973 年，普林斯顿经济学家伯顿·戈登·马基尔在他的著作《华尔街随机漫步》中提出，股价也遵循随机游走的规律。
随机游走并非完全无法预测。虽然每个点与下一个点之间的差异由随机过程决定，但该差异的方差是恒定的。这意味着我们可以估计每一步的置信区间。然而，由于均值为零，我们无法有任何信心地判断差异相对于当前值是正的还是负的。
自回归模型
我们在本章中已经看到如何使用线性模型基于预测变量的线性组合做出预测。在自回归模型中，我们通过使用该变量过去的值的线性组合来预测感兴趣的变量。
自回归模型将预测变量与其自身回归。为了实际观察这一过程，让我们看一下以下代码：
(defn ar [ys coefs sigma]
(let [e (s/sample-normal 1 :sd sigma)
y (apply + e (map * ys coefs))]
(cons y (lazy-seq
(ar (cons y ys) coefs sigma)))))
这与我们在几页前遇到的随机游走递归定义有很多相似之处。然而，这次我们通过将之前的ys与coefs相乘，生成每一个新的y。
我们可以通过调用我们新的ar函数，传入之前的ys和自回归模型的系数，来生成一个自回归序列：
(defn ex-9-16 []
(->> (ar [1] [2] 0)
(take 10)
(timeseries-plot)))
这将生成以下图表：
通过取初始值 1.0 和系数 2.0，且无噪声，我们正在创建一个指数增长曲线。序列中的每一个时间步都是二的幂。
自回归序列被称为自相关的。换句话说，每个点与其前面的点是线性相关的。在之前的例子中，这只是前一个值的两倍。系数的数量被称为自相关模型的阶数，通常用字母p表示。因此，前面的例子是一个p=1的自回归过程，或者是*AR(1)*过程。
通过增加p，可以生成更复杂的自回归序列。
(defn ex-9-17 []
(let [init (s/sample-normal 5)]
(->> (ar init [0 0 0 0 1] 0)
(take 30)
(timeseries-plot))))
例如，前面的代码生成了一个 5 阶自回归时间序列，或*AR(5)*序列。效果在序列中表现为一个周期为 5 个点的规律性循环。
我们可以将自回归模型与一些噪声结合起来，引入我们之前看到的随机游走组件。让我们将 sigma 增加到 0.2：
(defn ex-9-18 []
(let [init (s/sample-normal 5)]
(->> (ar init [0 0 0 0 1] 0.2)
(take 30)
(timeseries-plot))))
这将生成以下图表：
请注意，每五个点的典型“季节性”周期是如何被保留下来的，但也与一些噪声元素结合在一起。虽然这是模拟数据，但这个简单的模型已经开始接近在时间序列分析中经常出现的系列类型。
一阶滞后 AR 模型的一般方程如下：
其中c是某个常数，ε[t]是误差，y[t-1]是序列在上一个时间步的值，φ[1]是由希腊字母phi表示的系数。更一般地，自回归模型在p个滞后下的方程为：
由于我们的序列是平稳的，因此我们在代码中省略了常数项c。
在 AR 模型中确定自相关
就像线性回归可以建立多个自变量之间的（线性）相关性一样，自回归可以建立一个变量与自身在不同时间点之间的相关性。
就像在线性回归中，我们试图建立预测变量与响应变量之间的相关性一样，在时间序列分析中，我们试图建立时间序列与自身在一定数量滞后下的自相关。知道存在自相关的滞后数量可以让我们计算出自回归模型的阶数。
由此可得，我们想要研究序列在不同滞后下的自相关。例如，零滞后意味着我们将每个点与其本身进行比较（自相关为 1.0）。滞后 1 意味着我们将每个点与直接前面的点进行比较。自相关函数（ACF）是一个数据集与自身在给定滞后下的线性依赖关系。因此，ACF 由滞后k参数化。
Incanter 包含一个auto-correlation函数，可以返回给定序列和滞后的自相关。然而，我们定义了自己的autocorrelation函数，它将返回一系列滞后的autocorrelation：
(defn autocorrelation* [series variance n k]
(let [lag-product (->> (drop k series)
(map * series)
(i/sum))]
(cons (/ lag-product variance n)
(lazy-seq
(autocorrelation* series variance n (inc k))))))
(defn autocorrelation [series]
(autocorrelation* (s/sweep series)
(s/variance series)
(dec (count series)) 0))
在计算自相关之前，我们使用incanter.stats的sweep函数来从序列中去除均值。这意味着我们可以简单地将序列中的值与滞后* k *的值相乘，以确定它们是否有一起变化的趋势。如果有，它们的乘积将是正的；如果没有，它们的乘积将是负的。
该函数返回一个无限懒加载的自相关值序列，对应于滞后0...k的自相关。我们来定义一个函数，将这些值绘制成条形图。与timeseries-plot一样，该函数将接受一个有序的值序列：
(defn bar-plot [ys]
(let [xs (range (count ys))]
(-> (c/bar-chart xs ys
:x-label "Lag"
:y-label "Value")
(i/view))))
(defn ex-9-19 []
(let [init (s/sample-normal 5)
coefs [0 0 0 0 1]]
(->> (ar init coefs 0.2)
(take 100)
(autocorrelation)
(take 15)
(bar-plot))))
这生成了以下图表：
每隔 5 个滞后的峰值与我们的*AR(5)*系列生成器一致。随着噪声对信号的干扰，峰值逐渐减弱，测得的自相关也随之下降。
移动平均模型
AR 模型的一个假设是噪声是随机的，且均值和方差恒定。我们的递归 AR 函数从正态分布中抽取值来生成满足这些假设的噪声。因此，在 AR 过程中，噪声项是相互不相关的。
然而，在某些过程中，噪声项本身并不是不相关的。以一个示例来说明，假设有一个时间序列记录了每天的烧烤销售数量。我们可能会观察到每隔 7 天就有一个峰值，代表着顾客在周末更倾向于购买烧烤。偶尔，我们可能会观察到几周的销售总量下降，另外几周销售则会相应上升。我们可能推测这是由天气造成的，差的销售额对应于寒冷或雨天的时期，而好的销售额则对应于有利天气的时期。无论是什么原因，这种现象都会以均值明显变化的形式出现在我们的数据中。展现此行为的系列被称为移动平均（MA）模型。
一阶移动平均模型，记作MA(1)，为：
其中，μ 是系列的均值，ε[t] 是噪声值，θ[1] 是模型的参数。更一般地，对于q项，MA 模型可表示为：
因此，MA 模型本质上是当前系列值与当前及前一时刻（未观测到的）白噪声误差项或随机冲击的线性回归。假设每个点的误差项是相互独立的，并且来自同一（通常为正态）分布，具有零均值和恒定方差。
在 MA 模型中，我们假设噪声值本身是自相关的。我们可以这样建模：
(defn ma [es coefs sigma]
(let [e (s/sample-normal 1 :sd sigma)
y (apply + e (map * es coefs))]
(cons y (lazy-seq
(ma (cons e es) coefs sigma)))))
这里，es 是前一时刻的误差，coefs 是 MA 模型的参数，sigma 是误差的标准差。
注意该函数与之前定义的ar函数的不同。我们不再保留前一系列的ys，而是保留前一系列的es。接下来，让我们看看 MA 模型生成的系列：
(defn ex-9-20 []
(let [init (s/sample-normal 5)
coefs [0 0 0 0 1]]
(->> (ma init coefs 0.5)
(take 100)
(timeseries-plot))))
这会生成类似于以下的图形：
你可以看到，MA 模型的图表没有 AR 模型那样明显的重复。如果从更长的时间序列来看，你可以看到它如何重新引入漂移，因为一个随机冲击的反响在新的临时均值中得以延续。
确定 MA 模型中的自相关
你可能会想，是否可以通过自相关图来帮助识别 MA 过程。我们现在就来绘制一下。MA 模型可能更难识别，因此我们将在绘制自相关之前生成更多数据点。
(defn ex-9-21 []
(let [init (s/sample-normal 5)
coefs [0 0 0 0 1]]
(->> (ma init coefs 0.2)
(take 5000)
(autocorrelation)
(take 15)
(bar-plot))))
这将生成以下图表：
你可以在前面的图表中看到，MA 过程的阶数被清晰地显示出来，在滞后 5 的位置有一个明显的峰值。请注意，然而，与自回归过程不同，这里没有每 5 个滞后就出现一个周期性的峰值。这是 MA 过程的一个特征，因为该过程引入了漂移至均值，自相关在其他滞后下大幅减弱。
结合 AR 和 MA 模型
到目前为止，我们所考虑的 AR 和 MA 模型是生成自相关时间序列的两种不同但紧密相关的方法。它们并不是互相排斥的，实际上在建模实际时间序列时，你经常会遇到序列看起来是两者的混合情况。
我们可以将 AR 和 MA 过程结合成一个单一的 ARMA 模型，包含两组系数：自回归模型的系数和移动平均模型的系数。每个模型的系数数量不必相同，按照惯例，AR 模型的阶数用p表示，而 MA 模型的阶数用q表示。
(defn arma [ys es ps qs sigma]
(let [e (s/sample-normal 1 :sd sigma)
ar (apply + (map * ys ps))
ma (apply + (map * es qs))
y (+ ar ma e)]
(cons y (lazy-seq
(arma (cons y ys)
(cons e es)
ps qs sigma)))))
让我们绘制一个更长的数据点序列，看看会出现什么样的结构：
(defn ex-9-22 []
(let [ys (s/sample-normal 10 :sd 1.0)
es (s/sample-normal 10 :sd 0.2)
ps [0 0 0 0.3 0.5]
qs [0.2 0.8]]
(->> (arma ys es ps qs 0.2)
(take 500)
(timeseries-plot))))
注意我们如何为模型的 AR 部分和 MA 部分指定不同数量的参数：AR 部分为 5 个参数，MA 模型为 2 个参数。这被称为*ARMA(5,2)*模型。
之前的 ARMA 模型在更长的点序列上绘制，可以让 MA 项的效果变得明显。在这个尺度下，我们看不到 AR 成分的效果，因此让我们像之前一样通过自相关图来运行这个序列：
(defn ex-9-23 []
(let [ys (s/sample-normal 10 :sd 1.0)
es (s/sample-normal 10 :sd 0.2)
ps [0 0 0 0.3 0.5]
qs [0.2 0.8]]
(->> (arma ys es ps qs 0.2)
(take 500)
(autocorrelation)
(take 15)
(bar-plot))))
你应该看到一个类似以下的图表：
事实上，随着更多数据和 AR 与 MA 成分的加入，ACF 图并没有让序列的阶数变得更加清晰，相反，它不像我们之前看到的那些自相关图那么明确和清晰。自相关缓慢衰减至零，这使得我们即使提供了大量数据，也无法确定 AR 和 MA 过程的阶数。
这样做的原因是 MA 部分的影响压倒了 AR 部分。我们无法在数据中识别出周期性模式，因为它被一个移动平均所掩盖，这使得所有相邻的点看起来都存在相关性。解决这个问题的最佳方法是绘制偏自相关。
计算偏自相关
部分自相关函数（PACF）旨在解决在混合 ARMA 模型中识别周期性成分的问题。它被定义为给定所有中间观测值时，y[t] 和 y[t+k] 之间的相关系数。换句话说，它是在滞后 k 下的自相关，且该自相关并未被滞后 1 到滞后 k-1 的自相关所涵盖。
一阶滞后 1 的部分自相关被定义为等于一阶自相关。二阶滞后 2 的部分自相关等于：
这是两个时间周期间隔的值之间的相关性，y[t] 和 y[t-2]，条件是已知 y[t-1]。在平稳时间序列中，分母中的两个方差将是相等的。
三阶滞后 3 的部分自相关等于：
以此类推，适用于任何滞后。
自协方差
部分自相关的公式要求我们计算数据在某个滞后下的自协方差。这称为自协方差。我们在前几章已经学习了如何测量两组数据之间的协方差，也就是两项或多项属性一起变化的趋势。这个函数与我们之前在本章定义的自相关函数非常相似，计算从零开始的一系列滞后的自协方差：
(defn autocovariance* [series n k]
(let [lag-product (->> (drop k series)
(map * series)
(i/sum))]
(cons (/ lag-product n)
(lazy-seq
(autocovariance* series n (inc k))))))
(defn autocovariance [series]
(autocovariance* (s/sweep series) (count series) 0))
如前所述，返回值将是一个懒序列（lazy sequence）滞后值，所以我们只需要取我们所需的值。
使用 Durbin-Levinson 递推法的 PACF
由于需要考虑先前已解释的变化，计算部分自相关比计算自相关要复杂得多。Durbin-Levinson 算法提供了一种递归计算的方法。
注意
Durbin-Levinson 递推法，简称 Levinson 递推法，是一种计算涉及对角线元素为常数的矩阵（称为托普利茨矩阵）方程解的方法。更多信息请参考 en.wikipedia.org/wiki/Levinson_recursion。
Levinson 递推法的实现如下所示。数学推导超出了本书的范围，但递推函数的一般形状你现在应该已经熟悉了。在每次迭代中，我们通过上一轮部分自相关和自协方差的函数来计算部分自相关。
(defn pac* [pacs sigma prev next]
(let [acv (first next)
sum (i/sum (i/mult pacs (reverse prev)))
pac (/ (- acv sum) sigma)]
(cons pac
(lazy-seq
(pac* (->> (i/mult pacs pac)
(reverse)
(i/minus pacs)
(cons pac))
(* (- 1 (i/pow pac 2)) sigma)
(cons acv prev)
(rest next))))))
(defn partial-autocorrelation [series]
(let [acvs (autocovariance series)
acv1 (first acvs)
acv2 (second acvs)
pac (/ acv2 acv1)]
(concat [1.0 pac]
(pac* (vector pac)
(- acv1 (* pac acv2))
(vector acv2)
(drop 2 acvs)))))
如前所述，这个函数将创建一个无限的懒序列部分自相关，所以我们只能从中取出我们实际需要的数字。
绘制部分自相关
现在我们已经实现了一个计算时间序列部分自相关的函数，接下来就来绘制它们。我们将使用与之前相同的 ARMA 系数，以便进行差异比较。
(defn ex-9-24 []
(let [ys (s/sample-normal 10 :sd 1.0)
es (s/sample-normal 10 :sd 0.2)
ps [0 0 0 0.3 0.5]
qs [0.2 0.8]]
(->> (arma ys es ps qs 0.2)
(take 500)
(partial-autocorrelation)
(take 15)
(bar-plot))))
这应该生成一个类似于以下的条形图：
幸运的是，这与我们之前绘制的 ACF 图有所不同。在滞后 1 和滞后 2 处有很高的部分自相关。这表明正在进行一个*MA(2)过程。然后，滞后 5 之前部分自相关较低，这表明也存在一个AR(5)*模型。
使用 ACF 和 PACF 确定 ARMA 模型阶数
ACF 和 PACF 图之间的差异有助于选择最合适的时间序列模型。下表描述了理想化 AR 和 MA 序列的 ACF 和 PACF 图的特征。
| 模型 | ACF | PACF |
|---|---|---|
| AR(p) | 渐渐衰减 | 在p滞后之后截断 |
| MA(q) | 在* q *滞后之后截断 | 渐渐衰减 |
| ARMA(p,q) | 渐渐衰减 | 渐渐衰减 |
然而，我们常常面对的数据并不符合这些理想情况。给定一个实际的时间序列，尤其是没有足够多数据点的序列，通常不容易明确哪个模型最为合适。最好的做法通常是选择最简单的模型（即最低阶的模型），它能有效描述你的数据。
上面的插图展示了理想化*AR(1)序列的 ACF 和 PACF 示例图。接下来是理想化MA(1)*序列的 ACF 和 PACF 示例图。
图中的虚线表示显著性阈值。一般来说，我们无法创建一个完美捕捉时间序列中所有自相关的模型，显著性阈值帮助我们优先考虑最重要的部分。使用 5%显著性水平的显著性阈值的简单公式是：
在这里，n是时间序列中的数据点数量。如果 ACF 和 PACF 中的所有点都接近零，则数据基本上是随机的。
航空数据的 ACF 和 PACF
让我们回到之前开始考虑的航空数据，并绘制数据的前 25 个滞后的 ACF。
(defn ex-9-25 []
(->> (airline-passengers)
(difference)
(autocorrelation)
(take 25)
(bar-plot)))
这段代码生成了以下图表：
可以看到数据中有规律的峰值和谷值。第一个峰值出现在滞后 12；第二个出现在滞后 24。由于数据是按月记录的，这些峰值对应着年度季节性循环。因为我们的时间序列有 144 个数据点，显著性阈值大约为 或 0.17。
接下来，让我们看一下航空数据的部分自相关图：
(defn ex-9-26 []
(->> (airline-passengers)
(difference)
(partial-autocorrelation)
(take 25)
(bar-plot)))
这段代码生成了以下图表：
部分自相关图在滞后 12 处也有一个峰值。与自相关图不同，它在滞后 24 处没有峰值，因为周期性自相关已在滞后 12 时被考虑在内。
尽管这看起来表明 AR(12)模型会合适，但这会创建大量的系数需要学习，特别是在数据量相对较少的情况下。由于周期性循环是季节性的，我们应该通过第二次差分来去除它。
通过差分去除季节性
我们已经对数据进行了第一次差分，这意味着我们的模型被称为自回归积分滑动平均（ARIMA）模型。差分的级别由参数d给出，因此完整的模型阶数可以指定为ARIMA(p,d,q)。
我们可以对数据进行第二次差分，以去除数据中的强季节性。接下来我们来做这个：
(defn ex-9-27 []
(->> (airline-passengers)
(difference)
(difference 12)
(autocorrelation)
(take 15)
(bar-plot)))
首先，我们绘制自相关图：
接下来，部分自相关：
(defn ex-9-28 []
(->> (airline-passengers)
(difference)
(difference 12)
(partial-autocorrelation)
(take 15)
(bar-plot)))
这会生成以下图表：
强季节性周期解释了图表中大部分的显著性。在两个图表中，我们剩下了滞后 1 的负自相关，以及在 ACF 中滞后 9 的微弱显著自相关。一个一般的经验法则是，正自相关最好通过在模型中添加AR项来处理，而负自相关通常最好通过添加MA项来处理。
从之前的图表来看，似乎一个合理的模型是*MA(1)模型。这可能是一个足够好的模型，但我们可以借此机会演示如何通过尝试捕捉AR(9)*自相关来拟合大量的参数。
我们将考虑成本函数的替代方案——似然度，它衡量给定模型与数据的拟合程度。模型拟合得越好，似然度越大。因此，我们希望最大化似然度，这个目标也称为最大似然估计。
最大似然估计
在本书的多次场合，我们已经将优化问题表示为需要最小化的成本函数。例如，在第四章，分类中，我们使用 Incanter 最小化逻辑回归成本函数，同时建立一个逻辑回归分类器；在第五章，大数据中，我们在执行批量和随机梯度下降时使用梯度下降法最小化最小二乘成本函数。
优化也可以表示为要最大化的效益，有时用这种方式思考会更自然。最大似然估计的目标是通过最大化似然函数来找到模型的最佳参数。
假设给定模型参数β时，观测值x的概率表示为：
然后，似然度可以表示为：
似然是给定数据下参数的概率的度量。最大似然估计的目标是找到使观察到的数据最有可能的参数值。
计算似然
在计算时间序列的似然之前，我们将通过一个简单的例子来说明这一过程。假设我们投掷硬币 100 次，观察到 56 次正面，h，44 次反面，t。我们不假设硬币是公平的，P(h)=0.5（因此，稍微不平等的总数是由于偶然变动的结果），而是可以问观察到的值是否与 0.5 有显著差异。我们可以通过询问什么值的*P(h)*使得观察到的数据最有可能来做到这一点。
(defn ex-9-29 []
(let [data 56
f (fn [p]
(s/pdf-binomial data :size 100 :prob p))]
(-> (c/function-plot f 0.3 0.8
:x-label "P"
:y-label "Likelihood")
(i/view))))
在上面的代码中，我们使用二项分布来模拟投掷硬币的序列（回想一下在第四章中提到的分类，二项分布用于模拟二元结果发生的次数）。关键点是数据是固定的，我们在绘制给定不同参数下，观察到该数据的变化概率。以下图表展示了似然面：
正如我们可能预期的，二项分布中最可能的参数是p=0.56。这个人工设定的例子本可以通过手工计算更简单地得出，但最大似然估计的原理能够应对更为复杂的模型。
事实上，我们的 ARMA 模型就是一个复杂的模型。计算时间序列参数似然的数学超出了本书的范围。我们将使用 Clojure 库 Succession（github.com/henrygarner/succession）来计算时间序列的似然。
我们通常使用对数似然而不是似然。这仅仅是出于数学上的方便，因为对数似然：
可以重写为：
在这里，k是模型的参数个数。求多个参数的和比求其积更方便进行计算，因此第二个公式通常更受青睐。让我们通过绘制不同参数下的对数似然与简单的*AR(2)*时间序列的关系，来感受一下似然函数在一些测试数据上的表现。
(defn ex-9-30 []
(let [init (s/sample-normal 2)
coefs [0 0.5]
data (take 100 (ar init coefs 0.2))
f (fn [coef]
(log-likelihood [0 coef] 2 0 data))]
(-> (c/function-plot f -1 1
:x-label "Coefficient"
:y-label "Log-Likelihood")
(i/view))))
上面的代码生成了以下图表：
曲线的峰值对应于根据数据得到的参数最佳估计值。注意，前面图表中的峰值略高于 0.5：我们为模型添加的噪声使得最佳估计并非恰好为 0.5\。
估计最大似然
我们的 ARMA 模型的参数数量很大，因此为了确定最大似然，我们将使用一种在高维空间中表现良好的优化方法。该方法叫做 Nelder-Mead 方法，或称为 simplex 方法。在 n 维空间中，简单形体是一个具有 n+1 个顶点的 多面体。
注意
多面体是一个具有平面边的几何对象，可以存在于任意维度中。二维多边形是 2-多面体，三维多面体是 3-多面体。
简单形体优化的优势在于，它不需要在每个点计算梯度来下降（或上升）到更优位置。Nelder-Mead 方法通过推测在每个测试点测量的目标函数行为来进行优化。最差的点将被一个通过剩余点的质心反射创建的新点所替代。如果新点比当前最优点更好，我们将沿此线指数地拉伸简单形体。如果新点和之前的点没有显著改善，我们可能是跨越了一个低谷，因此我们会将简单形体收缩向可能更优的点。
以下图示展示了简单形体（作为三角形表示）如何反射和收缩以找到最优参数的示例。
简单形体总是表示为一个形状，其顶点数量比维度数多一个。对于二维优化，如前面的图所示，简单形体表示为三角形。对于任意 n 维空间，简单形体将表示为一个 n+1 个顶点的多边形。
注意
简单形体方法也被称为 变形虫方法，因为它看起来像是在朝着更优位置爬行。
简单形体优化方法没有在 Incanter 中实现，但可以在 Apache Commons Math 库中使用（commons.apache.org/proper/commons-math/）。要使用它，我们需要将我们的目标函数——对数似然——包装成库能够理解的表示方式。
使用 Apache Commons Math 进行 Nelder-Mead 优化
Apache Commons Math 是一个庞大且复杂的库。我们无法在这里覆盖太多细节。下一个示例仅用于展示如何将 Clojure 代码与该库提供的 Java 接口集成。
注意
有关 Apache Commons Math 库广泛优化功能的概述，可以参考 commons.apache.org/proper/commons-math/userguide/optimization.html。
Apache Commons Math 库期望我们提供一个要优化的 ObjectiveFunction。接下来，我们通过实现一个 MultivariateFunction 来创建它，因为我们的目标函数需要多个参数。我们的响应将是一个单一值：对数似然。
(defn objective-function [f]
(ObjectiveFunction. (reify MultivariateFunction
(value [_ v]
(f (vec v))))))
上述代码将返回一个ObjectiveFunction表示的任意函数f。MultivariateFunction期望接收一个参数向量v，我们将其直接传递给f。
在此基础上，我们使用一些 Java 互操作性调用optimize，在SimplexOptimizer上使用一些合理的默认值。我们对参数的InitialGuess只是一个零数组。NelderMeadSimplex必须使用每个维度的默认步长进行初始化，该步长可以是任何非零值。我们选择每个参数的步长为 0.2。
(defn arma-model [p q ys]
(let [m (+ p q)
f (fn [params]
(sc/log-likelihood params p q ys))
optimal (.optimize (SimplexOptimizer. 1e-10 1e-10)
(into-array
OptimizationData
[(MaxEval. 100000)
(objective-function f)
GoalType/MAXIMIZE
(->> (repeat 0.0)
(take m)
(double-array)
(InitialGuess.))
(->> (repeat 0.1)
(take m)
(double-array)
(NelderMeadSimplex.))]))
point (-> optimal .getPoint vec)
value (-> optimal .getValue)]
{:ar (take p point)
:ma (drop p point)
:ll value}))
(defn ex-9-31 []
(->> (airline-passengers)
(i/log)
(difference)
(difference 12)
(arma-model 9 1)))
我们的模型是一个大模型，拥有许多参数，因此优化过程会花费一段时间才能收敛。如果你运行上面的示例，你最终会看到类似下面的返回参数：
;; {:ar (-0.23769808471685377 -0.012617164166298971 ...),
;; :ma (-0.14754455658280236),
;; :ll 232.97813750669314}
这些是我们模型的最大似然估计。此外，响应中还包含了使用最大似然参数的模型的对数似然值。
使用赤池信息量准则识别更好的模型
在评估多个模型时，可能会出现最佳模型是最大似然估计值最大的那个模型的情况。毕竟，最大似然估计表明该模型是生成观测数据的最佳候选者。然而，最大似然估计没有考虑模型的复杂性，通常来说，简单的模型更受青睐。回想一下章节开始时我们的高阶多项式模型，它虽然有很高的R²，但却没有提供任何预测能力。
赤池信息量准则（AIC）是一种用于比较模型的方法，它奖励拟合优度（通过似然函数评估），但同时也包括了与参数数量相关的惩罚项。这种惩罚项能有效避免过拟合，因为增加模型的参数数量几乎总是会提高拟合优度。
AIC 可以通过以下公式计算：
在这里，k是模型的参数数量，L是似然函数。我们可以在 Clojure 中通过以下方式使用参数计数p和q来计算 AIC。
(defn aic [coefs p q ys]
(- (* 2 (+ p q 1))
(* 2 (log-likelihood coefs p q ys))))
如果我们要生成多个模型并选择最佳的，我们会选择具有最低 AIC 的那个模型。
时间序列预测
在参数估计已被定义后，我们终于可以使用我们的模型进行预测了。实际上，我们已经编写了大部分需要的代码：我们有一个arma函数，能够根据一些种子数据和模型参数p和q生成自回归滑动平均序列。种子数据将是我们从航空公司数据中测量到的y值，p和q的值是我们使用 Nelder-Mead 方法计算出来的参数。
让我们将这些数字代入我们的 ARMA 模型，并生成y的预测序列：
(defn ex-9-32 []
(let [data (i/log (airline-passengers))
diff-1 (difference 1 data)
diff-12 (difference 12 diff-1)
forecast (->> (arma (take 9 (reverse diff-12))
[]
(:ar params)
(:ma params) 0)
(take 100)
(undifference 12 diff-1)
(undifference 1 data))]
(->> (concat data forecast)
(i/exp)
(timeseries-plot))))
上述代码生成了以下图表：
到时间切片 144 为止的线条是原始序列。从此点之后的线条是我们的预测序列。预测看起来与我们预期的差不多：指数增长趋势持续下去，季节性波动的峰值和谷值模式也在继续。
实际上，预测几乎太规则了。与第 1 到 144 点的系列不同，我们的预测没有噪声。让我们加入一些噪声，使我们的预测更具现实性。为了确定噪声的合理性，我们可以查看过去预测中的误差。为了避免错误的积累，我们应该每次预测一个时间步长，并观察预测值与实际值之间的差异。
让我们用 0.02 的 sigma 值运行我们的 ARMA 函数：
(defn ex-9-33 []
(let [data (i/log (airline-passengers))
diff-1 (difference 1 data)
diff-12 (difference 12 diff-1)
forecast (->> (arma (take 9 (reverse diff-12))
[]
(:ar params)
(:ma params) 0.02)
(take 10)
(undifference 12 diff-1)
(undifference 1 data))]
(->> (concat data forecast)
(i/exp)
(timeseries-plot))))
上面的代码可能会生成如下图表：
现在我们可以感知到预测的波动性了。通过多次运行模拟，我们可以了解不同可能结果的多样性。如果我们能够确定预测的置信区间，包括噪声的上下期望值，这将非常有用。
蒙特卡洛模拟预测
尽管确实存在用于计算时间序列预期未来值的方法，并且可以得到置信区间，我们将在本节通过模拟来获得这些值。通过研究多个预测之间的变异性，我们可以为模型预测得出置信区间。
例如，如果我们运行大量的模拟，我们可以计算基于 95%时间内值范围的未来预测的 95%置信区间。这正是蒙特卡洛模拟的核心，它是一个常用于解决分析上难以处理问题的统计工具。
注意
蒙特卡洛方法是在曼哈顿计划期间开发并系统地使用的，这是美国在二战期间开发核武器的努力。约翰·冯·诺依曼和斯坦尼斯瓦夫·乌拉姆建议将其作为研究中子穿越辐射屏蔽的性质的一种手段，并以摩纳哥蒙特卡洛赌场命名该方法。
我们已经为时间序列预测的蒙特卡洛模拟奠定了所有基础。我们只需要多次运行模拟并收集结果。在以下代码中，我们运行 1,000 次模拟，并收集每个未来时间切片下所有预测的均值和标准差。通过创建两个新序列（一个上界，通过将标准差乘以 1.96 并加上标准差，另一个下界，通过将标准差乘以 1.96 并减去标准差），我们能够可视化系列未来值的 95%置信区间。
(defn ex-9-34 []
(let [data (difference (i/log (airline-passengers)))
init (take 12 (reverse data))
forecasts (for [n (range 1000)]
(take 20
(arma init [0.0028 0.0028]
(:ar params1)
(:ma params1)
0.0449)))
forecast-mean (map s/mean (i/trans forecasts))
forecast-sd (-> (map s/sd (i/trans forecasts))
(i/div 2)
(i/mult 1.96))
upper (->> (map + forecast-mean forecast-sd)
(concat data)
(undifference 0)
(i/exp))
lower (->> (map - forecast-mean forecast-sd)
(concat data)
(undifference 0)
(i/exp))
n (count upper)]
(-> (c/xy-plot (range n) upper
:x-label "Time"
:y-label "Value"
:series-label "Upper Bound"
:legend true)
(c/add-lines (range n) lower
:series-label "Lower Bound")
(i/view))))
这生成了以下图表：
上限和下限为我们对未来时间序列预测提供了置信区间。
总结
在本章中，我们考虑了分析离散时间序列的任务：在固定时间间隔内进行的顺序观察。我们看到，通过将序列分解为一组组件：趋势成分、季节性成分和周期性成分，可以使建模这种序列的挑战变得更加容易。
我们看到了 ARMA 模型如何将一个序列进一步分解为自回归（AR）和滑动平均（MA）组件，每个组件都以某种方式由序列的过去值决定。这种对序列的理解本质上是递归的，我们也看到了 Clojure 自然支持定义递归函数和懒序列的能力，如何有助于算法生成此类序列。通过将序列的每个值确定为前一个值的函数，我们实现了一个递归的 ARMA 生成器，它能够模拟一个已测量的序列并进行未来的预测。
我们还学习了期望最大化：一种将优化问题的解重新表述为根据数据生成最大可能性的解的方法。我们还看到了如何使用 Apache Commons Math 库，通过 Nelder-Mead 方法估计最大似然参数。最后，我们了解了如何通过将序列向前推进来进行预测，以及如何使用蒙特卡洛模拟来估计序列的未来误差。
在最后一章，我们将注意力从数据分析转向数据可视化。从某种程度上说，数据科学家的最重要挑战是沟通，我们将看到 Clojure 如何帮助我们以最有效的方式展示数据。
第十章. 可视化
| “数字有一个重要的故事要告诉我们。它们依赖于你，给它们一个清晰且有说服力的声音。” | ||
|---|---|---|
| --Stephen Few |
本书中的每一章都以某种方式使用了可视化，主要使用 Incanter。Incanter 是一个有效的工具，可以在工作中制作各种各样的图表，这些图表通常是我们在试图理解数据集时首先会使用的。这个初步阶段通常被称为探索性数据分析，在这个阶段，我们感兴趣的是总结数据的统计信息，例如数值数据的分布、类别数据的计数，以及数据中属性之间的相关性。
在找到一种有意义的方式来解读数据后，我们通常会希望将其传达给他人。最重要的沟通工具之一就是可视化，我们可能需要向没有强大分析背景的人传达细微或复杂的概念。在本章中，我们将使用 Quil 库——这是为视觉艺术家开发的软件的延伸——来制作吸引人的图形，帮助使数据生动起来。可视化和沟通设计是广泛而丰富的领域，我们在这里无法详细覆盖。相反，本章将提供两个案例研究，展示如何将 Clojure 的数据抽象和 Quil 的绘图 API 结合起来，产生良好的效果。
本章的开始，我们将回到第一章，统计学中使用的数据。我们将通过演示如何从俄罗斯选举数据构建一个简单的二维直方图来介绍 Quil。在掌握了 Quil 绘图的基本知识后，我们将展示如何通过一些基本的绘图指令，结合起来，呈现美国财富分配的引人注目的图示。
下载代码和数据
在本章中，我们将回到本书第一章中使用的数据：2011 年俄罗斯选举的数据。在第一章，统计学中，我们使用了带透明度的散点图来可视化选民投票率与胜选者选票百分比之间的关系。在本章中，我们将编写代码，将数据渲染为二维直方图。
我们还将使用有关美国财富分配的数据。这个数据非常小，以至于我们不需要下载任何东西：我们将直接在源代码中输入这些数据。
注意
本章的源代码可以在github.com/clojuredatascience/ch10-visualization找到。
本章的示例代码包含一个脚本，用于下载我们在第一章，统计学中使用的选举数据。一旦下载了源代码，你可以通过在项目根目录内运行以下命令来执行脚本：
script/download-data.sh
如果你之前下载了第一章的统计学数据，你可以将数据文件直接移动到本章的数据目录中（如果你愿意的话）。
探索性数据可视化
在任何数据科学项目的初期，通常会有一个迭代数据探索的阶段，你可以在这个阶段获得对数据的洞察。在本书中，Incanter 一直是我们的主要可视化工具。尽管它包含了大量的图表，但也会有一些场合，它不包含你想要表示的数据的理想图表。
注意
其他 Clojure 库正在提供探索性数据可视化功能。例如，查看clojurewerkz/envision github.com/clojurewerkz/envision和 Karsten Schmidt 的thi-ng/geom，地址为github.com/thi-ng/geom/tree/master/geom-viz。
例如，在第一章的统计学中，我们使用了带有透明度的散点图来可视化选民投票率与赢家得票比例之间的关系。这不是理想的图表，因为我们主要关注的是某个特定区域内点的密度。透明度有助于揭示数据的结构，但它不是一种明确的表示。一些点仍然太微弱，无法看清，或者数量太多，以至于它们看起来像一个点：
我们本可以通过二维直方图解决这些问题。这种类型的图表使用颜色来传达二维空间中高低密度的区域。图表被划分为一个网格，网格中的每个单元格表示两个维度的一个范围。点越多地落入网格的单元格中，该范围内的密度就越大。
表示二维直方图
直方图只是将连续分布表示为一系列箱子的方式。直方图在第一章的统计学中已经介绍，当时我们编写了一个分箱函数，将连续数据分为离散的箱子：
(defn bin [n-bins xs]
(let [min-x (apply min xs)
range-x (- (apply max xs) min-x)
max-bin (dec n-bins)
bin-fn (fn [x]
(-> (- x min-x)
(/ range-x)
(* n-bins)
(int)
(min max-bin)))]
(map bin-fn xs)))
这段代码将对连续的xs范围进行分箱，并根据n-bins参数将其分成不同的组。例如，将 0 到 19 的范围分为 5 个箱子，结果如下：
(defn ex-1-1 []
(bin 5 (range 20)))
;;(0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4)
bin函数返回每个数据点的箱子索引，而不是计数，因此我们使用 Clojure 的frequencies函数来确定落入该箱子的点的数量：
(defn ex-1-2 []
(frequencies (bin 5 (range 20))))
;;{0 4, 1 4, 2 4, 3 4, 4 4}
这是一个合理的一维直方图表示：它将计数的箱子映射出来。要表示二维直方图，我们只需要对xs和ys执行相同的计算。我们将向量函数映射到箱子的索引上，以便将每个点转换为[x-bin y-bin]的表示：
(defn histogram-2d [xsys n-bins]
(-> (map vector
(bin n-bins xs)
(bin n-bins ys))
(frequencies)))
这个函数返回一个以两个值的向量为键的映射。frequencies 函数现在将计数所有在 x 和 y 轴上都共享同一个区间的点：
(defn ex-10-3 []
(histogram-2d (range 20)
(reverse (range 20)) 5))
;;{[0 4] 4, [1 3] 4, [2 2] 4, [3 1] 4, [4 0] 4}
我们希望在直方图中绘制实际数据，因此让我们从 第一章加载俄罗斯数据，统计学。如果你已经将数据下载到示例代码的 data 目录中，可以运行以下代码：
(defn ex-10-4 []
(let [data (load-data :ru-victors)]
(histogram-2d (i/$ :turnout data)
(i/$ :victors-share data) 5)))
;; {[4 3] 6782, [2 2] 14680, [0 0] 3, [1 0] 61, [2 3] 2593,
;; [3 3] 8171, [1 1] 2689, [3 4] 1188, [4 2] 3084, [3 0] 64,
;; [4 1] 1131, [1 4] 13, [1 3] 105, [0 3] 6, [2 4] 193, [0 2] 10,
;; [2 0] 496, [0 4] 1, [3 1] 3890, [2 1] 24302, [4 4] 10771,
;; [1 2] 1170, [3 2] 13384, [0 1] 4, [4 0] 264}
我们可以看到直方图区间中的值范围巨大：从区间 [0 4] 中的 1 到区间 [2 1] 中的 24,302。这些计数值将是我们在直方图上绘制的密度值。
使用 Quil 进行可视化
Quil (github.com/quil/quil) 是一个 Clojure 库，提供了大量的灵活性来生成自定义的可视化效果。它封装了 Processing (processing.org/)，这是一个 Java 框架，已被视觉艺术家和设计师积极开发多年，旨在促进“视觉艺术中的软件素养和技术中的视觉素养”。
使用 Quil 进行的任何可视化都涉及创建一个 sketch。sketch 是处理程序的术语，指的是运行一个由绘图指令组成的程序。大多数 API 函数都可以从 quil.core 命名空间中调用。我们将其作为 q 包含在代码中。调用 q/sketch 并不传递任何参数时，将会弹出一个空窗口（尽管它可能会被其他窗口遮挡）。
绘制到 sketch 窗口
默认的窗口大小是 500px x 300px。我们希望我们的二维直方图是正方形的，因此将窗口的宽高都设置为 250px：
(q/sketch :size [250 250])
由于我们每个轴都有 5 个区间，因此每个区间将由一个 50px 宽和 50px 高的正方形表示。
Quil 提供了标准的二维形状原语用于绘制：点、线、弧、三角形、四边形、矩形和椭圆。要绘制一个矩形，我们调用 q/rect 函数，并指定 x 和 y 坐标，以及宽度和高度。
让我们在原点绘制一个宽度为 50px 的正方形。有几种方式可以向 Quil 提供绘图指令，但在本章中，我们将传递一个被称为 setup 的函数。这是一个没有参数的函数，我们将其传递给 sketch。我们的零参数函数仅仅是调用 rect，并传入位置 [0, 0] 和宽高为 50：
(defn ex-10-5 []
(let [setup #(q/rect 0 0 50 50)]
(q/sketch :setup setup
:size [250 250])))
该代码生成了以下图像：
根据你对计算机图形学的熟悉程度，矩形的位置可能与预期不同。
注意
通过传递半径作为第五个参数，也可以绘制带有圆角的矩形。通过传递第五到第八个参数的值，可以为每个角使用不同的半径。
在继续之前，我们需要理解 Quil 的坐标系统。
Quil 的坐标系统
Quil 使用的坐标系统与 Processing 和大多数其他计算机图形程序相同。如果你不熟悉绘图，可能会觉得起点位于显示屏左上角是直觉上不对的。y 轴向下延伸，x 轴向右延伸。
很明显，这不是大多数图表中 y 轴的方向，这意味着在绘图时，y 坐标通常需要被翻转。
一种常见的方法是从草图底部测量的 y 值中减去所需的值，这样的转换使得 y 为零时对应草图的底部。较大的 y 值则对应草图中更高的位置。
绘制网格
让我们通过一个简单的网格来实践这个。以下函数接受一个箱子的数量 n-bins 和一个 size 参数，size 表示为 [width height] 的向量：
defn draw-grid [{:keys [n-bins size]}]
(let [[width height] size
x-scale (/ width n-bins)
y-scale (/ height n-bins)
setup (fn []
(doseq [x (range n-bins)
y (range n-bins)
:let [x-pos (* x x-scale)
y-pos (- height
(* (inc y) y-scale))]]
(q/rect x-pos y-pos x-scale y-scale)))]
(q/sketch :setup setup :size size)))
从中，我们可以计算出 x-scale 和 y-scale，这是一个因子，使我们能够将箱子索引转换为 x 和 y 维度中的像素偏移。这些值被我们的 setup 函数使用，setup 函数遍历 x 和 y 箱子，为每个箱子放置一个矩形。
注意
请注意，我们在 doseq 内部执行了循环。我们的绘图指令作为副作用执行。如果我们不这样做，Clojure 的惰性求值将导致没有任何绘制操作。
之前的代码生成了以下图形：
在定义了前面的函数后，我们几乎已经创建了一个直方图。我们只需要为网格中的每个方格着色，颜色代表每个箱子的适当值。为了实现这一点，我们需要两个额外的函数：一个从数据中获取与箱子对应的值，另一个将这些值解释为颜色。
指定填充颜色
在 Quil 中填充颜色是通过 q/fill 函数实现的。我们指定的任何填充将在我们指定新的填充之前一直使用。
注意
Quil 中的许多函数会影响当前的绘图上下文，并且是 有状态的。例如，当我们指定一个填充值时，它将在后续的所有绘图指令中使用，直到填充值被更改。其他例子包括填充、笔触、缩放和字体。
以下代码是我们 draw-grid 函数的一个改编版本。draw-filled-grid 的新增部分是 fill-fn：用于给网格中的点着色的方式。fill-fn 函数应该是一个接受 x 和 y 箱子索引作为参数的函数，它应返回 Quil 可以用作填充的表示：
(defn draw-filled-grid [{:keys [n-bins size fill-fn]}]
(let [[width height] size
x-scale (/ width n-bins)
y-scale (/ height n-bins)
setup (fn []
(doseq [x (range n-bins)
y (range n-bins)
:let [x-pos (* x x-scale)
y-pos (- height
(* (inc y) y-scale))]]
(q/fill (fill-fn x y))
(q/rect x-pos y-pos x-scale y-scale)))]
(q/sketch :setup setup :size size)))
Quil 的填充函数接受多个参数：
-
一个参数：RGB 值（可以是一个数字或
q/color表示法） -
两个参数：与单个参数的情况一样，除了加上一个 alpha 透明度值
-
三个参数：颜色的红、绿、蓝分量，作为 0 到 255 之间的数字（包括 0 和 255）
-
四个参数：红色、绿色、蓝色和透明度的数值
我们很快就会看到如何使用 color 表示法，但现在，我们将使用一个简单的数字表示法来表示颜色：介于 0 和 255 之间的数字。当红色、绿色和蓝色的数值相同时（或者当调用 fill 函数时传递一个或两个参数），我们得到一种灰色。0 对应黑色，255 对应白色。
如果我们将每个柱状图中值的频率除以最大值，我们将得到一个介于 0 和 1.0 之间的数字。将其乘以 255 将得到一个 Quil 会转换为灰色的值。我们在以下 fill-fn 实现中做了这个，并将其传递给我们之前定义的 draw-filled-grid 函数：
(defn ex-10-6 []
(let [data (load-data :ru-victors)
n-bins 5
hist (histogram-2d (i/$ :turnout data)
(i/$ :victors-share data)
n-bins)
max-val (apply max (vals hist))
fill-fn (fn [x y]
(-> (get hist [x y] 0)
(/ max-val)
(* 255)))]
(draw-filled-grid {:n-bins n-bins
:size [250 250]
:fill-fn fill-fn})))
前面的代码生成了以下图形：
图表完成了我们想要的功能，但它是我们数据的非常粗略的表示。让我们增加箱子的数量来提高柱状图的分辨率：
(defn ex-10-7 []
(let [data (load-data :ru-victors)
n-bins 25
hist (histogram-2d (i/$ :turnout data)
(i/$ :victors-share data)
n-bins)
max-val (apply max (vals hist))
fill-fn (fn [x y]
(-> (get hist [x y] 0)
(/ max-val)
(* 255)))]
(draw-filled-grid {:n-bins n-bins
:size [250 250]
:fill-fn fill-fn})))
这段代码生成了以下图形：
在 x 和 y 轴上各有 25 个矩形，我们就能获得数据结构的更细粒度的图像。然而，副作用是由于大多数单元格的色调较暗，柱状图的细节变得难以辨认。部分问题在于右上角的值如此之高，以至于即便是中央区域（之前最亮的部分）现在也不过是一个灰色的模糊斑点。
这个问题有两个解决方案：
-
通过绘制 z-分数而不是实际值来减轻离群值的影响
-
通过使用更广泛的颜色范围来多样化视觉线索
我们将在下一节中学习如何将值转换为完整的颜色谱，但首先，让我们将柱状图的值转换为 z-分数。绘制 z-分数是以分布为意识的方式着色图表，这将大大减小右上角极端离群值的影响。使用 z-分数，我们将绘制每个单元格与均值之间的标准差数量。
为了实现这一点，我们需要知道两件事：柱状图中频率的均值和标准差：
(defn ex-10-8 []
(let [data (load-data :ru-victors)
n-bins 25
hist (histogram-2d (i/$ :turnout data)
(i/$ :victors-share data)
n-bins)
mean (s/mean (vals hist))
sd (s/sd (vals hist))
fill-fn (fn [x y]
(-> (get hist [x y] 0)
(- mean)
(/ sd)
(q/map-range -1 3 0 255)))]
(draw-filled-grid {:n-bins n-bins
:size [250 250]
:fill-fn fill-fn})))
前面的代码从柱状图中的每个值中减去均值，然后除以均值。这样会得到一个均值为零的值。1 将代表离均值一个标准差，2 将代表离均值两个标准差，依此类推。
Quil 暴露了一个有用的map-range函数，它可以将一个值范围映射到另一个值范围。例如，我们可以将所需的标准差范围（在前面的例子中是-1 到 3）映射到 0 到 255 的范围。这样，分布的四个标准差将对应于从黑到白的完整灰度范围。任何超过此范围的数据将被裁剪。
结果是一个更引人注目的灰度数据表示方式。使用z分数使得直方图的主体部分呈现出更多细节，我们现在能够察觉到尾部的更多变化。
然而，尽管如此，直方图仍然没有达到理想的清晰度，因为区分不同的灰度色调可能会很有挑战性。在某些不相邻的单元格之间，可能很难判断它们是否共享相同的值。
我们可以通过利用颜色来表示每个单元格，从而扩展可用的范围。这使得直方图更像热力图：“较冷”的颜色，如蓝色和绿色，代表低值，而“较热”的颜色，如橙色和红色，则代表热力图中最密集的区域。
颜色与填充
为了创建我们二维直方图的热力图版本，我们需要将* z *分数与颜色值进行映射。与显示离散的颜色调色板（例如 5 种颜色）不同，我们的热力图应该有一个平滑的调色板，包含光谱中的所有颜色。
注意
对于那些阅读印刷版书籍或黑白版本的读者，你可以从 Packt 出版社的网站下载彩色图像：www.packtpub.com/sites/default/files/downloads/Clojure_for_Data_Science_ColorImages.pdf。
这正是 Quil 函数q/lerp-color所做的。给定两种颜色和介于 0 与 1 之间的比例，lerp-color将返回两者之间插值的新颜色。0 的比例返回第一个颜色，1 的比例返回第二个颜色，而 0.5 的比例则返回两者之间的颜色：
(defn z-score->heat [z-score]
(let [colors [(q/color 0 0 255) ;; Blue
(q/color 0 255 255) ;; Turquoise
(q/color 0 255 0) ;; Green
(q/color 255 255 0) ;; Yellow
(q/color 255 0 0)] ;; Red
offset (-> (q/map-range z-score -1 3 0 3.999)
(max 0)
(min 3.999))]
(q/lerp-color (nth colors offset)
(nth colors (inc offset))
(rem offset 1))))
这段代码使用了按光谱顺序排列的颜色数组。我们用q/map-range来确定我们将要插值的两个颜色，并使用q/lerp-color与范围的小数部分进行插值。
我们已经实现了一个draw-filled-grid函数，它接受fill-fn来决定使用哪种颜色填充网格。现在，让我们将z-score->heat函数传递给它：
(defn ex-10-9 []
(let [data (load-data :ru-victors)
n-bins 25
hist (histogram-2d (i/$ :turnout data)
(i/$ :victors-share data)
n-bins)
mean (s/mean (vals hist))
sd (s/sd (vals hist))
fill-fn (fn [x y]
(-> (get hist [x y] 0)
(- mean)
(/ sd)
(z-score->heat)))]
(draw-filled-grid {:n-bins n-bins
:size [250 250]
:fill-fn fill-fn})))
这段代码生成了以下图形：
热力图展示了数据的更多内部结构。尤其是，虽然数据的强烈对角形状依然显现，但我们现在可以看到更多的变化。先前难以确定的细节（无论是因为区域过于密集还是过于稀疏）变得更加明显。
输出图像文件
现在我们对直方图感到满意，接下来我们希望输出高质量的版本。通过在设置函数中加入q/save命令，并传入文件名，Quil 将同时将图表输出到文件和屏幕上。图像创建的格式将取决于文件名的后缀：.tif为 TIFF 文件，.jpg为 JPEG 文件，.png为 PNG 文件，.tga为 TARGA 文件：
(defn draw-filled-grid [{:keys [n-bins size fill-fn]}]
(let [[width height] size
x-scale (/ width n-bins)
y-scale (/ height n-bins)
setup (fn []
(doseq [x (range n-bins)
y (range n-bins)
:let [x-pos (* x x-scale)
y-pos (- height
(* (inc y) y-scale))]]
(q/fill (fill-fn x y))
(q/rect x-pos y-pos x-scale y-scale))
(q/save "heatmap.png"))]
(q/sketch :setup setup :size size)))
我们还可以将输出结果保存为 PDF 格式，正如接下来的可视化展示。
用于沟通的可视化
作为数据科学家，在工作中我们可能需要与各种各样的人沟通。我们身边的同事和经理可能能够阅读并解读我们的 Incanter 图表，但这不会给 CEO 留下深刻印象。我们可能还需要与公众沟通。
无论是哪种情况，我们都应该专注于制作简单且有力的可视化图表，同时又不牺牲数据的完整性。缺乏统计学训练并不妨碍我们理解那些微妙而复杂的论点，我们应该尊重观众的智慧。作为数据科学家，我们面临的挑战是找到一种有效地传达信息的表现形式。
在本章的剩余部分，我们将制作一个可视化图表，旨在以简洁而真实的方式传达更复杂的数据集。
注意
我们将创建的可视化图表是来自于美国财富不平等视频中的一个图表版本，视频可以在www.youtube.com/watch?v=QPKKQnijnsM观看。该视频由匿名电影制作人 Politizane 制作，强有力的视觉冲击使其在 YouTube 上获得了超过 1600 万次点击。
像这样的图形展示常常有一个问题，我们的数据将来自多个不同的来源。
可视化财富分布
我们将使用的第一个数据集来自加利福尼亚大学圣克鲁兹分校的心理学与社会学研究教授 G. William Domhoff 的一篇文章。我们接下来引用的数字来自一篇名为财富、收入与权力的文章，地址为www2.ucsc.edu/whorulesamerica/power/wealth.html。
虽然整篇文章都值得阅读，但其中一张特别引人注目的图表是一张饼图，展示了 2010 年美国人财务净资产的分布：
这张饼图之所以引人注目，有几个原因。首先，40%以上的财富掌握在如此小的群体手中，这一概念令人难以理解。其次，饼图的每一块不仅代表着截然不同的财富数量，还代表着截然不同的人群数量：从占人口 1%的群体到占人口 80%的群体。饼图本就因阅读困难而著名，因此这张图表的挑战性更大。
注意
饼图通常不是表示数据的好方式，即使总数在概念上代表整体的一部分。作者兼程序员 Steve Fenton 已经记录了许多原因，并提供了适当的替代方法，详见www.stevefenton.co.uk/2009/04/pie-charts-are-bad/。
让我们看看如何重新解读这些数据，让它更易于理解。首先，提取出以下表格中所列的数字：
| 百分位 | 总金融财富，2010 |
|---|---|
| 0-79 | 5% |
| 80-89 | 11% |
| 90-95 | 13% |
| 96-99 | 30% |
| 100 | 42% |
相较于饼图，改进的小方法是将相同的数据表示为条形图。虽然人们通常很难成功地解读饼图中各个扇形的相对大小，条形图则不会遇到此类问题。下一个示例将简单地将之前的数据转换为条形图：
(defn ex-10-10 []
(let [categories ["0-79" "80-89" "90-95" "96-99" "100"]
percentage [5 11 13 30 42 ]]
(-> (c/bar-chart categories percentage
:x-label "Category"
:y-label "% Financial Wealth")
(i/view))))
这将返回以下图表：
这是对饼图的一种改进，因为它更容易比较各类别的相对大小。然而，仍然存在一个显著的问题：每个类别所代表的人数差异巨大。左侧的条形代表了 80%的人口，而右侧的条形则代表了 1%的人口。
如果我们想让这些数据更加易于理解，我们可以将总量划分为 100 个相等的单位，每个单位代表人口的一个百分位。每个条形的宽度可以根据其所代表的百分位数进行调整，同时保持其面积不变。由于每个百分位单位代表相同数量的人口，因此生成的图表可以让我们更容易地进行不同群体间的比较。
我们可以通过返回一个包含 100 个元素的序列来实现这一目标，每个元素代表一个人口百分位。序列中每个元素的值将是该百分位所占的财富比例。我们已经知道，前 1%的人口拥有 42%的财富，但其他群体将根据它们所代表的百分位数调整其值：
(def wealth-distribution
(concat (repeat 80 (/ 5 80))
(repeat 10 (/ 11 10))
(repeat 5 (/ 13 5))
(repeat 4 (/ 30 4))
(repeat 1 (/ 42 1))))
(defn ex-10-11 []
(let [categories (range (count wealth-distribution))]
(-> (c/bar-chart categories wealth-distribution
:x-label "Percentile"
:y-label "% Financial Wealth")
(i/view))))
这个示例生成了以下的条形图：
通过应用一个简单的转换，我们能够更好地理解真正的分布情况。每个条形现在代表着相同比例的人口，条形的面积代表该百分位所拥有的财富比例。
使用 Quil 使数据更生动
前一部分的变换结果是一个几乎过于直观显示极端值差异的图表：除了最大条形外，几乎无法解读任何内容。一种解决方案是使用对数刻度或对数-对数刻度来显示数字，正如我们在本书的其他地方所做的那样。如果这个图表的观众是统计学上有素养的人，那么这样做可能是最合适的，但我们假设我们的可视化面向的目标观众是普通大众。
之前展示的图表的问题是最右侧的条形图过大，压倒了其他所有条形。80%的面积仅由几个像素表示。在下一部分中，我们将利用 Quil 制作一个更好地利用空间，同时又能保持图表完整性的可视化。
绘制不同宽度的条形图
在接下来的几个部分中，我们将分阶段构建一个可视化。由于我们将绘制一个 Quil 草图，我们首先会定义一些常数，以便相对于草图的尺寸生成绘图指令。为了简洁起见，下一段代码省略了一些常数：
(def plot-x 56)
(def plot-y 60)
(def plot-width 757)
(def plot-height 400)
(def bar-width 7)
在这些完成后，我们可以开始以更易于理解的方式呈现条形图。以下代码采用财富分布并将除了最后一条条形以外的所有条形绘制为一系列矩形。y轴的刻度被计算出来，以便我们绘制的最大条形图填满整个绘图的高度：
(defn draw-bars []
(let [pc99 (vec (butlast wealth-distribution))
pc1 (last wealth-distribution)
y-max (apply max pc99)
y-scale (fn [x] (* (/ x y-max) plot-height))
offset (fn [i] (* (quot i 10) 7))]
(dotimes [i 99] ;; Draw the 99%
(let [bar-height (y-scale (nth pc99 i))]
(q/rect (+ plot-x (* i bar-width) (offset i))
(+ plot-y (- plot-height bar-height))
bar-width bar-height)))
(let [n-bars 5 ;; Draw the 1%
bar-height (y-scale (/ pc1 n-bars))]
(q/rect (+ plot-x (* 100 bar-width) (offset 100))
(+ plot-y (- plot-height bar-height))
(* bar-width n-bars) bar-height))))
到目前为止，我们绘制的条形图代表了 99%的比例。最后一条条形将代表人口的最后 1%。为了使其适应我们设计的垂直刻度，且不至于从草图的顶部消失，我们将相应地加宽该条形，同时保持其面积。因此，这条条形比其他条形短 5 倍——但也宽 5 倍：
(defn ex-10-12 []
(let [size [960 540]]
(q/sketch :size size
:setup draw-bars)))
示例输出如下图：
到目前为止，我们已经能够更清晰地看到最大条形之间的关系，但还不明显能看出这是一个图表。在接下来的部分中，我们将添加文本，以标明图表的主题和坐标轴的范围。
添加标题和坐标轴标签
专用可视化工具（如 Incanter）的一大便利之处是，坐标轴可以自动生成。Quil 在这里没有提供帮助，但由于条形宽度已知，因此我们并不难实现。在下面的代码中，我们将使用text、text-align、text-size等函数将文本写入我们的可视化中：
(defn group-offset [i]
(* (quot i 10) 7))
(defn draw-axis-labels []
(q/fill 0)
(q/text-align :left)
(q/text-size 12)
(doseq [pc (range 0 (inc 100) 10)
:let [offset (group-offset pc)
x (* pc bar-width)]]
(q/text (str pc "%") (+ plot-x x offset) label-y))
(q/text "\"The 1%\"" pc1-label-x pc1-label-y))
使用非专业图表库所失去的，我们在灵活性上得到了补偿。接下来，我们将编写一个函数，在文本上制作凸印风格的压印效果：
(defn emboss-text [text x y]
(q/fill 255)
(q/text text x y)
(q/fill 100)
(q/text text x (- y 2)))
(defn draw-title []
(q/text-size 35)
(q/text-leading 35)
(q/text-align :center :top)
(emboss-text "ACTUAL DISTRIBUTION\nOF WEALTH IN THE US"
title-x title-y))
我们使用emboss-text函数在图表的中心绘制一个大标题。注意，我们还指定了文本的对齐方式，可以选择从文本的顶部、底部、中心、左侧或右侧来测量位置：
(defn ex-10-13 []
(let [size [960 540]]
(q/sketch :size size
:setup #((draw-bars)
(draw-axis-labels)
(draw-title)))))
之前的例子生成了以下图形：
这个图表结合了条形图高度、面积和自定义文本可视化，使用标准图表应用程序很难实现。使用 Quil，我们有一个工具箱，可以轻松地将图形和数据自由混合。
通过插图提高清晰度
我们的图表已经有了一些进展，但目前看起来非常简洁。为了增加更多的视觉吸引力，可以使用图像。在示例项目的资源目录中，有两个 SVG 图像文件，一个是人形图标，另一个是来自维基百科的美国地图。
注意
维基百科提供了各种各样的 SVG 地图，这些地图都可以在灵活的创作共享许可证下使用。例如，commons.wikimedia.org/wiki/Category:SVG_maps_of_the_United_States 上的美国地图。
我们在本章使用的地图可以在commons.wikimedia.org/wiki/File:Blank_US_Map,_Mainland_with_no_States.svg找到，由 Lokal_Profil 在 CC-BY-SA-2.5 许可证下提供。
在 Quil 中使用 SVG 图像是一个两步过程。首先，我们需要使用q/load-shape将图像加载到内存中。此函数接受一个参数：要加载的 SVG 文件的路径。接下来，我们需要将图像实际绘制到屏幕上。这是通过使用q/shape函数来完成的，该函数需要一个* x ， y *位置以及可选的宽度和高度。如果我们使用的是基于像素的图像（如 JPEG 或 PNG），我们将使用相应的q/load-image和q/image函数：
(defn draw-shapes []
(let [usa (q/load-shape "resources/us-mainland.svg")
person (q/load-shape "resources/person.svg")
colors [(q/color 243 195 73)
(q/color 231 119 46)
(q/color 77 180 180)
(q/color 231 74 69)
(q/color 61 76 83)]]
(.disableStyle usa)
(.disableStyle person)
(q/stroke 0 50)
(q/fill 200)
(q/shape usa 0 0)
(dotimes [n 99]
(let [quintile (quot n 20)
x (-> (* n bar-width)
(+ plot-x)
(+ (group-offset n)))]
(q/fill (nth colors quintile))
(q/shape person x icons-y icon-width icon-height)))
(q/shape person
(+ plot-x (* 100 bar-width) (group-offset 100))
icons-y icon-width icon-height)))
在这段代码中，我们对usa和person形状都调用了.disableStyle。这是因为 SVG 文件可能包含嵌入的样式信息，比如填充颜色、描边颜色或边框宽度，这些都会影响 Quil 绘制形状的方式。为了完全控制我们的表现形式，我们选择禁用所有样式。
此外，请注意，我们只加载一次 person 形状，并通过dotimes绘制多次。我们根据用户所处的quintile设置颜色：
(defn ex-10-14 []
(let [size [960 540]]
(q/sketch :size size
:setup #((draw-shapes)
(draw-bars)
(draw-axis-labels)
(draw-title)))))
结果如下一张图所示：
这张图已经开始看起来像我们可以展示给人们而不感到羞愧的那种。人物图标帮助传达了每个条形代表一个人口百分位的概念。条形图现在看起来还不太吸引人。由于每个条形代表每个人的财富，所以我们把每个条形画成一叠钞票。虽然这可能看起来是一个过于字面化的解释，但实际上会更清晰地传达出 1%的条形实际上是其他所有条形的 5 倍宽。
向条形图中添加文字
到现在为止，我们已经可以将钞票画成一系列矩形了：
(defn banknotes [x y width height]
(q/no-stroke)
(q/fill 80 127 64)
(doseq [y (range (* 3 (quot y 3)) (+ y height) 3)
x (range x (+ x width) 7)]
(q/rect x y 6 2)))
上述代码的唯一复杂之处在于需要将起始的y位置调整为3的偶数倍。这将确保所有的钞票在y轴上的高度如何，最终都能与x轴对齐。这是从上到下绘制条形的副作用，而不是反向绘制。
我们将在以下示例中将之前的函数添加到我们的草图中：
(defn ex-10-15 []
(let [size [960 540]]
(q/sketch :size size
:setup #((draw-shapes)
(draw-banknotes)
(draw-axis-labels)
(draw-title)))))
这将生成以下图表：
现在这已经是一个相对完整的图表，代表了美国实际的财富分布。之前提供的 YouTube 视频链接的一个优点是，它将实际分布与几种其他分布进行了对比：人们预期的财富分布和他们理想的财富分布。
融入额外数据
哈佛商学院教授迈克尔·诺顿和行为经济学家丹·阿里利对 5000 多名美国人进行了研究，评估他们对财富分配的看法。当他们展示了各种财富分配的例子，并让他们判断哪个例子来自美国时，大多数人选择了一个比实际情况更加平衡的分布。当被问及理想的财富分配时，92%的人选择了一个更加公平的分配方式。
以下图形展示了本研究的结果：
前面的图形由 Mother Jones 发布，来源为www.motherjones.com/politics/2011/02/income-inequality-in-america-chart-graph，数据来源于www.people.hbs.edu/mnorton/norton%20ariely%20in%20press.pdf。
之前的图表很好地展示了人们对每个五分位的财富分布的感知与现实之间的相对差异。我们将把这些数据转化为一种替代表示方式，就像之前一样，我们可以将这些数据转化为表格表示。
从之前的图表和参考链接的论文中，我得到了以下按五分位划分的大致数据：
| 五分位 | 理想百分比 | 预期百分比 | 实际百分比 |
|---|---|---|---|
| 100th | 32.0% | 58.5% | 84.5% |
| 80th | 22.0% | 20.0% | 11.5% |
| 60th | 21.5% | 12.0% | 3.7% |
| 40th | 14.0% | 6.5% | 0.2% |
| 20th | 10.5% | 3.0% | 0.1% |
让我们取理想和预期的分布，并找到一种方法将它们绘制在我们现有的财富分布图上。我们的柱状图已经表示了不同百分位的相对财富，作为面积的大小。为了使两个数据集可比，我们也应该对这些数据做相同的处理。之前的表格通过已经将数据表示为五个大小相等的组，帮助我们简化了处理，因此我们无需像处理饼图数据时那样应用转换。
然而，让我们利用这个机会，进一步了解在 Quil 中绘制复杂形状，看看能否得到如下所示的数据展示：
该表格提供了我们要通过标签为A、B和C的形状表示的相对面积。为了绘制先前的形状，我们必须计算出高度x、y和z。这些将为我们提供可在图表上绘制的坐标。
区域A、B和C的宽度是w。因此，x和w的乘积将等于A的面积：
由此可得，x的高度仅为A的面积除以w。Y则稍微复杂一点，但也不算太复杂。B的三角形部分的面积等于：
因此：
我们可以通过相同的方式计算z：
扩展我们的定义，得到以下z的方程：
如果我们假设w为 1（所有我们的五分位宽度相等），那么我们可以得到以下方程，适用于任意数量的区间：
这可以通过一个简单的递归函数来表示。我们的第一个比例将被赋值为x。随后可以通过以下方式计算其他值：
(defn area-points [proportions]
(let [f (fn [prev area]
(-> (- area prev)
(* 2)
(+ prev)))
sum (reduce + proportions)]
(->> (reductions f (first proportions) proportions)
(map #(/ % sum)))))
reductions函数的行为与reduce完全相同，但会保留我们计算的中间步骤。我们不会只得到一个值，而是会得到一个值序列，对应于我们y-坐标的（比例）高度。
绘制复杂形状
前面定义的area-points函数将为我们提供一系列点以便绘制。然而，我们还没有涉及 Quil 中的函数，这些函数将允许我们绘制它们。要绘制直线，我们可以使用q/line函数。线函数将接受起点和终点坐标，并在它们之间绘制一条直线。我们可以通过这种方式构建一个面积图，但它没有填充。线条仅描述轮廓，我们不能像使用q/rect绘制直方图时那样构建带有填充色的形状。为了给形状填充颜色，我们需要一个一个顶点地构建它们。
要用 Quil 构建任意复杂的形状，首先调用q/begin-shape。这是一个有状态的函数，它让 Quil 知道我们想开始构建一系列顶点。随后调用的q/vertex将与我们正在构建的形状关联。最后，调用q/end-shape将完成形状的构建。我们将根据当前绘图上下文中指定的描边和填充样式来绘制它。
让我们通过绘制一些使用上一节定义的area-points函数的测试形状，看看它是如何工作的：
(defn plot-area [proportions px py width height]
(let [ys (area-points proportions)
points (map vector (range) ys)
x-scale (/ width (dec (count ys)))
y-scale (/ height (apply max ys))]
(q/stroke 0)
(q/fill 200)
(q/begin-shape)
(doseq [[x y] points]
(q/vertex (+ px (* x x-scale))
(- py (* y y-scale))))
(q/end-shape)))
(defn ex-10-16 []
(let [expected [3 6.5 12 20 58.5]
width 640
height 480
setup (fn []
(q/background 255)
(plot-area expected 0 height width height))]
(q/sketch :setup setup :size [width height])))
这个示例使用之前定义的area-points函数绘制了[3 6.5 12 20 58.5]序列。这是美国“预期”财富分布数据表中列出的百分比值序列。plot-area函数会调用begin-shape，遍历由area-points返回的*ys*序列，并调用end-shape。结果如下：
这还不是我们想要的效果。虽然我们要求填充形状，但我们并没有描述完整的待填充形状。Quil 不知道我们想如何封闭形状，因此它只是从最后一个点画到第一个点，横跨了我们图表的对角线。幸运的是，这个问题可以通过确保图表的底部两个角都有点来轻松解决：
(defn plot-full-area [proportions px py width height]
(let [ys (area-points proportions)
points (map vector (range) ys)
x-scale (/ width (dec (count ys)))
y-scale (/ height (apply max ys))]
(q/stroke 0)
(q/fill 200)
(q/begin-shape)
(q/vertex 0 height)
(doseq [[x y] points]
(q/vertex (+ px (* x x-scale))
(- py (* y y-scale))))
(q/vertex width height)
(q/end-shape)))
(defn ex-10-17 []
(let [expected [3 6.5 12 20 58.5]
width 640
height 480
setup (fn []
(q/background 255)
(plot-full-area expected 0 height width height))]
(q/sketch :setup setup :size [width height])))
plot-full-area函数在遍历*ys*序列之前和之后都会多次调用vertex。这些指定的点确保在调用end-shape之前，形状已完全描述。结果如下图所示：
这已经好多了，开始看起来像一个面积图了。在下一节中，我们将讨论如何使用曲线来描述更复杂的形状。虽然曲线不是我们面积图的必需元素，但它会让结果看起来更有吸引力。
绘制曲线
面积图看起来不错，但我们可以通过使用 Quil 的样条曲线来去除那些尖锐的角。与其通过添加顶点来构建形状，我们可以调用q/curve-vertex来平滑边缘之间的连接。
q/curve-vertex函数实现了一种已知的曲线绘制方法，叫做 Catmull-Rom 样条曲线。要绘制一条曲线，我们必须至少指定四个顶点：第一个和最后一个顶点将作为控制点，曲线将在中间的两个顶点之间绘制。
我们通过下面的图示来可视化 Catmull-Rom 样条曲线的工作原理，该图显示了由点a、b、c和d指定的路径：
在c点处的曲线切线与X平行：由a和b两点所描述的直线；在b点处的曲线切线与Y平行：由c和d两点所描述的直线。因此，要绘制一条曲线，我们需要确保在曲线的起点和终点添加这些额外的控制点。每个控制点都通过curve-vertex添加，我们会在迭代点之前先调用一次，然后在结束时再调用一次：
(defn smooth-curve [xs ys]
(let [points (map vector xs ys)]
(apply q/curve-vertex (first points))
(doseq [point points]
(apply q/curve-vertex point))
(apply q/curve-vertex (last points))))
现在我们已经定义了一个smooth-curve函数，我们将在接下来的两个函数中使用它，分别是smooth-stroke和smooth-area：
(defn smooth-stroke [xs ys]
(q/begin-shape)
(q/vertex (first xs) (first ys))
(smooth-curve (rest xs) (rest ys))
(q/end-shape))
(defn smooth-area [xs ys]
(q/begin-shape)
(q/vertex (first xs) (first ys))
(smooth-curve (rest xs) (rest ys))
(q/vertex (last xs) (first ys))
(q/end-shape))
smooth-stroke函数将通过为每个xs和ys创建顶点来绘制定义的形状。smooth-area函数通过闭合形状并避免之前看到的填充与形状对角线交叉的情况，扩展了这一点。将这两个函数结合起来的是plot-curve，该函数接受要绘制的xs和ys，以及用于绘制的填充颜色、边框颜色和边框粗细：
(defn plot-curve [xs ys fill-color
stroke-color stroke-weight]
(let [points (map vector xs ys)]
(q/no-stroke)
(q/fill fill-color)
(smooth-area xs ys)
(q/no-fill)
(q/stroke stroke-color)
(q/stroke-weight stroke-weight)
(smooth-stroke xs ys)))
让我们在之前绘制的相同预期值序列上调用plot-curve函数，并比较差异：
(defn plot-smooth-area [proportions px py width height]
(let [ys (cons 0 (area-points proportions))
points (map vector (range) ys)
x-scale (/ width (dec (count ys)))
y-scale (/ height (apply max ys) -1)]
(plot-curve (map (point->px px x-scale) (range (count ys)))
(map (point->px py y-scale) ys)
(q/color 200)
(q/color 0) 2)))
(defn ex-10-18 []
(let [expected [3 6.5 12 20 58.5]
width 640
height 480
setup (fn []
(q/background 255)
(plot-smooth-area expected 0 height
width height))]
(q/sketch :setup setup :size [width height])))
这个示例生成了如下图像：
曲线的效果微妙，但它为我们的图表提供了润色，否则图表会显得不完整。前面的图表展示了 Norton 和 Ariely 研究中预期的财富分布。在将其与之前创建的实际财富分布图结合之前，让我们看看它如何与同一研究中的理想财富分布结合。
绘制复合图表
之前的描述展示了如何创建一个适应区域的单一曲线图。正如我们定义的那样，plot-smooth-area函数将在我们绘制的每个区域中填充指定的高度。从绘图的角度看，这样做是合理的，但在尝试绘制两个可比较图表时，这就不太合适了：我们需要确保它们使用相同的缩放比例。
在接下来的代码块中，我们将根据两个图形中较大的一个来计算缩放比例，然后使用该比例绘制两个图形。这样可以确保我们绘制的所有系列彼此之间可以进行比较。组合图表将填满我们分配给它的宽度和高度：
(defn plot-areas [series px py width height]
(let [series-ys (map area-points series)
n-points (count (first series-ys))
x-scale (point->px px (/ width (dec n-points)))
xs (map x-scale (range n-points))
y-max (apply max (apply concat series-ys))
y-scale (point->px py (/ height y-max -1))]
(doseq [ys series-ys]
(plot-curve (cons (first xs) xs)
(map y-scale (cons 0 ys))
(q/color 255 100)
(q/color 255 200) 3))))
(defn ex-10-19 []
(let [expected [3 6.5 12 20 58.5]
ideal [10.5 14 21.5 22 32]
width 640
height 480
setup (fn []
(q/background 100)
(plot-areas [expected ideal] 0 height
width height))]
(q/sketch :setup setup :size [width height])))
我们使用 plot-areas 函数绘制了 expected 和 ideal 两个系列，且通过 background 函数将背景设为更深的颜色。在调用 plot-curve 时，我们指定使用半透明的白色作为填充颜色。下图展示了结果：
为了将这个图表与之前创建的实际数据图表结合，我们只需要调整它的尺度以便匹配。此图表右上角的最高点对应于 5% 的概率密度。我们实际图表上的 96-99^(th) 百分位代表总量的 7.5%，各自对应其图表。这意味着我们需要将之前的图表绘制为现有图表的 2/3 高度，才能使坐标轴可比。现在我们来做这个，并在此过程中为两个新系列添加标签：
(defn draw-expected-ideal []
(let [expected [3 6.5 12 20 58.5]
ideal [10.5 14 21.5 22 32]]
(plot-areas [expected ideal]
plot-x
(+ plot-y plot-height)
plot-width
(* (/ plot-height 0.075) 0.05))
(q/text-size 20)
(emboss-text "EXPECTED" 400 430)
(emboss-text "IDEAL" 250 430)))
最后，我们调用草图中的 draw-expected-ideal 函数以及之前定义的其他函数：
(defn ex-10-20 []
(let [size [960 540]]
(q/sketch :size size
:setup #((draw-shapes)
(draw-expected-ideal)
(draw-banknotes)
(draw-axis-labels)
(draw-title)))))
完成的结果在下图中展示：
希望你会同意，最终的图表既美观又富有信息性。最重要的是，我们是通过实际数据的绘图指令来生成这张图表。如果这张图表是手工制作的，其完整性将难以建立。
输出为 PDF
所有元素组合在一起，最终形成了可能被打印出来的图形。我们提供的绘图指令是基于矢量的——而非基于像素的——因此它将按需要的任何分辨率进行缩放而不失真。
与我们用 save 方法将直方图输出为基于像素的格式不同，我们将输出为 PDF 格式。PDF 格式能够保持我们的艺术作品的可扩展性，并且允许我们在任何所需的分辨率下输出。为此，我们通过传递 :pdf 关键字以及 :output-file 路径来配置草图以使用 PDF 渲染器。
(defn ex-10-21 []
(let [size [960 540]]
(q/sketch :size size
:setup #((draw-shapes)
(draw-expected-ideal)
(draw-banknotes)
(draw-axis-labels)
(draw-title))
:renderer :pdf
:output-file "wealth-distribution.pdf")))
最终的示例将把完成的 PDF 文件输出到项目目录的根目录。
总结
在本章中，我们看到通过使用简单的可视化—仅使用彩色矩形—如何从数据中获得有价值的洞察，并且如何结合 Clojure 核心函数与 Quil 的绘图 API，使我们能够生成传递信息的强大图形。
我们通过 Quil 库实现了这一切。Quil 的功能远不止我们在此展示的内容：它支持交互式动画，支持 ClojureScript 输出到网页，并且还可以进行 3D 渲染。可视化是一个庞大的主题，我们在本章中只能通过一些示例来激发你对它的兴趣。通过展示即使是使用矩形、曲线和 SVG 的基础绘图指令如何组合成强大的图形，我们希望能激发你创造自定义可视化的可能性。
这是《Clojure 数据科学》的最后一章。请务必访问本书的网站 clojuredatascience.com，获取更多关于所涵盖主题的信息和进一步阅读的链接。我们打算为数据科学家，尤其是 Clojure 程序员，提供一个持续更新的资源。
使用一种库也在迅速发展的编程语言，来传达数据科学这个多样且迅速发展的领域的广度和深度，确实是一个雄心勃勃的任务。尽管如此，我们希望《Clojure 数据科学》能让你对统计学、机器学习和大数据处理的一些基本概念有所了解。这些概念基础应该能为你提供帮助，即使技术选项——甚至可能是你的编程语言选择——未来还会继续演变。
