原文：
annas-archive.org/md5/4e1b7010caf1fbe3188c9c515fe244d4译者：飞龙
第三章：基于历史和实时数据的高频比特币价格预测
比特币是一种全球性的加密货币和数字支付系统，被认为是第一种去中心化的数字货币，因为该系统无需中央存储库或单一管理员即可运行。近年来，它在全球范围内获得了极大的关注和流行。
在本章中，我们将展示如何使用 Scala、Spark ML、Cryptocompare API 和比特币历史（以及实时）数据来预测未来一周、一个月等的价格，这将帮助我们为在线加密货币交易做出自动化决策。此外，我们还将展示如何生成一个简单的信号用于在线比特币交易。
简要来说，在这个端到端的项目中，我们将学习以下主题：
-
比特币，加密货币与在线交易
-
历史和实时价格数据收集
-
原型的高级管道
-
使用梯度提升树回归预测比特币价格
-
使用 Scala play 框架进行的演示预测和信号生成
-
未来展望——使用相同的技术处理其他数据集
比特币，加密货币与在线交易
比特币，作为按发行时间和市值（截至 2017 年 12 月）排名第一的加密货币，因其易于开始交易、能够保持伪匿名性以及历史上的剧烈增长（参见表 1和图 1获取一些统计数据）而吸引了投资者和交易员。这吸引了长期投资者；而其高度的波动性也吸引了日内交易者。
然而，很难预测比特币的长期价值，因为比特币背后的价值较为抽象。其价格主要反映市场认知，并且高度依赖新闻、法规、政府和银行的合作、平台的技术问题（如交易费用和区块大小）、机构投资者将比特币纳入其投资组合的兴趣等因素：
图 1：比特币及其剧烈的价格上涨
然而，从短期来看，比特币价格是市场活动的副产品，通常发生在一个被称为交易所的平台上（其中最著名的交易所包括 Bitstamp、Coinbase、Kraken 和 Bitfinex）。用户注册并完成KYC（了解你的客户）程序后，可以在其中交易比特币，兑换成法币如美元和欧元，也可以交易其他加密货币，称为山寨币（例如以太坊、莱特币和达世币等较为知名的山寨币）：
表 1 - 比特币历史价格变动
交易所维护订单簿——列出了所有买卖订单、其数量和价格——并在买卖双方匹配时执行交易。此外，交易所还保持并提供有关交易状态的统计信息，通常以 OCHL（开盘、最高、最低、收盘）和成交量表示，涵盖交易对的两种货币。对于这个项目，我们将使用 BTC/USD 加密货币对。
这些数据按周期汇总展示，从秒到天，甚至到月。为专业交易员和机构收集比特币数据的专用服务器也在运行。虽然不可能期望所有订单数据都能免费获取，但其中一些数据对公众开放并且可以使用。
最先进的比特币自动化交易
在传统证券的世界中，比如公司的股票，过去是由人类来进行分析、预测股票价格并进行交易的。今天，机器学习（ML）的发展和数据的日益丰富几乎已将人类从高频交易中排除，因为普通人无法捕捉并处理所有数据，而且情绪会影响决策；因此，这一领域现在主要由投资机构的自动化交易系统主导。
当前，比特币交易的交易量相对于传统交易所较低；金融机构传统上谨慎且风险厌恶，尚未涉足比特币交易（至少，尚未公开知晓）。其中一个原因是高费用和对加密货币法规的不确定性。
所以今天，大多数人购买和出售比特币，伴随着所有与非理性行为相关的后果，但也有人尝试自动化比特币交易。最著名的一次尝试是在 MIT 的一篇论文中提出的，另一次是在 2014 年由斯坦福研究人员发布的。很多事情发生了变化，考虑到这三年中比特币价格的大幅上涨，任何只是买入并持有的人都会对结果感到满意：
图 2：比特币买卖订单（截至 2017 年 11 月）
当然，确实有一些交易员使用机器学习进行交易，这类应用看起来前景广阔。到目前为止，从研究论文中确定的最佳方法如下：
训练
使用订单簿数据，而不是衍生的 OHLC + 成交量数据。因此，训练和预测时，使用的数据看起来是这样的：
-
将数据划分为一定
大小的时间序列（大小是一个需要调整的参数）。 -
将时间序列数据聚类为
K个聚类（K是一个需要调整的参数）。假设会出现一些具有自然趋势的聚类（如价格剧烈下降/上升等）。 -
对于每个聚类，训练回归模型和分类器，分别预测价格和价格变化。
预测
这种方法考虑了具有特定窗口大小的最新时间序列，并对模型进行训练。然后，它对数据进行如下分类：
-
采用用于训练的最新时间序列和窗口大小
-
进行分类——它属于哪个聚类？
-
使用该聚类的 ML 模型预测价格或价格变化
这个解决方案可以追溯到 2014 年，但它仍然提供了一定程度的鲁棒性。由于需要识别许多参数，并且没有容易获得的订单簿历史数据，本项目使用了更简单的方法和数据集。
原型的高级数据管道
本章的目标是开发一个原型系统，该系统将预测比特币价格的短期变化，使用历史数据训练算法，并使用实时数据预测和选择表现更好的算法。在本项目范围内，并不尝试预测实际的美元价格，而仅仅是预测其是否会上涨。因为比特币价格在某种程度上实际上并不只是价格本身，而是市场预期。这可以看作是交易者行为中的模式，从更高层次上讲，表现为历史价格本身。
图 3：原型的高级数据管道
当然，比特币有一个客观的价格；矿工愿意出售比特币以获取利润。因此，基准价格可以通过了解所有矿工为开采比特币必须支付的费用来估算，但这超出了本项目的范围。
从这个角度来看，与其试图预测美元价格，不如寻找价格上涨、下跌或保持不变的趋势，并根据这些趋势采取行动。第二个目标是构建一个实验工具，允许我们尝试不同的价格预测方法，并能轻松在实际数据上进行评估。代码必须具有灵活性、稳健性，并且易于扩展。
因此，总结来说，系统有三个主要组成部分：
-
用于将历史数据预处理成所需格式的 Scala 脚本
-
Scala 应用程序用于训练机器学习模型
-
Scala 网络服务预测未来价格
历史和实时价格数据收集
如前所述，我们将同时使用历史数据和实时数据。我们将使用来自 Kaggle 的比特币历史价格数据。对于实时数据，将使用 Cryptocompare API。
历史数据收集
为了训练机器学习算法，有一个公开的比特币历史价格数据数据集（版本 10），可以在 Kaggle 上找到。该数据集包含来自多个交易所的 BTC-USD 对的 1 分钟 OHLC 数据，可以从www.kaggle.com/mczielinski/bitcoin-historical-data/下载。
在项目初期，大多数数据可用的时间范围是从 2012 年 1 月 1 日到 2017 年 5 月 31 日；但是对于 Bitstamp 交易所，数据可用时间截至 2017 年 10 月 20 日（对于 Coinbase 也是如此，但该数据集稍后才提供）：
图 4：Kaggle 上的比特币历史数据集
请注意，您需要是注册用户并已登录才能下载该文件。我们使用的文件是bitstampUSD_1-min_data_2012-01-01_to_2017-10-20.csv*。*现在，让我们来获取我们拥有的数据。它有八列：
-
时间戳：自 1970 年 1 月 1 日以来的秒数。对于第一行，它是 1,325,317,920，对于第二行也是 1,325,317,920。（完整性检查！它们之间的差异是 60 秒）。
-
开盘价：该时间区间的开盘价格。为 4.39 美元。因此，它是发生在时间戳（第一个行中的时间戳为 1,325,317,920）之后的第一次交易的价格。
-
收盘价：该时间区间的收盘价格。
-
最高价：在该区间内执行的所有订单中的最高价格。
-
最低价：与最高价相同，但它是最低价格。
-
交易量 _(BTC)：在该时间区间内转移的所有比特币的总和。因此，选择该区间内发生的所有交易并将每笔交易的 BTC 值相加。
-
交易量 _(货币)：转移的所有美元的总和。
-
加权价格：这是根据 BTC 和 USD 的交易量得出的。通过将所有交易的美元数额除以所有比特币的交易量，我们可以得出该分钟的 BTC 加权平均价格。所以
加权价格=交易量 _(货币)/交易量 _(BTC)。
数据科学流程中最重要的部分之一是数据预处理——清理数据集并将其转换为适合我们需求的形式，这发生在数据收集之后（在某种意义上，数据收集是外包的；我们使用的是他人收集的数据）。
历史数据转换为时间序列
从我们的目标出发——预测价格变化的方向——我们可能会问自己，实际的美元价格是否有助于实现这一目标？ 历史上，比特币的价格通常是上涨的，所以如果我们尝试拟合线性回归，它将显示出进一步的指数增长（长期内是否会如此，还需要观察）。
假设和设计选择
这个项目的一个假设是：无论我们是在考虑 2016 年 11 月比特币的交易价格大约为 700 美元，还是 2017 年 11 月的价格在 6500 至 7000 美元区间，交易模式是相似的。现在，我们还有几个其他的假设，如下所述：
-
假设一：根据前面的说法，我们可以忽略实际价格，而是关注价格变化。作为这一变化的度量，我们可以取开盘价和收盘价之间的差值。如果为正，表示该分钟内价格上涨；如果为负，表示价格下跌；如果差值为 0，表示价格保持不变。
在下图中，我们可以看到，观察的第一分钟 Delta 为 -1.25，第二分钟为 -12.83，第三分钟为 -0.23。有时候，开盘价与上一分钟的收盘价差异较大（尽管在观察的三分钟内，Delta 都是负值，但在第三分钟，显示的价格实际上比收盘价略高）。但这种情况并不常见，通常情况下开盘价与上一分钟的收盘价变化不大。
-
假设二：接下来需要考虑的是在黑箱环境下预测价格变化。我们不使用新闻、Twitter 动态等其他信息来源来预测市场如何对它们做出反应。这是一个更为高级的话题。我们唯一使用的数据是价格和成交量。为了简化原型，我们可以只关注价格并构建时间序列数据。
时间序列预测是基于该参数过去的值来预测该参数的未来值。最常见的一个例子是温度预测。尽管有许多超级计算机使用卫星和传感器数据来预测天气，但简单的时间序列分析也能得出一些有价值的结果。例如，我们可以基于 T 时刻、T-60 秒、T-120 秒等时刻的价格来预测 T+60 秒时的价格。
-
假设三：数据集中并非所有数据都有价值。前 60 万条记录不具备信息量，因为价格变化稀少，成交量较小。这可能会影响我们训练的模型，从而导致最终结果变差。因此，数据集中前 60 万条记录会被剔除。
-
假设四：我们需要
标签化我们的数据，以便使用监督式机器学习算法。这是最简单的措施，且无需担心交易费用。
数据预处理
考虑到数据准备的目标，选择了 Scala 作为一种简便且交互性强的数据处理方式：
val priceDataFileName: String = "bitstampUSD_1-min_data_2012-01-01_to_2017-10-20.csv"
val spark = SparkSession
.builder()
.master("local[*]")
.config("spark.sql.warehouse.dir", "E:/Exp/")
.appName("Bitcoin Preprocessing")
.getOrCreate()
val data = spark.read.format("com.databricks.spark.csv").option("header", "true").load(priceDataFileName)
data.show(10)
>>>
图 5：比特币历史价格数据集一瞥
println((data.count(), data.columns.size))
(3045857, 8)
在前面的代码中，我们从 Kaggle 下载的文件中加载数据，并查看其内容。数据集中有3045857行数据，8列，之前已经描述过。接着我们创建了Delta列，包含了收盘价与开盘价之间的差值（也就是只考虑那些已经开始有意义的交易数据）：
val dataWithDelta = data.withColumn("Delta", data("Close") - data("Open"))
以下代码通过将Delta值为正的行标记为 1，其他行标记为0，为我们的数据打上标签：
import org.apache.spark.sql.functions._
import spark.sqlContext.implicits._
val dataWithLabels = dataWithDelta.withColumn("label", when($"Close" - $"Open" > 0, 1).otherwise(0))
rollingWindow(dataWithLabels, 22, outputDataFilePath, outputLabelFilePath)
这段代码将原始数据集转换为时间序列数据。它将WINDOW_SIZE行（在此实验中为22）的 Delta 值合并成一行。这样，第一行包含从t0到t21的 Delta 值，第二行包含从t1到t22的 Delta 值。然后我们创建相应的标签数组（1或0）。
最后，我们将 X 和 Y 保存到文件中，其中 612000 行数据被从原始数据集中截取；22 表示滚动窗口大小，2 个类别表示标签是二进制的 0 和 1：
val dropFirstCount: Int = 612000
def rollingWindow(data: DataFrame, window: Int, xFilename: String, yFilename: String): Unit = {
var i = 0
val xWriter = new BufferedWriter(new FileWriter(new File(xFilename)))
val yWriter = new BufferedWriter(new FileWriter(new File(yFilename)))
val zippedData = data.rdd.zipWithIndex().collect()
System.gc()
val dataStratified = zippedData.drop(dropFirstCount)//slice 612K
while (i < (dataStratified.length - window)) {
val x = dataStratified
.slice(i, i + window)
.map(r => r._1.getAsDouble).toList
val y = dataStratified.apply(i + window)._1.getAsInteger
val stringToWrite = x.mkString(",")
xWriter.write(stringToWrite + "n")
yWriter.write(y + "n")
i += 1
if (i % 10 == 0) {
xWriter.flush()
yWriter.flush()
}
}
xWriter.close()
yWriter.close()
}
在前面的代码段中：
val outputDataFilePath: String = "output/scala_test_x.csv"
val outputLabelFilePath: String = "output/scala_test_y.csv"
通过 Cryptocompare API 获取实时数据
对于实时数据，使用 Cryptocompare API（www.cryptocompare.com/api/#），更具体来说是 HistoMinute（www.cryptocompare.com/api/#-api-data-histominute-），它使我们可以访问最多过去七天的 OHLC 数据。API 的细节将在专门讨论实现的部分中说明，但 API 响应与我们的历史数据集非常相似，这些数据是通过常规的 HTTP 请求获取的。例如，来自 min-api.cryptocompare.com/data/histominute?fsym=BTC&tsym=USD&limit=23&aggregate=1&e=Bitstamp 的简单 JSON 响应具有以下结构：
{
"Response":"Success",
"Type":100,
"Aggregated":false,
"Data":
[{"time":1510774800,"close":7205,"high":7205,"low":7192.67,"open":7198, "volumefrom":81.73,"volumeto":588726.94},
{"time":1510774860,"close":7209.05,"high":7219.91,"low":7205,"open":7205, "volumefrom":16.39,"volumeto":118136.61},
... (other price data)
],
"TimeTo":1510776180,
"TimeFrom":1510774800,
"FirstValueInArray":true,
"ConversionType":{"type":"force_direct","conversionSymbol":""}
}
通过 Cryptocompare HistoMinute，我们可以获取每分钟的历史数据中的 open、high、low、close、volumefrom 和 volumeto。这些数据仅保存 7 天；如果需要更多数据，可以使用按小时或按日的路径。如果数据不可用（因为该货币没有在指定的货币中进行交易），它会使用 BTC 转换：
图 6：通过 Cryptocompare HistoMinute 获取的开盘价、最高价、最低价、收盘价和交易量值
现在，以下方法获取正确格式的 Cryptocompare API URL（www.cryptocompare.com/api/#-api-data-histominute-），这是一个完全形成的 URL，指定了所有的参数，如货币、限制和聚合等。它最终返回一个 future 对象，该对象将解析响应体并转化为数据模型，价格列表将在上层进行处理：
import javax.inject.Inject
import play.api.libs.json.{JsResult, Json}
import scala.concurrent.Future
import play.api.mvc._
import play.api.libs.ws._
import processing.model.CryptoCompareResponse
class RestClient @Inject() (ws: WSClient) {
def getPayload(url : String): Future[JsResult[CryptoCompareResponse]] = {
val request: WSRequest = ws.url(url)
val future = request.get()
implicit val context = play.api.libs.concurrent.Execution.Implicits.defaultContext
future.map {
response => response.json.validate[CryptoCompareResponse]
}
}
}
在前面的代码段中，CryptoCompareResponse 类是 API 的模型，它需要以下参数：
-
响应 -
类型 -
聚合 -
数据 -
FirstValueInArray -
TimeTo -
TimeFrom
现在，它具有以下签名：
case class CryptoCompareResponse(Response : String,
Type : Int,
Aggregated : Boolean,
Data : List[OHLC],
FirstValueInArray : Boolean,
TimeTo : Long,
TimeFrom: Long)
object CryptoCompareResponse {
implicit val cryptoCompareResponseReads = Json.reads[CryptoCompareResponse]
}
再次说明，前两个代码段中的 open-high-low-close（也称为 OHLC）是用于与 CryptoAPI 响应的 data 数组内部映射的模型类。它需要这些参数：
-
Time：时间戳（以秒为单位），例如1508818680。 -
Open：给定分钟间隔的开盘价。 -
High：最高价。 -
Low：最低价。 -
Close：区间结束时的价格。 -
Volumefrom：from货币的交易量。在我们的例子中是 BTC。 -
Volumeto：to货币的交易量，在我们的例子中是 USD。 -
将
Volumeto除以Volumefrom可以给我们带来 BTC 的加权价格。
现在，它具有以下签名：
case class OHLC(time: Long,
open: Double,
high: Double,
low: Double,
close: Double,
volumefrom: Double,
volumeto: Double)
object OHLC {
implicit val implicitOHLCReads = Json.reads[OHLC]
}
预测模型训练
在项目中，prediction.training包下有一个名为TrainGBT.scala的 Scala 对象。在启动之前，你需要指定/更改以下四个设置：
-
在代码中，你需要设置
spark.sql.warehouse.dir，将其指向你电脑中某个有足够存储空间的实际路径：set("spark.sql.warehouse.dir", "/home/user/spark") -
RootDir是主文件夹，所有文件和训练模型将存储在此处：rootDir = "/home/user/projects/btc-prediction/" -
确保
x文件名与前一步 Scala 脚本生成的文件名匹配：x = spark.read.format("com.databricks.spark.csv ").schema(xSchema).load(rootDir + "scala_test_x.csv") -
确保
y文件名与 Scala 脚本生成的文件名匹配：y_tmp=spark.read.format("com.databricks.spark.csv").schema(ySchema).load(rootDir + "scala_test_y.csv")
训练的代码使用了 Apache Spark 的 ML 库（以及为其所需的库）来训练分类器，这意味着这些库必须出现在你的class路径中，才能运行它。最简单的做法是（由于整个项目使用 SBT），从项目根文件夹运行，输入 sbt run-main prediction.training.TrainGBT，这将解析所有依赖并启动训练。
根据迭代次数和深度的不同，训练模型可能需要几个小时。现在，让我们通过梯度提升树模型的示例来看看训练是如何进行的。首先，我们需要创建一个SparkSession对象：
val spark = SparkSession
.builder()
.master("local[*]")
.config("spark.sql.warehouse.dir", ""/home/user/spark/")
.appName("Bitcoin Preprocessing")
.getOrCreate()
然后，我们定义x和y的数据模式。我们将列重命名为t0至t21，表示它是时间序列数据：
val xSchema = StructType(Array(
StructField("t0", DoubleType, true),
StructField("t1", DoubleType, true),
StructField("t2", DoubleType, true),
StructField("t3", DoubleType, true),
StructField("t4", DoubleType, true),
StructField("t5", DoubleType, true),
StructField("t6", DoubleType, true),
StructField("t7", DoubleType, true),
StructField("t8", DoubleType, true),
StructField("t9", DoubleType, true),
StructField("t10", DoubleType, true),
StructField("t11", DoubleType, true),
StructField("t12", DoubleType, true),
StructField("t13", DoubleType, true),
StructField("t14", DoubleType, true),
StructField("t15", DoubleType, true),
StructField("t16", DoubleType, true),
StructField("t17", DoubleType, true),
StructField("t18", DoubleType, true),
StructField("t19", DoubleType, true),
StructField("t20", DoubleType, true),
StructField("t21", DoubleType, true))
)
然后，我们读取为模式定义的文件。为了方便起见，我们在 Scala 中生成了两个单独的文件来存储数据和标签，因此这里我们需要将它们合并为一个 DataFrame：
import spark.implicits._
val y = y_tmp.withColumn("y", 'y.cast(IntegerType))
import org.apache.spark.sql.functions._
val x_id = x.withColumn("id", monotonically_increasing_id())
val y_id = y.withColumn("id", monotonically_increasing_id())
val data = x_id.join(y_id, "id")
下一步是 Spark 要求的——我们需要将特征向量化：
val featureAssembler = new VectorAssembler()
.setInputCols(Array("t0", "t1", "t2", "t3",
"t4", "t5", "t6", "t7",
"t8", "t9", "t10", "t11",
"t12", "t13", "t14", "t15",
"t16", "t17", "t18", "t19",
"t20", "t21"))
.setOutputCol("features")
我们将数据随机拆分为训练集和测试集，比例为 75%对 25%。我们设置种子，以确保每次运行训练时，数据拆分相同：
val Array(trainingData,testData) = dataWithLabels.randomSplit(Array(0.75, 0.25), 123)
然后，我们定义模型。它告诉我们哪些列是特征，哪些是标签，同时设置参数：
val gbt = new GBTClassifier()
.setLabelCol("label")
.setFeaturesCol("features")
.setMaxIter(10)
.setSeed(123)
创建一个pipeline步骤——特征向量组合和 GBT 运行：
val pipeline = new Pipeline()
.setStages(Array(featureAssembler, gbt))
定义评估器函数——模型如何知道自己是否表现良好。由于我们只有两个不平衡的类别，准确率是一个不好的衡量标准；ROC 曲线下的面积更合适：
val rocEvaluator = new BinaryClassificationEvaluator()
.setLabelCol("label")
.setRawPredictionCol("rawPrediction")
.setMetricName("areaUnderROC")
使用 K 折交叉验证来避免过拟合；每次迭代会去除五分之一的数据，利用其余数据训练模型，然后在这一五分之一的数据上进行测试：
val cv = new CrossValidator()
.setEstimator(pipeline)
.setEvaluator(rocEvaluator)
.setEstimatorParamMaps(paramGrid)
.setNumFolds(numFolds)
.setSeed(123)
val cvModel = cv.fit(trainingData)
在获得训练后的模型后（根据迭代次数和我们在paramGrid中指定的参数，可能需要一个小时或更长时间），我们接着在测试数据上计算预测结果：
val predictions = cvModel.transform(testData)
此外，评估预测的质量：
val roc = rocEvaluator.evaluate(predictions)
训练后的模型将被保存，以便预测服务后续使用：
val gbtModel = cvModel.bestModel.asInstanceOf[PipelineModel]
gbtModel.save(rootDir + "__cv__gbt_22_binary_classes_" + System.nanoTime() / 1000000 + ".model")
总结来说，模型训练的代码如下所示：
import org.apache.spark.{ SparkConf, SparkContext }
import org.apache.spark.ml.{ Pipeline, PipelineModel }
import org.apache.spark.ml.classification.{ GBTClassificationModel, GBTClassifier, RandomForestClassificationModel, RandomForestClassifier}
import org.apache.spark.ml.evaluation.{BinaryClassificationEvaluator, MulticlassClassificationEvaluator}
import org.apache.spark.ml.feature.{IndexToString, StringIndexer, VectorAssembler, VectorIndexer}
import org.apache.spark.ml.tuning.{CrossValidator, ParamGridBuilder}
import org.apache.spark.sql.types.{DoubleType, IntegerType, StructField, StructType}
import org.apache.spark.sql.SparkSession
object TrainGradientBoostedTree {
def main(args: Array[String]): Unit = {
val maxBins = Seq(5, 7, 9)
val numFolds = 10
val maxIter: Seq[Int] = Seq(10)
val maxDepth: Seq[Int] = Seq(20)
val rootDir = "output/"
val spark = SparkSession
.builder()
.master("local[*]")
.config("spark.sql.warehouse.dir", ""/home/user/spark/")
.appName("Bitcoin Preprocessing")
.getOrCreate()
val xSchema = StructType(Array(
StructField("t0", DoubleType, true),
StructField("t1", DoubleType, true),
StructField("t2", DoubleType, true),
StructField("t3", DoubleType, true),
StructField("t4", DoubleType, true),
StructField("t5", DoubleType, true),
StructField("t6", DoubleType, true),
StructField("t7", DoubleType, true),
StructField("t8", DoubleType, true),
StructField("t9", DoubleType, true),
StructField("t10", DoubleType, true),
StructField("t11", DoubleType, true),
StructField("t12", DoubleType, true),
StructField("t13", DoubleType, true),
StructField("t14", DoubleType, true),
StructField("t15", DoubleType, true),
StructField("t16", DoubleType, true),
StructField("t17", DoubleType, true),
StructField("t18", DoubleType, true),
StructField("t19", DoubleType, true),
StructField("t20", DoubleType, true),
StructField("t21", DoubleType, true)))
val ySchema = StructType(Array(StructField("y", DoubleType,
true)))
val x = spark.read.format("csv").schema(xSchema).load(rootDir +
"scala_test_x.csv")
val y_tmp =
spark.read.format("csv").schema(ySchema).load(rootDir +
"scala_test_y.csv")
import spark.implicits._
val y = y_tmp.withColumn("y", 'y.cast(IntegerType))
import org.apache.spark.sql.functions._
//joining 2 separate datasets in single Spark dataframe
val x_id = x.withColumn("id", monotonically_increasing_id())
val y_id = y.withColumn("id", monotonically_increasing_id())
val data = x_id.join(y_id, "id")
val featureAssembler = new VectorAssembler()
.setInputCols(Array("t0", "t1", "t2", "t3", "t4", "t5",
"t6", "t7", "t8", "t9", "t10", "t11",
"t12", "t13", "t14", "t15", "t16",
"t17","t18", "t19", "t20", "t21"))
.setOutputCol("features")
val encodeLabel = udf[Double, String] { case "1" => 1.0 case
"0" => 0.0 }
val dataWithLabels = data.withColumn("label",
encodeLabel(data("y")))
//123 is seed number to get same datasplit so we can tune
params
val Array(trainingData, testData) =
dataWithLabels.randomSplit(Array(0.75, 0.25), 123)
val gbt = new GBTClassifier()
.setLabelCol("label")
.setFeaturesCol("features")
.setMaxIter(10)
.setSeed(123)
val pipeline = new Pipeline()
.setStages(Array(featureAssembler, gbt))
// ***********************************************************
println("Preparing K-fold Cross Validation and Grid Search")
// ***********************************************************
val paramGrid = new ParamGridBuilder()
.addGrid(gbt.maxIter, maxIter)
.addGrid(gbt.maxDepth, maxDepth)
.addGrid(gbt.maxBins, maxBins)
.build()
val cv = new CrossValidator()
.setEstimator(pipeline)
.setEvaluator(new BinaryClassificationEvaluator())
.setEstimatorParamMaps(paramGrid)
.setNumFolds(numFolds)
.setSeed(123)
// ************************************************************
println("Training model with GradientBoostedTrees algorithm")
// ************************************************************
// Train model. This also runs the indexers.
val cvModel = cv.fit(trainingData)
cvModel.save(rootDir + "cvGBT_22_binary_classes_" +
System.nanoTime() / 1000000 + ".model")
println("Evaluating model on train and test data and
calculating RMSE")
// **********************************************************************
// Make a sample prediction
val predictions = cvModel.transform(testData)
// Select (prediction, true label) and compute test error.
val rocEvaluator = new BinaryClassificationEvaluator()
.setLabelCol("label")
.setRawPredictionCol("rawPrediction")
.setMetricName("areaUnderROC")
val roc = rocEvaluator.evaluate(predictions)
val prEvaluator = new BinaryClassificationEvaluator()
.setLabelCol("label")
.setRawPredictionCol("rawPrediction")
.setMetricName("areaUnderPR")
val pr = prEvaluator.evaluate(predictions)
val gbtModel = cvModel.bestModel.asInstanceOf[PipelineModel]
gbtModel.save(rootDir + "__cv__gbt_22_binary_classes_" +
System.nanoTime()/1000000 +".model")
println("Area under ROC curve = " + roc)
println("Area under PR curve= " + pr)
println(predictions.select().show(1))
spark.stop()
}
}
现在让我们看看训练的进展：
>>>
Area under ROC curve = 0.6045355104779828
Area under PR curve= 0.3823834607704922
因此，我们并未获得非常高的准确度，因为最佳 GBT 模型的 ROC 仅为 60.50%。尽管如此，如果我们调整超参数，准确度会更好。
然而，由于时间不足，我并未长时间进行训练迭代，但你绝对应该尝试一下。
Scala Play Web 服务
作为应用框架，Play2 被选择为一个易于配置且强大的框架。与 Spring（另一个流行框架）相比，Play2 在从零开始构建一个小型应用时所需时间更少。Play 自带 Guice 用于依赖注入，并且使用 SBT 作为包管理器：
-
Spark ML：选择 Spark ML 库是因为它是 Java 世界中维护得最好的库之一。许多库本身没有的算法是由第三方开发者实现的，并且可以在 Spark 上进行训练。Spark 的一个缺点是它相对较慢，因为它设计上是分布式的；它使用 Hadoop 并大量写入文件系统。
-
Akka：这使得实现 Actor 模式成为可能——拥有多个独立对象实例并在彼此之间并发传递消息，从而提高了系统的健壮性。
-
Anorm：这是一个基于 JDBC 的 SQL 操作库。Slick 是另一个选择，功能更强大，但由于 Akka 和 Slick 所需库之间的兼容性问题，最终选择了另一个库。
-
H2：一个数据库，作为 Play 和 Ruby-on-Rails 的默认数据库，易于启动，且可以将数据存储在本地数据库文件中，而无需安装数据库服务器。这提供了可移植性并加快了开发速度。在后期阶段，它可以被其他数据库替代，因为 Scala 代码与特定数据库无关；所有配置均在配置层面完成。
通过 Akka Actor 实现并发
通过使用 Akka Scala 库中的 actor 模型来实现并发。Actor 作为独立的实体，可以异步地向其他 Actor 传递消息。在这个项目中，有三个 Actor：SchedulerActor、PredictionActor 和 TraderActor：
-
SchedulerActor：请求价格数据，将其存入数据库，向PredictionActor发送价格消息，接收答案并将其传递给TraderActor。 -
PredictionActor：在接收到带有价格的消息后，使用最佳模型预测下一个价格（此模型需要在application.conf中选择；稍后我们会看到详细信息）。它将带有预测结果的消息传回SchedulerActor，并使用model文件夹中的其他模型对历史数据进行预测，最后使用最新的价格来评估预测。这样的预测结果将存储在数据库中。 -
TraderActor：在接收到预测消息后，使用rules（此时规则非常简单，*如果预测价格上涨则买入，*否则不操作），它将决策写入日志。它还可以向一个 URL 发送 HTTP 请求来触发这个决策。
Web 服务工作流
现在让我们更深入地了解代码是如何执行预测的。如前所述，每隔 60 秒，应用程序会触发从 Cryptocompare 获取数据，将价格存储到数据库中，并运行预测，保存回溯测试结果以验证预测质量。
在这一部分，我们将更深入地了解哪些 Scala 类在这个项目中扮演了重要角色以及它们是如何互相通信的。
JobModule
当应用程序启动时，一切从 JobModule 开始。它配置了 Scheduler 的创建，Scheduler 根据 application.conf 中的设置向 SchedulerActor 发送消息：
class JobModule extends AbstractModule with AkkaGuiceSupport {
def configure(): Unit = {
//configuring launch of price-fetching Actor
bindActorSchedulerActor
bind(classOf[Scheduler]).asEagerSingleton()
}
}
要启用此模块，在 application.conf 中需要添加以下行：
play.modules.enabled += "modules.jobs.JobModule"
Scheduler
Scheduler 从 application.conf 中获取频率常量，并使用 Actor 系统每隔 X 秒向 SchedulerActor 发送 update 消息（消息的内容不重要；SchedulerActor 对任何消息都有反应）：
class Scheduler @Inject()
(val system: ActorSystem, @Named("scheduler-actor") val schedulerActor: ActorRef, configuration: Configuration)(implicit ec: ExecutionContext) {
//constants.frequency is set in conf/application.conf file
val frequency = configuration.getInt("constants.frequency").get
var actor = system.scheduler.schedule(
0.microseconds, //initial delay: whether execution starts immediately after app launch
frequency.seconds, //every X seconds, specified above
schedulerActor,
"update")
}
SchedulerActor
相关的代码部分已经显示并解释过了。现在让我们看看如何获取价格数据：
def constructUrl(exchange: String): String =
{
"https://min-api.cryptocompare.com/data/histominute?fsym=BTC&tsym=USD&limit=23&aggregate=1&e=" + exchange
}
ConstructUrl 返回一个完全构造好的 URL，用于向 Cryptocompare API 发送请求。更多细节请参见与 API 相关的章节：
final val predictionActor = system.actorOf(Props(new PredictionActor(configuration, db)))
final val traderActor = system.actorOf(Props(new TraderActor(ws)))
创建 PredictionActor 和 TraderActor 的实例：
override def receive: Receive = {
Receive 方法在 actor 特质中定义并且必须实现。当有人向这个 actor（在我们这个例子中是 Scheduler）传递消息时，它会被触发：
case _ =>
val futureResponse=restClient.getPayload(constructUrl(exchange))
在前面的代码中，case _ => 表示我们对任何类型和内容的消息做出反应。首先，通过前面指定的 URL 异步调用 Cryptocompare API。这个过程借助 RestClient 完成，它返回一个带有响应 JSON 的 Future。在接收到响应后（在 futureResponse 的 complete 回调中），.json 被映射到自定义的 case 类 CryptoCompareResponse：
case class CryptoCompareResponse(Response: String, Type: Int, Aggregated: Boolean, Data: List[OHLC], FirstValueInArray: Boolean, TimeTo: Long,TimeFrom: Long)
这个 case 类类似于 POJO（Plain Old Java Object），不需要编写构造函数和 getter/setter：
object CryptoCompareResponse {
implicit val cryptoCompareResponseReads = Json.reads[CryptoCompareResponse]
}
这个伴随对象用于将 JSON 映射到这个类中。CryptocompareResponse 对象存储了 API 的输出——OHLC 数据的列表、数据的时间范围以及其他与我们无关的内容。OHLC 类对应实际的价格数据：
case class OHLC(time: Long, open: Double,
high: Double,
low: Double,
close: Double,
volumefrom: Double,
volumeto: Double)
数据准备好后，通过调用 storePriceData(cryptoCompareResponse) 将价格存储到数据库中。首先，它使用 Anorm 的 BatchSQL 批量插入到 PRICE_STAGING 表中，然后根据时间戳去重后再插入到 PRICE 表中，因为我们收到的价格数据可能会有重叠：
val batch = BatchSql(
"""|INSERT INTO PRICE_STAGING(TIMESTAMP,EXCHANGE,PRICE_OPEN,PRICE_CLOSED,VOLUME_BTC,
VOLUME_USD)| VALUES({timestamp}, {exchange}, {priceOpen}, {priceClosed}, {volumeBTC}, {volumeUSD})""".stripMargin,transformedPriceDta.head,transformedPriceDta.tail:_*)
val res: Array[Int] = batch.execute() // array of update count
val reInsert = SQL(
"""
|INSERT INTO PRICE(TIMESTAMP, EXCHANGE, PRICE_OPEN, PRICE_CLOSED, VOLUME_BTC, VOLUME_USD)
|SELECT TIMESTAMP, EXCHANGE, PRICE_OPEN, PRICE_CLOSED, VOLUME_BTC, VOLUME_USD
|FROM PRICE_STAGING AS s
|WHERE NOT EXISTS (
|SELECT *
|FROM PRICE As t
|WHERE t.TIMESTAMP = s.TIMESTAMP
|)
""".stripMargin).execute()
Logger.debug("reinsert " + reInsert)
在存入数据库后，SchedulerActor将 OHLC 数据转换为（时间戳，增量）元组，其中增量是（closePrice-openPrice）。因此，格式适用于机器学习模型。转换后的数据作为消息传递给PredictionActor，并明确等待响应。这是通过使用?操作符实现的。我们向预测actor提出请求：
(predictionActor ? CryptoCompareDTOToPredictionModelTransformer.tranform(cryptoCompareResponse)).mapTo[CurrentDataWithShortTermPrediction].map {
它的响应被映射到CurrentDataWithShortTermPrediction类，并通过!操作符传递给TraderActor。与?不同，!操作符不要求响应：
predictedWithCurrent =>
traderActor ! predictedWithCurrent}
这基本上是SchedulerActor的操作流程。我们从 Cryptocompare API 读取数据，存储到数据库中，发送给PredictionActor并等待其响应。然后我们将其响应转发给TraderActor。
现在让我们看看PredictionActor内部发生了什么：
PredictionActor和预测步骤
该 Scala web 应用程序每分钟从 Cryptocompare API 获取 Bitstamp 交易所最新的比特币价格数据，使用训练好的机器学习分类器预测下一分钟的价格变动方向，并通知用户决策结果。
现在，要启动它，可以从项目目录中使用sbt run（或者在需要时使用$ sudo sbt run）。现在让我们看看application.conf文件的内容：
# This is the main configuration file for the application.
# Secret key
# The secret key is used to secure cryptographics functions.
# If you deploy your application to several instances be sure to use the same key!
application.secret="%APPLICATION_SECRET%"
# The application languages
application.langs="en"
# Global object class
# Define the Global object class for this application.
# Default to Global in the root package.sb
# application.global=Global
# Router
# Define the Router object to use for this application.
# This router will be looked up first when the application is starting up,
# so make sure this is the entry point.
# Furthermore, it's assumed your route file is named properly.
# So for an application router like `my.application.Router`,
# you may need to define a router file `conf/my.application.routes`.
# Default to Routes in the root package (and conf/routes)
# application.router=my.application.Routes
# Database configuration
# You can declare as many datasources as you want.
# By convention, the default datasource is named `default`
rootDir = "<path>/Bitcoin_price_prediction/"
db.default.driver = org.h2.Driver
db.default.url = "jdbc:h2: "<path>/Bitcoin_price_prediction/DataBase"
db.default.user = user
db.default.password = ""
play.evolutions.db.default.autoApply = true
# Evolutions
# You can disable evolutions if needed
# evolutionplugin=disabled
# Logger
# You can also configure logback (http://logback.qos.ch/),
# by providing an application-logger.xml file in the conf directory.
# Root logger:
logger.root=ERROR
# Logger used by the framework:
logger.play=INFO
# Logger provided to your application:
logger.application=DEBUG
#Enable JobModule to run scheduler
play.modules.enabled += "modules.jobs.JobModule"
#Frequency in seconds to run job. Might make sense to put 30 seconds, for recent data
constants.frequency = 30
ml.model_version = "gbt_22_binary_classes_32660767.model"
现在你可以理解，有几个变量需要根据你的平台和选择进行配置/更改：
- 将
rootDir目录更改为你在TrainGBT中使用的目录：
rootDir = "<path>/ Bitcoin_price_prediction"
- 指定数据库文件的名称：
db.default.url = "jdbc:h2: "<path>/Bitcoin_price_prediction/DataBase"
- 指定用于实际预测的模型版本：
ml.model_version = "gbt_22_binary_classes_32660767.model"
请注意，具有此名称的文件夹必须位于rootDir内部。因此，在rootDir中创建一个名为models的文件夹，并将所有训练模型的文件夹复制到其中。
该类还实现了actor特性并重写了receive方法。最佳实践是，在伴随对象中定义actor可以接收的类型，从而为其他类建立接口：
object PredictionActor {
def props = Props[PredictionActor]
case class PriceData(timeFrom: Long,
timeTo: Long,
priceDelta: (Long, Double)*)
}
首先，PredictionActor从models文件夹加载模型列表并加载etalon模型：
val models: List[(Transformer, String)] =
SubDirectoryRetriever.getListOfSubDirectories(modelFolder)
.map(modelMap => (PipelineModel.load(modelMap("path")),modelMap("modelName")))
.toList
首先，我们从models文件夹中提取子目录列表，并从每个子目录中加载训练好的PipeLine模型。以类似的方式加载etalon模型，但我们已经知道它的目录。以下是如何在receive方法中处理PriceData类型消息的方式：
override def receive: Receive = {
case data: PriceData =>
val priceData = shrinkData(data, 1, 22)
val (predictedLabelForUnknownTimestamp, details) =
predictionService.predictPriceDeltaLabel(priceData,productionModel)
预测的标签（字符串）和分类详细信息会被记录，是否有可能看到每个类别的概率分布？如果actor收到其他类型的消息，则会显示错误并且不执行任何操作。然后，结果被发送回SchedulerActor，并存储在predictedWithCurrent变量中，如前面的代码所示：
sender() ! CurrentDataWithShortTermPrediction(predictedLabelForUnknownTimestamp, data)
sender是一个ActorRef引用，指向正在处理的消息的发送者对象，因此我们可以使用!运算符将消息返回给它。然后，对于我们最开始加载的每个模型，我们会预测 1 分钟前的数据（总共 23 行中的第 0-21 行），并获取我们所知道的最新一分钟的实际价格变化：
models.foreach { mlModel =>
val (predictedLabel, details) =predictionService.predictPriceDeltaLabel(shrinkData(data, 0, 21), mlModel._1)
val actualDeltaPoint = data.priceDelta.toList(22)
对于每个模型，我们在数据库中存储以下信息：模型名称、每次测试预测的时间戳、模型预测的标签以及实际的价格变化（delta）。这些信息稍后会用于生成关于模型表现的报告：
storeShortTermBinaryPredictionIntoDB( mlModel._2, actualDeltaPoint._1,
predictedLabel, actualDeltaPoint._2)
TraderActor
TraderActor接收预测结果，并根据标签写入日志信息。它可以触发对指定端点的 HTTP 请求：
override def receive: Receive = {
case data: CurrentDataWithShortTermPrediction =>
Logger.debug("received short-term prediction" + data)
data.prediction match {
case "0" => notifySellShortTerm()
case "1" => notifyHoldShortTerm()
}
预测价格并评估模型
ShortTermPredictionServiceImpl是实际执行预测的类，使用给定的模型和数据。首先，它通过调用transformPriceData(priceData: PriceData)方法将PriceData转换为 Spark DataFrame，该 DataFrame 的模式与用于训练的数据模式相对应。接下来，调用model.transform(dataframe)方法；我们提取需要的变量，写入调试日志并返回给调用者：
override def predictPriceDeltaLabel(priceData: PriceData, mlModel: org.apache.spark.ml.Transformer): (String, Row) = {
val df = transformPriceData(priceData)
val prediction = mlModel.transform(df)
val predictionData = prediction.select("probability", "prediction", "rawPrediction").head()
(predictionData.get(1).asInstanceOf[Double].toInt.toString, predictionData)
}
在运行时，应用程序会收集有关预测输出的数据：预测标签和实际的价格变化（delta）。这些信息用于构建根网页，显示如TPR（真正率）、FPR（假正率）、TNR（真负率）和FNR（假负率）等统计数据，这些内容在之前已经描述过。
这些统计数据是从SHORT_TERM_PREDICTION_BINARY表中动态计算得出的。基本上，通过使用CASE-WHEN构造，我们添加了新的列：TPR、FPR、TNR 和 FNR。它们的定义如下：
-
如果预测标签为 1 且价格变化（delta）>
0，则 TPR 值为1，否则为0 -
如果预测标签为 1 且价格变化（delta）<=
0，则 FPR 值为1，否则为0 -
如果预测标签为 0 且价格变化（delta）<=
0，则 TNR 值为1，否则为0 -
如果预测标签为 0 且价格变化（delta）>
0，则 FNR 值为1，否则为0
然后，所有记录按照模型名称进行分组，TPR、FPR、TNR 和 FNR 进行求和，得到每个模型的总数。以下是负责此操作的 SQL 代码：
SELECT MODEL, SUM(TPR) as TPR, SUM(FPR) as FPR, SUM(TNR) as TNR,
SUM(FNR) as FNR, COUNT(*) as TOTAL FROM (SELECT *,
case when PREDICTED_LABEL='1' and ACTUAL_PRICE_DELTA > 0
then 1 else 0 end as TPR,
case when PREDICTED_LABEL='1' and ACTUAL_PRICE_DELTA <=0
then 1 else 0 end as FPR,
case when PREDICTED_LABEL='0' and ACTUAL_PRICE_DELTA <=0
then 1 else 0 end as TNR,
case when PREDICTED_LABEL='0' and ACTUAL_PRICE_DELTA > 0
then 1 else 0 end as FNR
FROM SHORT_TERM_PREDICTION_BINARY)
GROUP BY MODEL
使用 Scala Play 框架进行演示预测
现在我们已经看到了这个项目的所有步骤，是时候查看现场演示了。我们将把整个应用程序打包成一个 Scala Play 的 Web 应用程序。在查看演示之前，先让我们将项目启动并运行。然而，了解一些关于使用 Scala Play 的 RESTful 架构基础会非常有帮助。
为什么使用 RESTful 架构？
好吧，Play 的架构默认是 RESTful 的。其核心基于模型-视图-控制器（MVC）模式。每个入口点，配合 HTTP 动词，都会映射到一个控制器函数。控制器使得视图可以是网页、JSON、XML 或几乎任何其他形式。
Play 的无状态架构使得水平扩展成为可能，非常适合处理大量的请求，而不需要在它们之间共享资源（如会话）。它走在响应式编程（Reactive Programming）趋势的前沿，在这种编程模式下，服务器是事件驱动的，并通过并行处理来应对现代网站日益增长的需求。
在某些配置下，Play 启用完全异步和非阻塞的 I/O，贯穿整个应用程序。其目的是通过高效的线程管理和并行处理，在网络上实现更高的可扩展性，同时避免 JavaScript 解决方案常见的回调地狱问题。
AngularJS 是一个基于 JavaScript 的开源前端 Web 应用框架，主要由 Google 和一群个人及公司维护，旨在解决开发单页应用时遇到的诸多挑战。
现在问题是，为什么选择 AngularJS 呢？好吧，HTML 非常适合声明静态文档，但当我们尝试使用它来声明动态视图时，效果就不太理想。AngularJS 让您可以扩展 HTML 词汇以适应您的应用程序。最终的开发环境非常富有表现力、易于阅读，并且开发迅速。
另一个问题是，难道没有其他替代方案吗？好吧，其他框架通过抽象掉 HTML、CSS 和/或 JavaScript，或者通过提供一种命令式的方式来操作 DOM，解决了 HTML 的不足。但这些都没有解决根本问题，即 HTML 并未为动态视图设计。
最后，关于可扩展性如何呢？好吧，AngularJS 是一个用于构建最适合您应用程序开发框架的工具集。它是完全可扩展的，并且与其他库兼容。每个功能都可以被修改或替换，以适应您独特的开发工作流和功能需求。继续阅读，了解更多信息。
项目结构
包装后的 Scala Web ML 应用具有以下目录结构：
图 7：Scala ML Web 应用目录结构
在上述结构中，bitcoin_ml 文件夹包含所有后端和前端代码。models 文件夹包含所有训练好的模型。一个示例训练模型位于 gbt_22_binary_classes_32660767 文件夹中。最后，数据库文件和痕迹分别存放在 DataBase.mv.db 和 DataBase.trace.db 文件中。
那么我们来看看包含实际代码的 bitcoin_ml 文件夹的子文件夹结构：
图 8： bitcoin_ml 目录结构
在前面的图示中，conf 文件夹包含 Scala Web 应用配置文件 application.conf，该文件包含必要的配置（如已显示）。所有依赖项在 build.sbt 文件中定义，如下所示：
libraryDependencies ++= Seq(jdbc, evolutions,
"com.typesafe.play" %% "anorm" % "2.5.1",
cache, ws, specs2 % Test, ws)
unmanagedResourceDirectories in Test <+= baseDirectory(_ / "target/web/public/test")
resolvers += "scalaz-bintray" at "https://dl.bintray.com/scalaz/releases"
resolvers ++= Seq(
"apache-snapshots" at "http://repository.apache.org/snapshots/")
routesGenerator := InjectedRoutesGenerator
val sparkVersion = "2.2.0"
libraryDependencies += "org.apache.spark" %% "spark-mllib" % sparkVersion
libraryDependencies += "org.apache.hadoop" % "hadoop-mapreduce-client-core" % "2.7.2"
libraryDependencies += "org.apache.hadoop" % "hadoop-common" % "2.7.2"
libraryDependencies += "commons-io" % "commons-io" % "2.4"
libraryDependencies += "org.codehaus.janino" % "janino" % "3.0.7" //fixing "java.lang.ClassNotFoundException: de.unkrig.jdisasm.Disassembler" exception
libraryDependencies ++= Seq(
"com.typesafe.slick" %% "slick" % "3.1.1",
"org.slf4j" % "slf4j-nop" % "1.6.4"
)
坦白说，在最初写作时，我并没有想到将这个应用包装成 Scala Play Web 应用。因此，过程有些无序。然而，不必担心，想了解更多关于后端和前端的内容，请参考 第七章，使用 Q 学习和 Scala Play 框架进行期权交易。
运行 Scala Play Web 应用
要运行应用程序，只需按照以下步骤操作：
-
从
www.kaggle.com/mczielinski/bitcoin-historical-data下载历史比特币数据。然后解压并提取.csv文件。 -
打开你喜欢的 IDE（例如，Eclipse/IntelliJ）并创建 Maven 或 SBT 项目。
-
运行
Preprocess.scala脚本将历史数据转换为时间序列。此脚本应生成两个.csv文件（即scala_test_x.csv和scala_test_y.csv）。 -
然后使用之前生成的文件，训练
GradientBoostedTree模型（使用TrainGBT.scala脚本）。 -
保存最好的（即交叉验证的）
Pipeline模型，该模型包含所有管道的步骤。 -
然后从 Packt 仓库或 GitHub 下载 Scala Play 应用程序和所有文件（即
Bitcoin_price_prediction），可在书中找到（见书中内容）。 -
接下来将训练好的模型复制到
Bitcoin_price_prediction/models/目录下。 -
然后：
$ cd Bitcoin_price_prediction/bitcoin_ml/conf/并更新application.conf文件中的参数值，如前所示。 -
最后，使用
$ sudo sbt run命令运行该项目。
启动应用后，使用 $ sudo sbt run，应用程序将从 models 文件夹读取所有模型，etalon 模型由 ml.model_version 指定。每 30 秒（在 application.conf 文件中指定为 constants.frequency = 30），从 Cryptocompare API 获取最新的价格数据。使用 etalon 模型进行预测，结果会通过控制台日志消息的形式展示给用户，并有可能触发指定端点的 HTTP 请求。
之后，models 文件夹中的所有模型都会用于基于前 22 分钟的数据进行预测，并使用最新的价格数据作为当前分钟的预测质量检查方式。每个模型所做的所有预测都会存储在一个数据库文件中。当用户访问 http://localhost:9000 时，系统会展示一个包含预测摘要的表格给用户：
-
模型名称
-
TPR（实际上不是比率，在这种情况下只是原始计数）- 模型预测价格会增加的次数，以及这些预测中有多少次是正确的
-
FPR，模型预测价格会增加，但价格却下降或保持不变的次数
-
TNR，模型预测价格不增加并且正确的次数
-
FNR，即模型预测价格不升高时的错误次数
-
模型进行预测的总次数
好的，接下来是，使用$ sudo sbt run（在终端中）启动应用程序：
图 9：基于历史价格和实时数据生成的模型样本信号
前面的图显示了我们基于历史价格和实时数据生成的模型样本信号。此外，我们还可以看到模型的原始预测。当你尝试通过浏览器访问http://localhost:9000时，你应该能看到这个（不过，随着时间推移，计数会增加）：
图 10：使用 Scala Play2 框架的模型表现
在前面的图中，模型的表现并不令人满意，但我建议您使用最合适的超参数并进行更多次的训练，例如，训练 10,000 次。此外，在接下来的章节中，我尝试提供一些更深入的见解和改进指南。
最后，如果你计划在进行了一些扩展（如果有的话）后部署此应用程序，那么我建议你快速查看第七章中的最后一节，使用 Q-Learning 和 Scala Play 框架进行期权交易，在该节中，你将找到关于如何将应用程序部署为 Web 应用的服务器指南。
总结
本章实现了一个完整的机器学习流水线，从收集历史数据，到将其转换为适合测试假设的格式，训练机器学习模型，并在实时数据上运行预测，还可以评估多种不同的模型并选择最佳模型。
测试结果表明，与原始数据集类似，约 60 万个分钟（2.4 百万中的一部分）可以被分类为价格上涨（收盘价高于开盘价）；该数据集可以视为不平衡数据集。尽管随机森林通常在不平衡数据集上表现良好，但 ROC 曲线下的面积为 0.74，仍然不是最佳值。由于我们需要减少假阳性（减少在触发购买时价格下跌的次数），我们可能需要考虑一种更加严格的模型来惩罚这种错误。
尽管分类器获得的结果不能用于盈利交易，但它为测试新方法提供了基础，可以相对快速地进行测试。在这里，列出了一些进一步发展的可能方向：
-
实现一开始讨论的管道：将时间序列数据转换成几个簇，并为每个簇训练回归模型/分类器；然后将最近的数据分类到某个簇中，并使用为该簇训练的预测模型。根据定义，机器学习是从数据中推导模式，因此可能不会有一种适合比特币历史所有数据的模式；这就是为什么我们需要理解市场可能处于不同的阶段，而每个阶段有其独特的模式。
-
比特币价格预测的一个主要挑战可能是，训练数据（历史数据）在随机划分为训练集和测试集时，可能与测试数据的分布不一致。由于 2013 年和 2016 年间价格模式发生了变化，它们可能属于完全不同的分布。这可能需要对数据进行人工检查，并制作一些信息图表。可能已经有人做过这项研究。
-
其中一个主要的尝试方向是训练两个一对多分类器：一个用于预测价格上涨超过 20 美元，另一个用于预测价格下跌超过 20 美元；这样就可以相应地采取做多/做空的策略。
-
也许，预测下一分钟的变化并不是我们需要的；我们更愿意预测平均价格。因为开盘价可能比上一分钟的收盘价高很多，而下一分钟的收盘价可能稍微低于开盘价，但仍高于当前价格，这样的交易就会是有利可图的。所以，如何精确标记数据也是一个悬而未决的问题。
-
尝试使用不同的时间序列窗口大小（甚至 50 分钟可能适合），并使用 ARIMA 时间序列预测模型，因为它是最广泛使用的算法之一。然后，尝试预测价格变化，不是预测下一分钟的价格，而是预测接下来的 2-3 分钟的价格。此外，还可以尝试加入交易量。
-
如果价格在接下来的三个分钟中的至少一个时间段内上涨了 20 美元，标记数据为价格上涨，这样我们就可以从交易中获利。
-
目前，
Scheduler与 Cryptocompare 的分钟数据不同步。这意味着我们可以在下一个分钟的任何时刻获取关于 12:00:00 - 12:00:59 的分钟间隔数据，或者是 12:01:00 到 12:01:59 的数据。在后一种情况下，进行交易是没有意义的，因为我们是基于已经过时的数据做出的预测。 -
与其每分钟基于旧的数据进行预测来积累
actor的预测结果，不如直接获取最大可用的 HistoMinute 数据（七天），然后使用用于历史数据的 Scala 脚本将其拆分成时间序列数据，并对七天的数据进行预测。将这个过程作为定时任务每天运行一次，这样可以减少数据库和PredictionActor的负载。 -
与通常的数据集相比，在比特币中，历史数据的行是按日期升序排列的，这意味着：
-
最新的数据可能与今天的价格更相关，而“少即是多”的方法也适用；采用数据的一个较小子集可能会带来更好的性能。
-
数据的子抽样方法可能会影响结果（如将数据分为训练集和测试集）。
-
最后，尝试使用 LSTM 网络以获得更好的预测准确性（可以参考第十章获得一些提示）。
-
对基因组序列中变异的理解有助于我们识别容易罹患常见疾病的个体，解决罕见疾病问题，并从一个更大的群体中找到相应的特定人群。尽管经典的机器学习技术可以帮助研究人员识别相关变量的群体（簇），但这些方法在处理像全人类基因组这样的大型和高维数据集时，其准确性和有效性会有所下降。另一方面，深度神经网络架构（深度学习的核心）能更好地利用大规模数据集来构建复杂的模型。
在下一章中，我们将看到如何在《千人基因组计划》的大规模基因组数据上应用 K-means 算法，旨在对人群规模上的基因型变异进行聚类。然后，我们将训练一个基于 H2O 的深度学习模型，用于预测地理民族。最后，将使用基于 Spark 的随机森林算法来提高预测准确性。
第四章：人口规模的聚类与民族预测
理解基因组序列的变异有助于我们识别易患常见疾病的人群、治愈罕见疾病，并从更大的群体中找到对应的目标群体。尽管经典的机器学习技术可以帮助研究人员识别相关变量的群体（即簇），但这些方法在处理如整个基因组这类大规模和高维度数据集时，准确性和有效性会下降。
另一方面，深度神经网络 (DNNs) 形成了深度学习 (DL)的核心，并提供了建模复杂、高层次数据抽象的算法。它们能够更好地利用大规模数据集来构建复杂的模型。
本章中，我们将应用 K-means 算法对来自 1000 基因组计划分析的大规模基因组数据进行聚类，旨在在人群规模上聚类基因型变异。最后，我们训练一个基于 H2O 的 DNN 模型和一个基于 Spark 的随机森林模型，用于预测地理民族。本章的主题是给我你的基因变异数据，我会告诉你你的民族。
尽管如此，我们将配置 H2O，以便在接下来的章节中也可以使用相同的设置。简而言之，我们将在这个端到端项目中学习以下内容：
-
人口规模的聚类与地理民族预测
-
1000 基因组计划——一个深入的人的基因变异目录
-
算法与工具
-
使用 K-means 进行人口规模的聚类
-
使用 H2O 进行民族预测
-
使用随机森林进行民族预测
人口规模聚类与地理民族
下一代基因组测序 (NGS)减少了基因组测序的开销和时间，以前所未有的方式产生了大数据。相比之下，分析这些大规模数据在计算上非常昂贵，且逐渐成为关键瓶颈。随着 NGS 数据中样本数量和每个样本特征的增加，这对大规模并行数据处理提出了需求，从而对机器学习解决方案和生物信息学方法带来了前所未有的挑战。在医疗实践中使用基因组信息需要高效的分析方法来应对来自数千人及其数百万变异的数据。
其中一项最重要的任务是分析基因组特征，以便将个体归类为特定的民族群体，或分析核苷酸单倍型与疾病易感性的关系。来自 1000 基因组计划的数据作为分析全基因组单核苷酸多态性 (SNPs)的主要来源，旨在预测个体的祖先背景，涉及大陆和区域的起源。
基因变异的机器学习
研究表明，来自亚洲、欧洲、非洲和美洲的群体可以根据其基因组数据进行区分。然而，准确预测单倍群和起源大陆（即，地理、民族和语言）的难度更大。其他研究表明，Y 染色体谱系可以在地理上定位，为（地理）聚类人类基因型中的人类等位基因提供证据。
因此，个体的聚类与地理来源和祖先有关系。由于种族也依赖于祖先，聚类与更传统的种族概念之间也存在关联，但这种关联并不完全准确，因为基因变异是按照概率原则发生的。因此，它并不遵循在不同种族间的连续分布，而是呈现出交叉或重叠的现象。
因此，确定祖先，甚至是种族，可能对生物医学有一定的用处，但任何直接评估与疾病相关的基因变异，最终将提供更为准确和有益的信息。
各种基因组学项目提供的数据集，如癌症基因组图谱（TCGA）、国际癌症基因组联盟（ICGC）、1000 基因组计划以及个人基因组计划（PGP），都包含了大规模数据。为了快速处理这些数据，已提出基于 ADAM 和 Spark 的解决方案，并且这些解决方案现在广泛应用于基因组数据分析研究。
Spark 形成了最有效的数据处理框架，并且提供了内存集群计算的基本构件，例如用于反复查询用户数据的功能。这使得 Spark 成为机器学习算法的优秀候选框架，其性能超过了基于 Hadoop 的 MapReduce 框架。通过使用来自 1000 基因组计划的基因变异数据集，我们将尝试回答以下问题：
-
人类的基因变异在不同群体之间的地理分布是怎样的？
-
我们能否利用个体的基因组信息，将其归类为特定的群体，或从其核苷酸单倍型中推导出疾病易感性？
-
个体的基因组数据是否适合预测其地理来源（即，个体所属的群体）？
在本项目中，我们以一种可扩展且更高效的方式解决了前述问题。特别地，我们研究了如何应用 Spark 和 ADAM 进行大规模数据处理，使用 H2O 对整个群体进行 K-means 聚类以确定群体内外的组别，以及通过调节更多超参数来进行基于 MLP 的监督学习，以更准确地根据个体的基因组数据预测该个体所属的群体。现在不必担心；我们将在后续部分提供关于这些技术的工作细节。
然而，在开始之前，让我们简要了解一下 1000 基因组项目的数据集，以便为您提供一些关于为什么跨技术互操作如此重要的理由。
1000 基因组项目数据集描述
1000 基因组项目的数据是一个非常庞大的人类遗传变异目录。该项目旨在确定在研究的人群中频率超过 1%的遗传变异。数据已经公开，并通过公共数据仓库向全球科学家自由访问。此外，1000 基因组项目的数据广泛用于筛选在遗传性疾病个体的外显子数据中发现的变异，以及在癌症基因组项目中的应用。
变异调用格式（VCF）中的基因型数据集提供了人类个体（即样本）及其遗传变异的数据，此外，还包括全球的等位基因频率，以及各超人群的等位基因频率。数据指明了每个样本所属的人群地区，这些信息用于我们方法中的预测类别。特定的染色体数据（以 VCF 格式呈现）可能包含其他信息，表明样本的超人群或所使用的测序平台。对于多等位基因变异，每个替代等位基因频率（AF）以逗号分隔的列表形式呈现，如下所示：
1 15211 rs78601809 T G 100 PASS AC=3050;
AF=0.609026;
AN=5008;
NS=2504;
DP=32245;
EAS_AF=0.504;
AMR_AF=0.6772;
AFR_AF=0.5371;
EUR_AF=0.7316;
SAS_AF=0.6401;
AA=t|||;
VT=SNP
等位基因频率（AF）是通过等位基因计数（AC）与等位基因总数（AN）的商计算得出的，而 NS 则是具有数据的样本总数，而_AF表示特定区域的 AF 值。
1000 基因组计划始于 2008 年；该联盟由 400 多名生命科学家组成，第三阶段于 2014 年 9 月完成，共涵盖了来自 26 个人群（即种族背景）的2,504个个体。总共识别出了超过 8800 万个变异（8470 万个单核苷酸多态性（SNPs）、360 万个短插入/缺失（indels）和 6 万个结构变异），这些变异被确定为高质量的单倍型。
简而言之，99.9%的变异是由 SNPs 和短插入/缺失（indels）组成的。较不重要的变异——包括 SNPs、indels、缺失、复杂的短替代以及其他结构变异类别——已被去除以进行质量控制。因此，第三阶段发布的数据保留了 8440 万个变异。
这 26 个人群中的每一个大约有 60 到 100 个来自欧洲、非洲、美洲（南美和北美）以及亚洲（南亚和东亚）的个体。这些人群样本根据其主要祖先的来源被分为超人群组：东亚人群（CHB，JPT，CHS，CDX，和KHV）、欧洲人群（CEU，TSI，FIN，GBR，和IBS）、非洲人群（YRI，LWK，GWD，MSL，ESN，ASW，和ACB）、美洲人群（MXL，PUR，CLM，和PEL）以及南亚人群（GIH，PJL，BEB，STU，和ITU）。具体内容请参考图 1：
图 1：来自 1000 基因组项目发布版 3 的地理种族群体（来源 www.internationalgenome.org/）
发布的数据集提供了来自 2,504 名健康成年人的数据（年龄 18 岁及以上，第三阶段项目）；只有至少 70 碱基对 (bp) 的读取被使用，直到有更先进的解决方案可用为止。所有来自所有样本的基因组数据被整合，以便将所有变异归因于某个区域。然而，值得注意的是，特定的单倍型可能不会出现在某个特定区域的基因组中；也就是说，多样本方法允许将变异归因于个体的基因型，即使这些变异未被该样本的测序读取覆盖。
换句话说，提供的是重叠的读取数据，并且单一样本基因组未必已被整合。所有个体均使用以下两种技术进行测序：
-
全基因组测序 (平均深度 = 7.4x，其中 x 表示在给定参考 bp 上，平均可能对齐的读取数量)
-
靶向外显子组测序 (平均深度 = 65.7x)
此外，个体及其直系亲属（如成人后代）通过高密度 SNP 微阵列进行基因分型。每个基因型包含所有 23 对染色体，并且一个单独的面板文件记录了样本和种群信息。表 1 给出了 1000 基因组项目不同发布版的概览：
表 1 – 1000 基因组项目基因型数据集统计 （来源： www.internationalgenome.org/data ）
| 1000 基因组发布版 | 变异 | 个体 | 种群 | 文件格式 |
|---|---|---|---|---|
| 第 3 阶段 | 第 3 阶段 | 2,504 | 26 | VCF |
| 第 1 阶段 | 3,790 万 | 1,092 | 14 | VCF |
| 试点 | 1,480 万 | 179 | 4 | VCF |
计算五个超级种群组中的等位基因频率（AF）：EAS=东亚，EUR=欧洲，AFR=非洲，AMR=美洲，SAS=南亚，这些频率来自等位基因数（AN，范围 = [0, 1]）。
请参阅面板文件的详细信息：ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel。
算法、工具和技术
来自 1000 基因组项目发布版 3 的大规模数据贡献了 820 GB 的数据。因此，使用 ADAM 和 Spark 来以可扩展的方式预处理和准备数据（即训练、测试和验证集），以供 MLP 和 K-means 模型使用。Sparkling water 用于在 H2O 和 Spark 之间转换数据。
然后，使用 K-means 聚类和多层感知机（MLP，使用 H2O）进行训练。对于聚类和分类分析，需要使用样本 ID、变异 ID 以及替代等位基因的计数，这些基因型信息来自每个样本，我们所使用的大多数变异为 SNP 和插入缺失（indels）。
现在，我们应该了解每个工具的最基本信息，如 ADAM、H2O 以及关于算法的一些背景信息，如 K-means、MLP 用于聚类和分类人口群体。
H2O 和 Sparkling water
H2O 是一个机器学习的人工智能平台。它提供了丰富的机器学习算法和一个基于网页的数据处理用户界面，既有开源版本，也有商业版本。使用 H2O，可以通过多种语言（如 Java、Scala、Python 和 R）开发机器学习和深度学习应用：
图 2：H2O 计算引擎及其可用功能（来源：https://h20.ai/）
它还具有与 Spark、HDFS、SQL 和 NoSQL 数据库接口的能力。简而言之，H2O 可在 Hadoop/Yarn、Spark 或笔记本电脑上与 R、Python 和 Scala 配合使用。另一方面，Sparkling water 将 H2O 的快速、可扩展机器学习算法与 Spark 的功能结合起来。它通过 Scala/R/Python 驱动计算，并利用 H2O flow 用户界面。简而言之，Sparkling water = H2O + Spark。
在接下来的几章中，我们将探索 H2O 和 Sparkling water 的广泛丰富功能；然而，我认为提供一张涵盖所有功能领域的图表会更有帮助：
图 3：可用算法和支持的 ETL 技术概览（来源：https://h20.ai/）
这是从 H2O 网站整理出的功能和技术列表。它可以用于数据清洗、使用数据建模以及评分结果模型：
-
流程
-
模型
-
评分工具
-
数据概况
-
广义线性模型（GLM）
-
预测
-
总结统计
-
决策树
-
混淆矩阵
-
聚合、过滤、分箱和推导列
-
梯度提升机（GBM）
-
AUC
-
切片、对数变换和匿名化
-
K-means
-
命中率
-
变量创建
-
异常检测
-
PCA/PCA 评分
-
深度学习（DL）
-
多模型评分
-
训练和验证抽样计划
-
朴素贝叶斯
-
网格搜索
以下图表展示了如何清晰地描述 H2O Sparkling water 如何被用来扩展 Apache Spark 的功能。H2O 和 Spark 都是开源系统。Spark MLlib 包含大量功能，而 H2O 在此基础上扩展了许多额外功能，包括深度学习（DL）。它提供了数据转换、建模和评分的工具，正如我们在 Spark ML 中看到的那样。它还提供了一个基于网页的用户界面来进行交互：
图 4：Sparkling water 扩展 H2O，并与 Spark 互操作（来源：https://h20.ai/）
以下图表展示了 H2O 如何与 Spark 集成。正如我们所知，Spark 有主节点和工作节点；工作节点创建执行器来执行实际的工作。以下是运行基于 Sparkling water 的应用程序的步骤：
-
Spark 的提交命令将 Sparkling water JAR 发送到 Spark 主节点
-
Spark 主节点启动工作节点并分发 JAR 文件
-
Spark 工作节点启动 Executor JVM 以执行工作
-
Spark 执行器启动一个 H2O 实例
H2O 实例嵌入在 Executor JVM 中，因此它与 Spark 共享 JVM 堆内存。当所有 H2O 实例启动后，H2O 将形成一个集群，并提供 H2O flow 网页界面：
图 5：Sparkling Water 如何融入 Spark 架构（来源：http://blog.cloudera.com/blog/2015/10/how-to-build-a-machine-learning-app-using-sparkling-water-and-apache-spark/）
前面的图解释了 H2O 如何融入 Spark 架构以及如何启动，但数据共享又该如何处理呢？现在的问题是：数据如何在 Spark 和 H2O 之间传递？下面的图解释了这一点：
图 6：Spark 和 H2O 之间的数据传递机制
为了更清晰地查看前面的图，已经为 H2O 和 Sparkling Water 创建了一个新的 H2O RDD 数据结构。它是基于 H2O 框架顶部的一个层，每一列代表一个数据项，并独立压缩，以提供最佳的压缩比。
ADAM 用于大规模基因组数据处理
分析 DNA 和 RNA 测序数据需要大规模的数据处理，以根据其上下文解释数据。优秀的工具和解决方案已经在学术实验室中开发出来，但往往在可扩展性和互操作性方面存在不足。因此，ADAM 是一个基因组分析平台，采用了 Apache Avro、Apache Spark 和 Parquet 构建的专用文件格式。
然而，像 ADAM-Spark 这样的规模化数据处理解决方案可以直接应用于测序管道的输出数据，也就是说，在质量控制、比对、读段预处理和变异定量之后，使用单一样本数据进行处理。举几个例子，DNA 测序的 DNA 变异、RNA 测序的读数计数等。
更多内容见bdgenomics.org/以及相关的出版物：Massie, Matt 和 Nothaft, Frank 等人，《ADAM：用于云规模计算的基因组格式与处理模式》，UCB/EECS-2013-207，加利福尼亚大学伯克利分校电子工程与计算机科学系。
在我们的研究中，使用 ADAM 实现了支持 VCF 文件格式的可扩展基因组数据分析平台，从而将基于基因型的 RDD 转化为 Spark DataFrame。
无监督机器学习
无监督学习是一种机器学习算法，用于通过推理未标记数据集（即由没有标签的输入数据组成的训练集）来对相关数据对象进行分组并发现隐藏的模式。
让我们看一个现实生活中的例子。假设你在硬盘上有一个装满非盗版且完全合法的 MP3 文件的大文件夹。如果你能建立一个预测模型，帮助你自动将相似的歌曲分组并组织成你喜欢的类别，比如乡村、说唱和摇滚，如何？
这是一种将项目分配到组中的行为，以便 MP3 文件能够以无监督的方式添加到相应的播放列表中。对于分类，我们假设你提供了一个正确标记的数据训练集。不幸的是，当我们在现实世界中收集数据时，我们并不总是有这种奢侈的条件。
例如，假设我们想将大量的音乐分成有趣的播放列表。如果我们没有直接访问它们的元数据，我们如何可能将歌曲分组呢？一种可能的方法是混合使用各种机器学习技术，但聚类通常是解决方案的核心：
图 7：聚类数据样本一览
换句话说，无监督学习算法的主要目标是探索输入数据中未知/隐藏的模式，这些数据是没有标签的。然而，无监督学习也包含其他技术，以探索性方式解释数据的关键特征，寻找隐藏的模式。为了克服这一挑战，聚类技术被广泛应用于基于某些相似性度量，无监督地对未标记数据点进行分组。
人群基因组学与聚类
聚类分析是将数据样本或数据点划分并放入相应的同质类或聚类中。因此，聚类的简单定义可以被看作是将对象组织成一些成员在某种方式上相似的组的过程，如下所示。
这样，聚类就是将一些在某种程度上彼此相似的对象集合在一起，并与属于其他聚类的对象不相似的集合。如果给定的是遗传变异集合，聚类算法会基于相似性将这些对象放入一个组中——也就是人口组或超级人口组。
K 均值算法是如何工作的？
聚类算法，如 K 均值算法，定位数据点组的质心。然而，为了使聚类更加准确和有效，算法会评估每个点与聚类质心的距离。
最终，聚类的目标是确定一组未标记数据中的内在分组。例如，K 均值算法尝试将相关的数据点聚类到预定义的三个（即k = 3）聚类中，如图 8所示：
图 8：典型聚类算法的结果及聚类中心的表示
在我们的案例中，结合使用 Spark、ADAM 和 H2O 的方法能够处理大量不同的数据点。假设我们有 n 个数据点（x[i]，i=1, 2… n，举例来说，是遗传变异），这些数据点需要被分为k个簇。然后，K-means 将每个数据点分配给一个簇，目标是找到簇的位置μ[i]，i=1…k，以最小化数据点到簇的距离。从数学上讲，K-means 试图通过解一个方程来实现这一目标——也就是一个优化问题：
在前面的方程中，c[i]是分配给簇i的数据点集合，d(x,μ[i])=∥x−μ[i]∥[2]²是需要计算的欧几里得距离。该算法通过最小化簇内平方和（即WCSS），计算数据点与 k 个簇中心之间的距离，其中c[i]是属于簇i的点的集合。
因此，我们可以理解，使用 K-means 进行的总体聚类操作并非一项简单的任务，而是一个 NP 难优化问题。这也意味着 K-means 算法不仅尝试找到全局最小值，而且经常陷入不同的解中。K-means 算法通过交替执行两个步骤进行：
-
簇分配步骤：将每个观测值分配给使得均值产生最小WCSS的簇。平方和是平方欧几里得距离。
-
质心更新步骤：计算新簇中观测值的均值作为质心。
简而言之，K-means 训练的整体方法可以用下图描述：
图 9：K-means 算法过程的总体方法
用于地理族群预测的 DNN
多层感知器（MLP）是一个示例 DNN，它是一个前馈神经网络；即神经元之间只有不同层之间的连接。它有一个（通过）输入层，一个或多个线性阈值单元（LTUs）（称为隐藏层），以及一个 LTU 的最终层（称为输出层）。
每一层（输出层除外）都涉及一个偏置神经元，并且与下一层完全连接，形成一个完全连接的二分图。信号只从输入流向输出，即单向(前馈)
直到最近，MLP 是通过反向传播训练算法进行训练的，但现在优化版本（即梯度下降）使用反向模式自动微分；即神经网络通过 SGD 和反向传播作为梯度计算技术进行训练。在 DNN 训练中，为了解决分类问题，使用了两层抽象：
-
梯度计算：使用反向传播
-
优化层级：使用 SGD、ADAM、RMSPro 和 Momentum 优化器来计算之前计算的梯度
在每个训练周期中，算法将数据输入到网络中，并计算每个神经元在连续层中的状态和输出。然后，算法衡量网络中的输出误差，即期望输出与当前输出之间的差距，并计算最后一个隐藏层中每个神经元对该误差的贡献。
逐步地，输出误差通过所有隐藏层反向传播到输入层，并在反向传播过程中计算所有连接权重的误差梯度：
图 10：由输入层、ReLU 和 softmax 组成的现代 MLP
对于多类别分类任务，输出层通常由一个共享的 softmax 函数决定（更多内容请参见图 2），与单独的激活函数不同，每个输出神经元提供对应类别的估计概率。
此外，我们将使用树集成方法，例如随机森林来进行分类。目前，我认为我们可以跳过随机森林的基础介绍，因为我们已经在第一章《分析保险赔付严重性》、第二章《分析和预测电信流失》和第三章《基于历史数据的高频比特币价格预测》中详细讲解过它。好了，现在是时候开始了。不过，在动手之前，确保你的编程环境已经准备好总是一个好习惯。
配置编程环境
在本节中，我们将介绍如何配置我们的编程环境，以便能够与 Spark、H2O 和 ADAM 互操作。请注意，在笔记本电脑或台式机上使用 H2O 会消耗大量资源。因此，请确保你的笔记本至少有 16GB 的 RAM 和足够的存储空间。
无论如何，我打算在 Eclipse 上将这个项目设置为 Maven 项目。不过，你也可以尝试在 SBT 中定义相同的依赖项。让我们在 pom.xml 文件中定义属性标签，以便构建一个适合 Maven 的项目：
<properties>
<spark.version>2.2.1</spark.version>
<scala.version>2.11.12</scala.version>
<h2o.version>3.16.0.2</h2o.version>
<sparklingwater.version>2.2.6</sparklingwater.version>
<adam.version>0.23.0</adam.version>
</properties>
然后，我们可以使用 Spark 2.2.1 的最新版本（任何 2.x 版本甚至更高版本都应该能正常工作）：
<dependency>
<groupId>org.apache.spark</groupId>
<artifactId>spark-core_2.11</artifactId>
<version>${spark.version}</version>
</dependency>
然后，我们需要声明 H2O 和 Sparkling Water 的依赖项，确保它们与属性标签中指定的版本匹配。较新版本可能也能正常工作，你可以尝试：
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>sparkling-water-core_2.11</artifactId>
<version>2.2.6</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>sparkling-water-examples_2.11</artifactId>
<version>2.2.6</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>h2o-core</artifactId>
<version>${h2o.version}</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>h2o-scala_2.11</artifactId>
<version>${h2o.version}</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>h2o-algos</artifactId>
<version>${h2o.version}</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>h2o-app</artifactId>
<version>${h2o.version}</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>h2o-persist-hdfs</artifactId>
<version>${h2o.version}</version>
</dependency>
<dependency>
<groupId>ai.h2o</groupId>
<artifactId>google-analytics-java</artifactId>
<version>1.1.2-H2O-CUSTOM</version>
</dependency>
最后，让我们定义 ADAM 及其依赖项：
<dependency>
<groupId>org.bdgenomics.adam</groupId>
<artifactId>adam-core_2.11</artifactId>
<version>0.23.0</version>
</dependency>
当我在 Windows 机器上尝试时，我还需要安装 joda-time 依赖项。让我们来做这件事（但根据你的平台，可能不需要）：
<dependency>
<groupId>joda-time</groupId>
<artifactId>joda-time</artifactId>
<version>2.9.9</version>
</dependency>
一旦你在 Eclipse 中创建了一个 Maven 项目（无论是通过 IDE 手动创建还是使用 $ mvn install），所有必需的依赖项将会被下载！我们现在可以开始编码了！
等等！如何在浏览器中查看 H2O 的 UI 界面？为此，我们必须手动下载 H2O 的 JAR 文件并将其作为常规的.jar文件运行。简而言之，这是一个三步过程：
-
从
www.h2o.ai/download/下载最新稳定版H[2]O。然后解压，它包含了开始所需的一切。 -
从你的终端/命令提示符中，使用
java -jar h2o.jar运行.jar文件。 -
在浏览器中输入
http://localhost:54321：
图 11：H2O FLOW 的 UI 界面
这显示了最新版本（即截至 2018 年 1 月 19 日的 h2o-3.16.0.4 版）H2O 的可用功能。然而，我不打算在这里解释所有内容，所以让我们停止探索，因为我相信目前为止，这些关于 H2O 和 Sparkling Water 的知识已经足够。
数据预处理和特征工程
我已经说明所有的 24 个 VCF 文件贡献了 820 GB 的数据。因此，我决定只使用 Y 染色体的遗传变异来使演示更加清晰。其大小约为 160 MB，这样不会带来巨大的计算挑战。你可以从 ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/下载所有的 VCF 文件以及面板文件。
让我们开始吧。我们从创建SparkSession开始，这是 Spark 应用程序的网关：
val spark:SparkSession = SparkSession
.builder()
.appName("PopStrat")
.master("local[*]")
.config("spark.sql.warehouse.dir", "C:/Exp/")
.getOrCreate()
然后让我们告诉 Spark VCF 文件和面板文件的路径：
val genotypeFile = "<path>/ALL.chrY.phase3_integrated_v2a.20130502.genotypes.vcf"
val panelFile = "<path>/integrated_call_samples_v3.20130502.ALL.panel "
我们使用 Spark 处理面板文件，以访问目标种群数据并识别种群组。我们首先创建一个我们想要预测的种群集合：
val populations = Set("FIN", "GBR", "ASW", "CHB", "CLM")
然后我们需要创建一个样本 ID → 种群的映射，以便筛选出我们不感兴趣的样本：
def extract(file: String,
filter: (String, String) => Boolean): Map[String, String] = {
Source
.fromFile(file)
.getLines()
.map(line => {
val tokens = line.split(Array('t', ' ')).toList
tokens(0) -> tokens(1)
}).toMap.filter(tuple => filter(tuple._1, tuple._2))
}
val panel: Map[String, String] = extract(
panelFile,
(sampleID: String, pop: String) => populations.contains(pop))
请注意，面板文件会生成所有个体的样本 ID、种群组、族群、超种群组以及性别，具体如下所示：
图 12：样本面板文件的内容
然后加载 ADAM 基因型，并筛选基因型，使我们仅保留那些在我们感兴趣的种群中的基因型：
val allGenotypes: RDD[Genotype] = sc.loadGenotypes(genotypeFile).rdd
val genotypes: RDD[Genotype] = allGenotypes.filter(genotype => {
panel.contains(genotype.getSampleId)
})
接下来的工作是将Genotype对象转换为我们自己的SampleVariant对象，以尽量节省内存。然后，将genotype对象转换为一个SampleVariant对象，其中只包含我们进一步处理所需的数据：样本 ID（唯一标识特定样本）；变异 ID（唯一标识特定遗传变异）；以及替代等位基因的计数（仅在样本与参考基因组不同的情况下）。
准备样本变异的签名如下所示；它接受sampleID、variationId和alternateCount：
case class SampleVariant(sampleId: String,
variantId: Int,
alternateCount: Int)
好的！让我们从genotype文件中找到variantID。variantId是一个String类型，由名称、起始位置和染色体上的终止位置组成：
def variantId(genotype: Genotype): String = {
val name = genotype.getVariant.getContigName
val start = genotype.getVariant.getStart
val end = genotype.getVariant.getEnd
s"$name:$start:$end"
}
一旦我们有了variantID，就应该寻找替代计数。在 genotype 文件中，那些没有等位基因参考的对象大致上是遗传变体：
def alternateCount(genotype: Genotype): Int = {
genotype.getAlleles.asScala.count(_ != GenotypeAllele.REF)
}
最后，我们构建一个简单的变体对象。为此，我们需要对样本 ID 进行内联，因为它们在 VCF 文件中会频繁出现：
def toVariant(genotype: Genotype): SampleVariant = {
new SampleVariant(genotype.getSampleId.intern(),
variantId(genotype).hashCode(),
alternateCount(genotype))
}
很棒！我们已经能够构建简单的变体了。现在，下一项具有挑战性的任务是准备 variantsRDD，然后我们才能创建 variantsBySampleId RDD：
val variantsRDD: RDD[SampleVariant] = genotypes.map(toVariant)
然后，我们必须按样本 ID 对变体进行分组，以便逐个样本处理变体。之后，我们可以获得将用于查找缺失变体的样本总数。最后，我们必须按变体 ID 对变体进行分组，并过滤掉那些在某些样本中缺失的变体：
val variantsBySampleId: RDD[(String, Iterable[SampleVariant])] =
variantsRDD.groupBy(_.sampleId)
val sampleCount: Long = variantsBySampleId.count()
println("Found " + sampleCount + " samples")
val variantsByVariantId: RDD[(Int, Iterable[SampleVariant])] =
variantsRDD.groupBy(_.variantId).filter {
case (_, sampleVariants) => sampleVariants.size == sampleCount
}
现在，让我们创建一个变体 ID → 替代计数大于零的样本计数的映射。然后我们过滤掉那些不在所需频率范围内的变体。这里的目标只是减少数据集中的维度数量，以便更容易训练模型：
val variantFrequencies: collection.Map[Int, Int] = variantsByVariantId
.map {
case (variantId, sampleVariants) =>
(variantId, sampleVariants.count(_.alternateCount > 0))
}.collectAsMap()
在基于变体 ID 对样本进行分组并过滤掉没有支持的变体之前，样本（或个体）的总数已确定，以简化数据预处理并更好地应对大量变体（总计 8440 万个）。
图 13 显示了 1000 基因组项目中基因型变体集合的概念视图，并展示了从相同数据中提取特征，以训练我们的 K-means 和 MLP 模型的过程：
图 13：1000 基因组项目中基因型变体集合的概念视图
指定的范围是任意的，选择它是因为它包含了一个合理数量的变体，但不至于太多。具体来说，对于每个变体，都计算了替代等位基因的频率，排除了少于 12 个替代等位基因的变体，分析中留下了约 300 万个变体（来自 23 个染色体文件）：
val permittedRange = inclusive(11, 11)
val filteredVariantsBySampleId: RDD[(String, Iterable[SampleVariant])] =
variantsBySampleId.map {
case (sampleId, sampleVariants) =>
val filteredSampleVariants = sampleVariants.filter(
variant =>
permittedRange.contains(
variantFrequencies.getOrElse(variant.variantId, -1)))
(sampleId, filteredSampleVariants)
}
一旦我们有了filteredVariantsBySampleId，下一个任务是对每个样本 ID 的变体进行排序。每个样本现在应该具有相同数量的排序变体：
val sortedVariantsBySampleId: RDD[(String, Array[SampleVariant])] =
filteredVariantsBySampleId.map {
case (sampleId, variants) =>
(sampleId, variants.toArray.sortBy(_.variantId))
}
println(s"Sorted by Sample ID RDD: " + sortedVariantsBySampleId.first())
现在，RDD 中的所有项应该具有相同顺序的相同变体。最终任务是使用sortedVariantsBySampleId构建一个包含区域和替代计数的 Row 类型 RDD：
val rowRDD: RDD[Row] = sortedVariantsBySampleId.map {
case (sampleId, sortedVariants) =>
val region: Array[String] = Array(panel.getOrElse(sampleId, "Unknown"))
val alternateCounts: Array[Int] = sortedVariants.map(_.alternateCount)
Row.fromSeq(region ++ alternateCounts)
}
因此，我们只需使用第一个来构建训练数据框架的头部：
val header = StructType(
Seq(StructField("Region", StringType)) ++
sortedVariantsBySampleId
.first()
._2
.map(variant => {
StructField(variant.variantId.toString, IntegerType)
}))
干得好！到目前为止，我们已经得到了我们的 RDD 和 StructType 头信息。现在，我们可以通过最小的调整/转换来使用 H2O 和 Spark 的深度/机器学习算法。整个端到端项目的流程如下图所示：
图 14：整体方法的管道流程
模型训练与超参数调优
一旦我们有了rowRDD和表头，接下来的任务是使用表头和rowRDD构建 Schema DataFrame 的行：
val sqlContext = spark.sqlContext
val schemaDF = sqlContext.createDataFrame(rowRDD, header)
schemaDF.printSchema()
schemaDF.show(10)
>>>
图 15：包含特征和标签（即 Region）列的训练数据集快照
在前面的 DataFrame 中，仅显示了少数列，包括标签，以便适合页面。
基于 Spark 的 K-means 用于大规模人口聚类
在前面的章节中，我们已经了解了 K-means 的工作原理。所以我们可以直接进入实现部分。由于训练是无监督的，我们需要删除标签列（即Region）：
val sqlContext = sparkSession.sqlContext
val schemaDF = sqlContext.createDataFrame(rowRDD, header).drop("Region")
schemaDF.printSchema()
schemaDF.show(10)
>>>
图 16：没有标签（即 Region）的 K-means 训练数据集快照
现在，我们已经在第一章，《分析保险严重性索赔》和第二章，《分析与预测电信流失》中看到，Spark 期望监督训练时有两列（即特征和标签），而无监督训练时它只期望包含特征的单列。由于我们删除了标签列，现在需要将整个变量列合并为一个features列。因此，我们将再次使用VectorAssembler()转换器。首先，我们选择要嵌入向量空间的列：
val featureCols = schemaDF.columns
然后，我们实例化VectorAssembler()转换器，指定输入列和输出列：
val assembler =
new VectorAssembler()
.setInputCols(featureCols)
.setOutputCol("features")
val assembleDF = assembler.transform(schemaDF).select("features")
现在让我们看看它的效果：
assembleDF.show()
>>>
图 17：K-means 特征向量的快照
由于我们的数据集维度非常高，我们可以使用一些降维算法，如 PCA。因此，我们通过实例化一个PCA()转换器来实现，如下所示：
val pca =
new PCA()
.setInputCol("features")
.setOutputCol("pcaFeatures")
.setK(50)
.fit(assembleDF)
然后我们转换组合后的 DataFrame（即已组合的）和前 50 个主成分。你可以调整这个数量。最后，为了避免歧义，我们将pcaFeatures列重命名为features：
val pcaDF = pca.transform(assembleDF)
.select("pcaFeatures")
.withColumnRenamed("pcaFeatures", "features")
pcaDF.show()
>>>
图 18：前 50 个主成分的快照，作为最重要的特征
很棒！一切顺利。最后，我们准备好训练 K-means 算法了：
val kmeans =
new KMeans().setK(5).setSeed(12345L)
val model = kmeans.fit(pcaDF)
所以我们通过计算组内平方误差和（WSSSE）来评估聚类效果：
val WSSSE = model.computeCost(pcaDF)
println("Within-Cluster Sum of Squares for k = 5 is" + WSSSE)
>>>
确定最佳聚类数
像 K-means 这样的聚类算法的优点是，它们可以对具有无限特征数的数据进行聚类。当你拥有原始数据并且希望了解数据中的模式时，它们是非常有用的工具。然而，在进行实验之前决定聚类数可能并不成功，有时甚至会导致过拟合或欠拟合的问题。
另一方面，K-means、二分 K-means 和高斯混合这三种算法的一个共同点是，聚类数必须事先确定，并作为参数传递给算法。因此，非正式地说，确定聚类数是一个独立的优化问题，需要解决。
现在，我们将使用基于肘部法则的启发式方法。从 K=2 个聚类开始，然后通过增加 K 值运行 K-means 算法，观察成本函数 WCSS 的变化：
val iterations = 20
for (i <- 2 to iterations) {
val kmeans = new KMeans().setK(i).setSeed(12345L)
val model = kmeans.fit(pcaDF)
val WSSSE = model.computeCost(pcaDF)
println("Within-Cluster Sum of Squares for k = " + i + " is " +
WSSSE)
}
在某些时刻，可以观察到成本函数有一个大幅下降，但随着k值增加，改进变得非常微小。正如聚类分析文献中所建议的，我们可以选择在 WCSS 的最后一次大幅下降之后的k值作为最优值。现在，让我们来看一下在 2 到 20 之间的不同聚类数的 WCSS 值，例如：
Within-Cluster Sum of Squares for k = 2 is 453.161838161838
Within-Cluster Sum of Squares for k = 3 is 438.2392344497606
Within-Cluster Sum of Squares for k = 4 is 390.2278787878787
Within-Cluster Sum of Squares for k = 5 is 397.72112098427874
Within-Cluster Sum of Squares for k = 6 is 367.8890909090908
Within-Cluster Sum of Squares for k = 7 is 362.3360347662672
Within-Cluster Sum of Squares for k = 8 is 347.49306362861336
Within-Cluster Sum of Squares for k = 9 is 327.5002901103624
Within-Cluster Sum of Squares for k = 10 is 327.29376873556436
Within-Cluster Sum of Squares for k = 11 is 315.2954156954155
Within-Cluster Sum of Squares for k = 12 is 320.2478696814693
Within-Cluster Sum of Squares for k = 13 is 308.7674242424241
Within-Cluster Sum of Squares for k = 14 is 314.64784054938576
Within-Cluster Sum of Squares for k = 15 is 297.38523698523704
Within-Cluster Sum of Squares for k = 16 is 294.26114718614707
Within-Cluster Sum of Squares for k = 17 is 284.34890572390555
Within-Cluster Sum of Squares for k = 18 is 280.35662525879917
Within-Cluster Sum of Squares for k = 19 is 272.765762015762
Within-Cluster Sum of Squares for k = 20 is 272.05702362771336
现在让我们讨论如何利用肘部法则来确定聚类数。如下面所示，我们计算了成本函数 WCSS，作为 K-means 算法在选定人群组的 Y 染色体基因变异上的聚类数函数。
可以观察到，当k = 9时出现了一个相对较大的下降（尽管这并不是一个剧烈的下降）。因此，我们选择将聚类数设置为 10，如图 10所示：
图 19：聚类数与 WCSS 的关系
使用 H2O 进行种族预测
到目前为止，我们已经看到如何聚类基因变异。我们还使用了肘部法则，找到了最优k值和初步的聚类数。现在我们应该探索另一个我们一开始计划的任务——即种族预测。
在之前的 K-means 部分中，我们准备了一个名为schemaDF的 Spark 数据框。这个数据框不能直接用于 H2O。然而，我们需要进行额外的转换。我们使用asH2OFrame()方法将 Spark 数据框转换为 H2O 框架：
val dataFrame = h2oContext.asH2OFrame(schemaDF)
现在，有一件重要的事情你需要记住，当使用 H2O 时，如果没有将标签列转换为类别型，它会将分类任务视为回归任务。为了避免这种情况，我们可以使用 H2O 的toCategoricalVec()方法。由于 H2O 框架是弹性的，我们可以进一步更新相同的框架：
dataFrame.replace(dataFrame.find("Region"),
dataFrame.vec("Region").toCategoricalVec()).remove()
dataFrame.update()
现在我们的 H2O 框架已准备好训练基于 H2O 的深度学习模型（即 DNN，或者更具体地说，是深度 MLP）。然而，在开始训练之前，让我们使用 H2O 内置的FrameSplitter()方法随机将数据框拆分为 60%的训练集、20%的测试集和 20%的验证集：
val frameSplitter = new FrameSplitter(
dataFrame, Array(.8, .1), Array("training", "test", "validation")
.map(Key.make[Frame]),null)
water.H2O.submitTask(frameSplitter)
val splits = frameSplitter.getResult
val training = splits(0)
val test = splits(1)
val validation = splits(2)
太棒了！我们的训练集、测试集和验证集已经准备好，现在让我们为我们的深度学习模型设置参数：
// Set the parameters for our deep learning model.
val deepLearningParameters = new DeepLearningParameters()
deepLearningParameters._train = training
deepLearningParameters._valid = validation
deepLearningParameters._response_column = "Region"
deepLearningParameters._epochs = 200
deepLearningParameters._l1 = 0.01
deepLearningParameters._seed = 1234567
deepLearningParameters._activation = Activation.RectifierWithDropout
deepLearningParameters._hidden = ArrayInt
在前面的设置中，我们已经指定了一个具有三层隐藏层的 MLP，隐藏层神经元数分别为 128、256 和 512。因此，总共有五层，包括输入层和输出层。训练将迭代最多 200 次。由于隐藏层神经元数量较多，我们应该使用 dropout 来避免过拟合。为了获得更好的正则化效果，我们使用了 L1 正则化。
前面的设置还表明，我们将使用训练集来训练模型，同时使用验证集来验证训练结果。最后，响应列为Region。另一方面，种子用于确保结果的可复现性。
一切准备就绪！现在，让我们来训练深度学习模型：
val deepLearning = new DeepLearning(deepLearningParameters)
val deepLearningTrained = deepLearning.trainModel
val trainedModel = deepLearningTrained.get
根据你的硬件配置，这可能需要一段时间。因此，现在是时候休息一下，喝杯咖啡了！一旦我们得到训练好的模型，我们可以查看训练误差：
val error = trainedModel.classification_error()
println("Training Error: " + error)
>>>
Training Error: 0.5238095238095238
不幸的是，训练结果并不理想！不过，我们应该尝试不同的超参数组合。尽管误差较高，但让我们不要过于担心，先评估一下模型，计算一些模型指标，并评估模型质量：
val trainMetrics = ModelMetricsSupport.modelMetricsModelMetricsMultinomial
val met = trainMetrics.cm()
println("Accuracy: "+ met.accuracy())
println("MSE: "+ trainMetrics.mse)
println("RMSE: "+ trainMetrics.rmse)
println("R2: " + trainMetrics.r2)
>>>
Accuracy: 0.42105263157894735
MSE: 0.49369297490740655
RMSE: 0.7026328877211816
R2: 0.6091597281983032
准确率不高！不过，你应该尝试其他 VCF 文件，并调整超参数。例如，在减少隐藏层神经元数量、使用 L2 正则化和 100 轮训练后，我的准确率提高了大约 20%：
val deepLearningParameters = new DeepLearningParameters()
deepLearningParameters._train = training
deepLearningParameters._valid = validation
deepLearningParameters._response_column = "Region"
deepLearningParameters._epochs = 100
deepLearningParameters._l2 = 0.01
deepLearningParameters._seed = 1234567
deepLearningParameters._activation = Activation.RectifierWithDropout
deepLearningParameters._hidden = ArrayInt
>>>
Training Error: 0.47619047619047616
Accuracy: 0.5263157894736843
MSE: 0.39112548936806274
RMSE: 0.6254002633258662
R2: 0.690358987583617
另一个改进的线索在这里。除了这些超参数外，使用基于 H2O 的深度学习算法的另一个优势是我们可以获取相对的变量/特征重要性。在前面的章节中，我们已经看到，在 Spark 中使用随机森林算法时，也可以计算变量的重要性。
因此，基本思路是，如果模型表现不佳，值得删除一些不重要的特征并重新进行训练。现在，在有监督训练过程中，我们可以找到特征的重要性。我观察到以下特征重要性：
图 20：使用 H2O 计算的相对特征重要性
现在的问题是，为什么不尝试删除这些特征，然后重新训练并观察准确性是否提高？这个问题留给读者思考。
使用随机森林进行种族预测
在上一节中，我们已经看到了如何使用 H2O 进行种族预测。然而，我们未能取得更好的预测准确性。因此，H2O 目前还不够成熟，无法计算所有必要的性能指标。
那么，为什么不尝试基于 Spark 的树集成技术，比如随机森林或梯度提升树（GBT）呢？因为我们已经看到在大多数情况下，随机森林（RF）表现出更好的预测精度，所以我们就从这个开始试试。
在 K-means 部分，我们已经准备好了名为schemaDF的 Spark DataFrame。因此，我们可以简单地将变量转换成我们之前描述的特征向量。然而，为了实现这一点，我们需要排除标签列。我们可以使用drop()方法来做到这一点，如下所示：
val featureCols = schemaDF.columns.drop(1)
val assembler =
new VectorAssembler()
.setInputCols(featureCols)
.setOutputCol("features")
val assembleDF = assembler.transform(schemaDF).select("features", "Region")
assembleDF.show()
在此阶段，您可以进一步降低维度，并使用 PCA 或任何其他特征选择算法提取最主要的成分。然而，我将把这个留给您决定。由于 Spark 期望标签列是数值型的，我们必须将民族名称转换为数值。我们可以使用StringIndexer()来完成这一操作。这是非常简单的：
val indexer =
new StringIndexer()
.setInputCol("Region")
.setOutputCol("label")
val indexedDF = indexer.fit(assembleDF)
.transform(assembleDF)
.select("features", "label")
然后我们随机划分数据集进行训练和测试。在我们的例子中，假设我们使用 75%作为训练数据，其余作为测试数据：
val seed = 12345L
val splits = indexedDF.randomSplit(Array(0.75, 0.25), seed)
val (trainDF, testDF) = (splits(0), splits(1))
由于这是一个小数据集，考虑到这一点，我们可以缓存训练集和测试集，以便更快地访问：
trainDF.cache
testDF.cache
val rf = new RandomForestClassifier()
.setLabelCol("label")
.setFeaturesCol("features")
.setSeed(1234567L)
现在让我们创建一个paramGrid，用于在决策树的maxDepth参数中进行搜索，以获得最佳模型：
val paramGrid =
new ParamGridBuilder()
.addGrid(rf.maxDepth, 3 :: 5 :: 15 :: 20 :: 25 :: 30 :: Nil)
.addGrid(rf.featureSubsetStrategy, "auto" :: "all" :: Nil)
.addGrid(rf.impurity, "gini" :: "entropy" :: Nil)
.addGrid(rf.maxBins, 3 :: 5 :: 10 :: 15 :: 25 :: 35 :: 45 :: Nil)
.addGrid(rf.numTrees, 5 :: 10 :: 15 :: 20 :: 30 :: Nil)
.build()
val evaluator = new MulticlassClassificationEvaluator()
.setLabelCol("label")
.setPredictionCol("prediction")
接着，我们设置 10 折交叉验证，以获得优化且稳定的模型。这将减少过拟合的可能性：
val numFolds = 10
val crossval =
new CrossValidator()
.setEstimator(rf)
.setEvaluator(evaluator)
.setEstimatorParamMaps(paramGrid)
.setNumFolds(numFolds)
好的，现在我们准备进行训练了。那么让我们使用最佳超参数设置来训练随机森林模型：
val cvModel = crossval.fit(trainDF)
现在我们已经得到了交叉验证和最佳模型，为什么不使用测试集来评估模型呢？为什么不呢？首先，我们为每个实例计算预测 DataFrame。然后我们使用MulticlassClassificationEvaluator()来评估性能，因为这是一个多类分类问题。
此外，我们还计算了性能指标，如准确率、精确率、召回率和F1值。请注意，使用 RF 分类器时，我们可以获得加权精确率和加权召回率：
val predictions = cvModel.transform(testDF)
predictions.show(10)
>>>
图 21：使用随机森林的原始预测概率、真实标签和预测标签
val metric =
new MulticlassClassificationEvaluator()
.setLabelCol("label")
.setPredictionCol("prediction")
val evaluator1 = metric.setMetricName("accuracy")
val evaluator2 = metric.setMetricName("weightedPrecision")
val evaluator3 = metric.setMetricName("weightedRecall")
val evaluator4 = metric.setMetricName("f1")
现在让我们计算分类的准确率、精确率、召回率、F1值以及测试数据上的错误率：
val accuracy = evaluator1.evaluate(predictions)
val precision = evaluator2.evaluate(predictions)
val recall = evaluator3.evaluate(predictions)
val f1 = evaluator4.evaluate(predictions)
最后，我们打印出性能指标：
println("Accuracy = " + accuracy);
println("Precision = " + precision)
println("Recall = " + recall)
println("F1 = " + f1)
println(s"Test Error = ${1 - accuracy}")
>>>
Accuracy = 0.7196470196470195
Precision = 0.7196470196470195
Recall = 0.7196470196470195
F1 = 0.7196470196470195
Test Error = 0.28035298035298046
是的，结果证明它的表现更好。这有点出乎意料，因为我们原本希望从深度学习模型中获得更好的预测准确性，但并没有得到。如我之前所说，我们仍然可以尝试使用 H2O 的其他参数。无论如何，现在我们看到使用随机森林大约提高了 25%。不过，可能仍然可以进一步改进。
总结
在本章中，我们展示了如何与一些大数据工具（如 Spark、H2O 和 ADAM）进行交互，以处理大规模基因组数据集。我们应用了基于 Spark 的 K-means 算法，分析了来自 1000 基因组计划的数据，旨在对人群规模的基因型变异进行聚类。
然后，我们应用基于 H2O 的深度学习（DL）算法和基于 Spark 的随机森林模型来预测地理族群。此外，我们还学习了如何安装和配置 H2O 进行深度学习。这些知识将在后续章节中使用。最后，也是最重要的，我们学会了如何使用 H2O 计算变量重要性，以便选择训练集中最重要的特征。
在下一章中，我们将看到如何有效地使用潜在狄利克雷分配（LDA）算法来发现数据中的有用模式。我们将比较其他主题建模算法以及 LDA 的可扩展性。同时，我们还将利用自然语言处理（NLP）库，如斯坦福 NLP。
第五章：主题建模 - 对大规模文本的深入理解
主题建模（TM）是一种广泛用于从大量文档中挖掘文本的技术。通过这些主题，可以总结和组织包含主题词及其相对权重的文档。这个项目将使用的数据集是纯文本格式，未经过结构化处理。
我们将看到如何有效地使用潜在狄利克雷分配（LDA）算法来发现数据中的有用模式。我们将比较其他 TM 算法及 LDA 的可扩展性。此外，我们还将利用自然语言处理（NLP）库，如斯坦福 NLP。
简而言之，在这个端到端的项目中，我们将学习以下主题：
-
主题建模与文本聚类
-
LDA 算法是如何工作的？
-
使用 LDA、Spark MLlib 和标准 NLP 进行主题建模
-
其他主题模型与 LDA 的可扩展性测试
-
模型部署
主题建模与文本聚类
在 TM 中，主题是通过一组词汇来定义的，每个词汇在该主题下都有一个出现的概率，不同的主题有各自的词汇集合及其相应的概率。不同的主题可能共享一些词汇，而一个文档可能与多个主题相关联。简而言之，我们有一组文本数据集，即一组文本文件。现在，挑战在于使用 LDA 从数据中发现有用的模式。
有一种流行的基于 LDA 的 TM 方法，其中每个文档被视为多个主题的混合，每个文档中的词汇被认为是从文档的主题中随机抽取的。这些主题被视为隐藏的，必须通过分析联合分布来揭示，从而计算给定观察变量和文档中的词语的条件分布（主题）。TM 技术广泛应用于从大量文档中挖掘文本的任务。这些主题随后可以用来总结和组织包含主题词及其相对权重的文档（见图 1）：
图 1：TM 概述（来源：Blei, D.M.等，概率主题模型，ACM 通信，55(4)，77-84，2012）
如前图所示，主题的数量远小于与文档集合相关联的词汇量，因此主题空间的表示可以被看作是一种降维过程：
图 2：TM 与文本聚类的对比
与 TM 相比，在文档聚类中，基本思想是根据一种广为人知的相似性度量，将文档分组。为了进行分组，每个文档都由一个向量表示，该向量表示文档中词汇的权重。
通常使用词频-逆文档频率（也称为 TF-IDF 方案）来进行加权。聚类的最终结果是一个簇的列表，每个文档出现在某一个簇中。TM 和文本聚类的基本区别可以通过下图来说明：
LDA 算法是如何工作的？
LDA 是一个主题模型，用于从一组文本文档中推断主题。LDA 可以被看作是一个聚类算法，其中主题对应于聚类中心，文档对应于数据集中的实例（行）。主题和文档都存在于特征空间中，其中特征向量是词频向量（词袋）。LDA 不是通过传统的距离估计聚类，而是使用基于文本文档生成的统计模型的函数（见 图 3）：
图 3：LDA 算法在一组文档上的工作原理
具体来说，我们希望讨论人们在大量文本中最常谈论的话题。自 Spark 1.3 发布以来，MLlib 支持 LDA，这是文本挖掘和自然语言处理领域中最成功使用的主题模型（TM）技术之一。
此外，LDA 还是第一个采用 Spark GraphX 的 MLlib 算法。以下术语是我们正式开始 TM 应用之前值得了解的：
-
"word" = "term"：词汇表中的一个元素 -
"token"：出现在文档中的术语实例 -
"topic"：表示某一概念的词汇的多项式分布
在 Spark 中开发的基于 RDD 的 LDA 算法是一个为文本文档设计的主题模型。它基于原始的 LDA 论文（期刊版）：Blei, Ng, 和 Jordan，Latent Dirichlet Allocation，JMLR，2003。
此实现通过 setOptimizer 函数支持不同的推理算法。EMLDAOptimizer 使用 期望最大化（EM）在似然函数上进行聚类学习，并提供全面的结果，而 OnlineLDAOptimizer 使用迭代的小批量采样进行在线变分推理，并且通常对内存友好。
EM 是一种迭代方式，用于逼近最大似然函数。在实际应用中，当输入数据不完整、缺失数据点或存在隐藏潜在变量时，最大似然估计可以找到 最佳拟合 模型。
LDA 输入一组文档，作为词频向量，并使用以下参数（通过构建器模式设置）：
-
K：主题的数量（即聚类中心的数量）（默认值是 10）。 -
ldaOptimizer：用于学习 LDA 模型的优化器，可以是EMLDAOptimizer或OnlineLDAOptimizer（默认是EMLDAOptimizer）。 -
Seed：用于可重复性的随机种子（虽然是可选的）。 -
docConcentration：文档主题分布的 Dirichlet 参数。较大的值会鼓励推断出的分布更平滑（默认值是Vectors.dense(-1)）。 -
topicConcentration: Drichilet 参数，用于先验主题在术语（单词）分布上的分布。较大的值确保推断分布更加平滑（默认为 -1）。 -
maxIterations: 迭代次数的限制（默认为 20）。 -
checkpointInterval: 如果使用检查点（在 Spark 配置中设置），此参数指定创建检查点的频率。如果maxIterations很大，使用检查点可以帮助减少磁盘上的洗牌文件大小，并有助于故障恢复（默认为 10）。
图 4：主题分布及其外观
让我们看一个例子。假设篮子里有 n 个球，有 w 种不同的颜色。现在假设词汇表中的每个术语都有 w 种颜色之一。现在还假设词汇表中的术语分布在 m 个主题中。现在篮子中每种颜色出现的频率与对应术语在主题 φ 中的权重成比例。
然后，LDA 算法通过使每个球的大小与其对应术语的权重成比例来包含一个术语加权方案。在 图 4 中，n 个术语在一个主题中具有总权重，例如主题 0 到 3。图 4 展示了从随机生成的 Twitter 文本中的主题分布。
现在我们已经看到通过使用 TM，我们可以在非结构化的文档集合中找到结构。一旦 发现 结构，如 图 4 所示，我们可以回答几个问题，如下所示：
-
文档 X 是关于什么的？
-
文档 X 和 Y 有多相似？
-
如果我对主题 Z 感兴趣，我应该先阅读哪些文档？
在接下来的部分中，我们将看到一个使用基于 Spark MLlib 的 LDA 算法的 TM 示例，以回答前面的问题。
使用 Spark MLlib 和 Stanford NLP 进行主题建模
在本小节中，我们使用 Spark 表示了一种半自动化的 TM 技术。在模型重用和部署阶段，我们将使用从 GitHub 下载的数据集训练 LDA，位于 github.com/minghui/Twitter-LDA/tree/master/data/Data4Model/test。然而，稍后在本章节中，我们将使用更知名的文本数据集。
实施
以下步骤展示了从数据读取到打印主题的 TM 过程，以及它们的术语权重。以下是 TM 管道的简短工作流程：
object topicmodelingwithLDA {
def main(args: Array[String]): Unit = {
val lda =
new LDAforTM()
// actual computations are done here
val defaultParams = Params().copy(input = "data/docs/") //Loading parameters for training
lda.run(defaultParams)
// Training the LDA model with the default parameters.
}
}
我们还需要导入一些相关的包和库：
import edu.stanford.nlp.process.Morphology
import edu.stanford.nlp.simple.Document
import org.apache.log4j.{Level, Logger}
import scala.collection.JavaConversions._
import org.apache.spark.{SparkConf, SparkContext}
import org.apache.spark.ml.Pipeline
import org.apache.spark.ml.feature._
import org.apache.spark.ml.linalg.{Vector => MLVector}
import org.apache.spark.mllib.clustering.{DistributedLDAModel, EMLDAOptimizer, LDA, OnlineLDAOptimizer, LDAModel}
import org.apache.spark.mllib.linalg.{ Vector, Vectors }
import org.apache.spark.rdd.RDD
import org.apache.spark.sql.{Row, SparkSession}
实际上，TM 的计算是在 LDAforTM 类中完成的。Params 是一个案例类，用于加载用于训练 LDA 模型的参数。最后，我们通过 Params 类设置的参数来训练 LDA 模型。现在我们将详细解释每个步骤，逐步源代码：
第一步 - 创建一个 Spark 会话
让我们通过以下方式创建一个 Spark 会话：定义计算核心数量、SQL 仓库和应用程序名称。
val spark = SparkSession
.builder
.master("local[*]")
.config("spark.sql.warehouse.dir", "C:/data/")
.appName(s"LDA")
.getOrCreate()
步骤 2 - 创建词汇表和标记计数以训练 LDA 模型，在文本预处理后进行
run() 方法接受 params，如输入文本、预定义的词汇表大小和停用词文件：
def run(params: Params)
然后，它将开始为 LDA 模型进行文本预处理，如下所示（即在 run 方法内部）：
// Load documents, and prepare them for LDA.
val preprocessStart = System.nanoTime()
val (corpus, vocabArray, actualNumTokens) = preprocess(params.input, params.vocabSize, params.stopwordFile)
Params case 类用于定义训练 LDA 模型的参数。代码如下所示：
//Setting the parameters before training the LDA model
case class Params(var input: String = "", var ldaModel: LDAModel = null,
k: Int = 5,
maxIterations: Int = 100,
docConcentration: Double = 5,
topicConcentration: Double = 5,
vocabSize: Int = 2900000,
stopwordFile: String = "data/docs/stopWords.txt",
algorithm: String = "em",
checkpointDir: Option[String] = None,
checkpointInterval: Int = 100)
为了获得更好的结果，你需要通过反复试验来设置这些参数。或者，你可以选择交叉验证以获得更好的性能。如果你想检查当前的参数，可以使用以下代码：
if (params.checkpointDir.nonEmpty) {
spark.sparkContext.setCheckpointDir(params.checkpointDir.get)
}
preprocess 方法用于处理原始文本。首先，让我们使用 wholeTextFiles() 方法读取整个文本，如下所示：
val initialrdd = spark.sparkContext.wholeTextFiles(paths).map(_._2)
initialrdd.cache()
在前面的代码中，paths 是文本文件的路径。然后，我们需要根据 lemma 文本准备从原始文本中提取的形态学 RDD，如下所示：
val rdd = initialrdd.mapPartitions { partition =>
val morphology = new Morphology()
partition.map { value => helperForLDA.getLemmaText(value, morphology) }
}.map(helperForLDA.filterSpecialCharacters)
在这里，helperForLDA 类中的 getLemmaText() 方法提供了经过过滤特殊字符后的 lemma 文本，如 ("""[! @ # $ % ^ & * ( ) _ + - − , " ' ; : . ` ? --])，作为正则表达式，使用 filterSpecialCharacters() 方法。方法如下所示：
def getLemmaText(document: String, morphology: Morphology) = {
val string =
new StringBuilder()
val value =
new Document(document).sentences().toList.flatMap {
a =>
val words = a.words().toList
val tags = a.posTags().toList
(words zip tags).toMap.map {
a =>
val newWord = morphology.lemma(a._1, a._2)
val addedWoed =
if (newWord.length > 3) {
newWord
}
else { "" }
string.append(addedWoed + " ")
}
}
string.toString()
}
需要注意的是，Morphology() 类通过移除仅有屈折变化的部分（而非派生形态学）来计算英语单词的基本形式。也就是说，它只处理名词复数、代词的格、动词词尾，而不涉及比较级形容词或派生名词等内容。getLemmaText() 方法接收文档及其对应的形态学信息，并最终返回已词干化的文本。
这个来自斯坦福 NLP 小组。要使用它，你需要在主类文件中添加以下导入：edu.stanford.nlp.process.Morphology。在 pom.xml 文件中，你需要包含以下条目作为依赖项：
<dependency>
<groupId>edu.stanford.nlp</groupId>
<artifactId>stanford-corenlp</artifactId>
<version>3.6.0</version>
</dependency>
<dependency>
<groupId>edu.stanford.nlp</groupId>
<artifactId>stanford-corenlp</artifactId>
<version>3.6.0</version>
<classifier>models</classifier>
</dependency>
filterSpecialCharacters() 方法如下所示：
def filterSpecialCharacters(document: String) = document.replaceAll("""[! @ # $ % ^ & * ( ) _ + - − , " ' ; : . ` ? --]""", " ")
一旦我们得到了移除特殊字符的 RDD，就可以创建 DataFrame 来构建文本分析管道：
rdd.cache()
initialrdd.unpersist()
val df = rdd.toDF("docs")
df.show()
DataFrame 只包含文档标签。DataFrame 的快照如下：
图 5：来自输入数据集的原始文本
现在，如果你仔细观察前面的 DataFrame，你会发现我们仍然需要对其进行分词。此外，DataFrame 中包含停用词，如 this、with 等，因此我们还需要将它们移除。首先，让我们使用 RegexTokenizer API 对其进行分词，如下所示：
val tokenizer = new RegexTokenizer()
.setInputCol("docs")
.setOutputCol("rawTokens")
现在让我们按照以下方式移除所有停用词：
val stopWordsRemover = new StopWordsRemover()
.setInputCol("rawTokens")
.setOutputCol("tokens")
stopWordsRemover.setStopWords(stopWordsRemover.getStopWords ++ customizedStopWords)
此外，我们还需要应用计数向量，以便从标记中找到仅重要的特征。这将有助于将管道阶段链接在一起。我们按照以下方式进行：
val countVectorizer = new CountVectorizer()
.setVocabSize(vocabSize)
.setInputCol("tokens")
.setOutputCol("features")
当先验词典不可用时，可以使用CountVectorizer作为估算器来提取词汇并生成CountVectorizerModel。换句话说，CountVectorizer用于将一组文本文档转换为标记（即术语）计数的向量。CountVectorizerModel为文档提供稀疏表示，这些文档会基于词汇表然后输入到 LDA 模型中。从技术上讲，当调用fit()方法进行拟合时，CountVectorizer会选择按术语频率排序的前vocabSize个词。
现在，通过链接转换器（tokenizer、stopWordsRemover和countVectorizer）来创建管道，代码如下：
val pipeline = new Pipeline().setStages(Array(tokenizer, stopWordsRemover, countVectorizer))
现在，让我们将管道拟合并转换为词汇表和标记数：
val model = pipeline.fit(df)
val documents = model.transform(df).select("features").rdd.map {
case Row(features: MLVector) => Vectors.fromML(features)
}.zipWithIndex().map(_.swap)
最后，返回词汇和标记计数对，代码如下：
(documents, model.stages(2).asInstanceOf[CountVectorizerModel].vocabulary, documents.map(_._2.numActives).sum().toLong) Now let's see the statistics of the training data:
println() println("Training corpus summary:")
println("-------------------------------")
println("Training set size: " + actualCorpusSize + " documents")
println("Vocabulary size: " + actualVocabSize + " terms")
println("Number of tockens: " + actualNumTokens + " tokens")
println("Preprocessing time: " + preprocessElapsed + " sec")
println("-------------------------------")
println()
>>>
Training corpus summary:
-------------------------------
Training set size: 19 documents
Vocabulary size: 21611 terms
Number of tockens: 75784 tokens
Preprocessing time: 46.684682086 sec
第 3 步 - 在训练之前实例化 LDA 模型
在开始训练 LDA 模型之前，让我们实例化 LDA 模型，代码如下：
val lda = new LDA()
第 4 步 - 设置 NLP 优化器
为了从 LDA 模型获得更好的优化结果，我们需要设置一个包含 LDA 算法的优化器，它执行实际的计算，并存储算法的内部数据结构（例如图形或矩阵）及其他参数。
在这里，我们使用了EMLDAOPtimizer优化器。你也可以使用OnlineLDAOptimizer()优化器。EMLDAOPtimizer存储一个数据+参数图形，以及算法参数。其底层实现使用 EM（期望最大化）。
首先，让我们通过添加(1.0 / actualCorpusSize)以及一个非常低的学习率（即 0.05）到MiniBatchFraction来实例化EMLDAOptimizer，以便在像我们这样的小数据集上收敛训练，代码如下：
val optimizer = params.algorithm.toLowerCase
match {
case "em" =>
new EMLDAOptimizer
// add (1.0 / actualCorpusSize) to MiniBatchFraction be more robust on tiny datasets.
case "online" =>
new OnlineLDAOptimizer().setMiniBatchFraction(0.05 + 1.0 / actualCorpusSize)
case _ =>
thrownew IllegalArgumentException("Only em, online are supported but got
${params.algorithm}.")
}
现在，使用 LDA API 中的setOptimizer()方法设置优化器，代码如下：
lda.setOptimizer(optimizer)
.setK(params.k)
.setMaxIterations(params.maxIterations)
.setDocConcentration(params.docConcentration)
.setTopicConcentration(params.topicConcentration)
.setCheckpointInterval(params.checkpointInterval)
第 5 步 - 训练 LDA 模型
让我们开始使用训练语料库训练 LDA 模型，并跟踪训练时间，代码如下：
val startTime = System.nanoTime()
ldaModel = lda.run(corpus)
val elapsed = (System.nanoTime() - startTime) / 1e9
println("Finished training LDA model. Summary:")
println("Training time: " + elapsed + " sec")
现在，此外，我们可以保存训练好的模型以供将来重用，保存的代码如下：
//Saving the model for future use
params.ldaModel.save(spark.sparkContext, "model/LDATrainedModel")
请注意，一旦完成训练并获得最优训练效果，部署模型之前需要取消注释前面的行。否则，它将在模型重用阶段因抛出异常而被停止。
对于我们拥有的文本，LDA 模型训练耗时 6.309715286 秒。请注意，这些时间代码是可选的。我们仅提供它们供参考，以便了解训练时间：
第 6 步 - 准备感兴趣的主题
准备前 5 个主题，每个主题包含 10 个术语。包括术语及其对应的权重：
val topicIndices = ldaModel.describeTopics(maxTermsPerTopic = 10)
println(topicIndices.length)
val topics = topicIndices.map {
case (terms, termWeights) => terms.zip(termWeights).map {
case (term, weight) => (vocabArray(term.toInt), weight)
}
}
第 7 步 - 主题建模
打印出前 10 个主题，展示每个主题的权重最高的词汇。同时，还包括每个主题的总权重，代码如下：
var sum = 0.0
println(s"${params.k} topics:")
topics.zipWithIndex.foreach {
case (topic, i) =>
println(s"TOPIC $i")
println("------------------------------")
topic.foreach {
case (term, weight) =>
term.replaceAll("\s", "")
println(s"$termt$weight")
sum = sum + weight
}
println("----------------------------")
println("weight: " + sum)
println()
现在让我们看看 LDA 模型在主题建模方面的输出：
5 topics:
TOPIC 0
------------------------------
come 0.0070183359426213635
make 0.006893251344696077
look 0.006629265338364568
know 0.006592594912464674
take 0.006074234442310174
little 0.005876330712306203
think 0.005153843469004155
time 0.0050685675513282525
hand 0.004524837827665401
well 0.004224698942533204
----------------------------
weight: 0.05805596048329406
TOPIC 1
------------------------------
thus 0.008447268016707914
ring 0.00750959344769264
fate 0.006802070476284118
trojan 0.006310545607626158
bear 0.006244268350438889
heav 0.005479939900136969
thro 0.005185211621694439
shore 0.004618008184651363
fight 0.004161178536600401
turnus 0.003899151842042464
----------------------------
weight: 0.11671319646716942
TOPIC 2
------------------------------
aladdin 7.077183389325728E-4
sultan 6.774311890861097E-4
magician 6.127791175835228E-4
genie 6.06094509479989E-4
vizier 6.051618911188781E-4
princess 5.654756758514474E-4
fatima 4.050749957608771E-4
flatland 3.47788388834721E-4
want 3.4263963705536023E-4
spaceland 3.371784715458026E-4
----------------------------
weight: 0.1219205386824187
TOPIC 3
------------------------------
aladdin 7.325869707607238E-4
sultan 7.012354862373387E-4
magician 6.343184784726607E-4
genie 6.273921840260785E-4
vizier 6.264266945018852E-4
princess 5.849046214967484E-4
fatima 4.193089052802858E-4
flatland 3.601371993827707E-4
want 3.5398019331108816E-4
spaceland 3.491505202713831E-4
----------------------------
weight: 0.12730997993615964
TOPIC 4
------------------------------
captain 0.02931475169407467
fogg 0.02743105575940755
nautilus 0.022748371008515483
passepartout 0.01802140608022664
nemo 0.016678258146358142
conseil 0.012129894049747918
phileas 0.010441664411654412
canadian 0.006217638883315841
vessel 0.00618937301246955
land 0.00615311666365297
----------------------------
weight: 0.28263550964558276
从上述输出中，我们可以看到输入文档的主题五占据了最大权重，权重为 0.28263550964558276。该主题涉及的词汇包括 captain、fogg、nemo、vessel 和 land 等。
第 8 步 - 测量两个文档的似然度
现在为了获取一些其他统计数据，例如文档的最大似然或对数似然，我们可以使用以下代码：
if (ldaModel.isInstanceOf[DistributedLDAModel]) {
val distLDAModel = ldaModel.asInstanceOf[DistributedLDAModel]
val avgLogLikelihood = distLDAModel.logLikelihood / actualCorpusSize.toDouble
println("The average log likelihood of the training data: " +
avgLogLikelihood)
println()
}
上述代码计算了 LDA 模型的平均对数似然度，作为 LDA 分布式版本的一个实例：
The average log likelihood of the training data: -209692.79314860413
有关似然度测量的更多信息，感兴趣的读者可以参考en.wikipedia.org/wiki/Likelihood_function。
现在假设我们已经计算了文档 X 和 Y 的前述度量。然后我们可以回答以下问题：
- 文档 X 和 Y 有多相似？
关键在于，我们应该尝试从所有训练文档中获取最低的似然值，并将其作为前述比较的阈值。最后，回答第三个也是最后一个问题：
- 如果我对主题 Z 感兴趣，应该先阅读哪些文档？
最小化的答案：通过仔细查看主题分布和相关词汇的权重，我们可以决定首先阅读哪个文档。
LDA 的可扩展性与其他主题模型的比较
在整个端到端的项目中，我们使用了 LDA，它是最流行的文本挖掘主题建模算法之一。我们还可以使用更强大的主题建模算法，如 概率潜在语义分析（pLSA）、帕金科分配模型（PAM）和 层次狄利克雷过程（HDP）算法。
然而，pLSA 存在过拟合问题。另一方面，HDP 和 PAM 是更复杂的主题建模算法，用于处理复杂的文本挖掘任务，例如从高维文本数据或非结构化文本文档中挖掘主题。最后，非负矩阵分解是另一种在文档集合中寻找主题的方法。无论采用哪种方法，所有主题建模算法的输出都是一个包含相关词汇簇的主题列表。
上述示例展示了如何使用 LDA 算法作为独立应用程序来执行主题建模。LDA 的并行化并不简单，许多研究论文提出了不同的策略。关键障碍在于所有方法都涉及大量的通信。
根据 Databricks 网站上的博客（databricks.com/blog/2015/03/25/topic-modeling-with-lda-mllib-meets-graphx.html），以下是实验过程中使用的数据集以及相关的训练和测试集的统计数据：
-
训练集大小：460 万份文档
-
词汇表大小：110 万个词条
-
训练集大小：11 亿个词条（约 239 个词/文档）
-
100 个主题
-
例如，16 个工作节点的 EC2 集群，可以选择 M4.large 或 M3.medium，具体取决于预算和需求
对于上述设置，经过 10 次迭代，平均时间结果为 176 秒/迭代。从这些统计数据可以看出，LDA 对于非常大的语料库也是非常具有可扩展性的。
部署训练好的 LDA 模型
对于这个小型部署，我们使用一个现实生活中的数据集：PubMed。包含 PubMed 术语的示例数据集可以从以下链接下载：nlp.stanford.edu/software/tmt/tmt-0.4/examples/pubmed-oa-subset.csv。这个链接实际上包含一个 CSV 格式的数据集，但其文件名为奇怪的4UK1UkTX.csv。
更具体地说，该数据集包含了一些生物学文章的摘要、它们的出版年份以及序列号。以下图像给出了一些示例：
图 6：示例数据集的快照
在以下代码中，我们已经保存了训练好的 LDA 模型以备未来使用，如下所示：
params.ldaModel.save(spark.sparkContext, "model/LDATrainedModel")
训练好的模型将被保存到之前提到的位置。该目录将包含模型和训练本身的数据以及元数据，如下图所示：
图 7：训练和保存的 LDA 模型的目录结构
如预期的那样，数据文件夹中包含了一些 parquet 文件，这些文件包含全球主题、它们的计数、标记及其计数，以及主题与相应的计数。现在，接下来的任务是恢复相同的模型，如下所示：
//Restoring the model for reuse
val savedLDAModel = DistributedLDAModel.load(spark.sparkContext, "model/LDATrainedModel/")
//Then we execute the following workflow:
val lda = new LDAforTM()
// actual computations are done here
// Loading the parameters to train the LDA model
val defaultParams = Params().copy(input = "data/4UK1UkTX.csv", savedLDAModel)
lda.run(defaultParams)
// Training the LDA model with the default parameters.
spark.stop()
>>>
Training corpus summary:
-------------------------------
Training set size: 1 documents
Vocabulary size: 14670 terms
Number of tockens: 14670 tokens
Preprocessing time: 12.921435786 sec
-------------------------------
Finished training LDA model.
Summary:
Training time: 23.243336895 sec
The average log likelihood of the training data: -1008739.37857908
5 topics:
TOPIC 0
------------------------------
rrb 0.015234818404037585
lrb 0.015154125349208018
sequence 0.008924621534990771
gene 0.007391453509409655
cell 0.007020265462594214
protein 0.006479622004524878
study 0.004954523307983932
show 0.0040023453035193685
site 0.0038006126784248945
result 0.0036634344941610534
----------------------------
weight: 0.07662582204885438
TOPIC 1
------------------------------
rrb 1.745030693927338E-4
lrb 1.7450110447001028E-4
sequence 1.7424254444446083E-4
gene 1.7411236867642102E-4
cell 1.7407234230511066E-4
protein 1.7400587965300172E-4
study 1.737407317498879E-4
show 1.7347354627656383E-4
site 1.7339989737227756E-4
result 1.7334522348574853E-4
---------------------------
weight: 0.07836521875668061
TOPIC 2
------------------------------
rrb 1.745030693927338E-4
lrb 1.7450110447001028E-4
sequence 1.7424254444446083E-4
gene 1.7411236867642102E-4
cell 1.7407234230511066E-4
protein 1.7400587965300172E-4
study 1.737407317498879E-4
show 1.7347354627656383E-4
site 1.7339989737227756E-4
result 1.7334522348574853E-4
----------------------------
weight: 0.08010461546450684
TOPIC 3
------------------------------
rrb 1.745030693927338E-4
lrb 1.7450110447001028E-4
sequence 1.7424254444446083E-4
gene 1.7411236867642102E-4
cell 1.7407234230511066E-4
protein 1.7400587965300172E-4
study 1.737407317498879E-4
show 1.7347354627656383E-4
site 1.7339989737227756E-4
result 1.7334522348574853E-4
----------------------------
weight: 0.08184401217233307
TOPIC 4
------------------------------
rrb 1.745030693927338E-4
lrb 1.7450110447001028E-4
sequence 1.7424254444446083E-4
gene 1.7411236867642102E-4
cell 1.7407234230511066E-4
protein 1.7400587965300172E-4
study 1.737407317498879E-4
show 1.7347354627656383E-4
site 1.7339989737227756E-4
result 1.7334522348574853E-4
----------------------------
weight: 0.0835834088801593
做得好！我们成功地重新使用了模型并进行了相同的预测。但由于数据的随机性，我们观察到略有不同的预测结果。让我们看看完整的代码以便更清晰地了解：
package com.packt.ScalaML.Topicmodeling
import org.apache.spark.sql.SparkSession
import org.apache.spark.mllib.clustering.{DistributedLDAModel, LDA}
object LDAModelReuse {
def main(args: Array[String]): Unit = {
val spark = SparkSession
.builder
.master("local[*]")
.config("spark.sql.warehouse.dir", "data/")
.appName(s"LDA_TopicModelling")
.getOrCreate()
//Restoring the model for reuse
val savedLDAModel = DistributedLDAModel.load(spark.sparkContext, "model/LDATrainedModel/")
val lda = new LDAforTM()
// actual computations are done here
val defaultParams = Params().copy(input = "data/4UK1UkTX.csv", savedLDAModel)
//Loading params
lda.run(defaultParams)
// Training the LDA model with the default parameters.
spark.stop()
}
}
总结
在本章中，我们已经看到如何高效地使用和结合 LDA 算法与自然语言处理库，如斯坦福 NLP，从大规模文本中发现有用的模式。我们还进行了 TM 算法与 LDA 可扩展性之间的对比分析。
最后，对于一个现实的示例和用例，有兴趣的读者可以参考以下博客文章：blog.codecentric.de/en/2017/01/topic-modeling-codecentric-blog-articles/。
Netflix 是一家美国娱乐公司，由 Reed Hastings 和 Marc Randolph 于 1997 年 8 月 29 日在加利福尼亚州斯科茨谷成立。它专门提供流媒体媒体和按需视频服务，通过在线和 DVD 邮寄的方式提供。在 2013 年，Netflix 扩展到电影和电视制作，以及在线分发。Netflix 为其订阅者使用基于模型的协同过滤方法进行实时电影推荐。
在下一章，我们将看到两个端到端的项目：一个基于项目的协同过滤电影相似度测量，以及一个基于模型的电影推荐引擎，利用 Spark 为新用户推荐电影。我们将看到如何在这两个可扩展的电影推荐引擎中实现ALS与矩阵分解的互操作。
扫一扫
6483
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
