机器学习实践秘籍-二-

机器学习实践秘籍（二）

原文：annas-archive.org/md5/15b58dbdde87ea8cb4894e5d29bf0db5

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：模型选择和正则化

在本章中，我们将介绍以下内容：

• 收敛方法 - 每日燃烧的卡路里

• 维度缩减方法 - Delta 的飞机机队

• 主成分分析 - 理解世界美食

引言

子集选择：在机器学习中，监督分类的主要挑战之一是使用标记示例来诱导一个将对象分类到有限已知类别的模型。数值或名义特征向量用于描述各种示例。在特征子集选择问题中，学习算法面临的问题是在选择一些特征子集上集中注意力，同时忽略其余部分。

当拟合线性回归模型时，我们感兴趣的变量子集是最好描述数据的。在寻找变量集时，可以采用多种不同的策略来选择最佳子集。如果有 m 个变量，并且最佳回归模型由 p 个变量组成，p≤m，那么选择最佳子集的更通用方法可能是尝试所有可能的 p 个变量的组合，并选择最适合数据的模型。

然而，存在 m! p!(m−p)! 种可能的组合，随着 m 值的增加而增加，例如，m = 20 和 p = 4 会产生 4,845 种可能的组合。此外，通过使用更少的特点，我们可以降低获取数据成本并提高分类模型的易理解性。

收敛方法：收敛回归指的是回归情况下的估计或预测的收敛方法；当回归变量之间存在多重共线性时很有用。在数据集相对于研究的协变量数量较小的情况下，收敛技术可以提高预测。常见的收敛方法如下：

• 线性收敛因子--以相同的因子收缩所有系数

• 岭回归--惩罚最大似然，惩罚因子添加到似然函数中，使得系数根据每个协变量的方差单独收缩

• Lasso--通过在标准化协变量的系数绝对值之和上设置约束，将一些系数收缩到零

收敛方法保留预测变量的一部分，同时丢弃其余部分。子集选择产生一个可解释的模型，并且可能比全模型产生更低的预测误差，而不会降低全模型的预测误差。收敛方法更连续，并且不会像高变异性那样受到很大影响。当线性回归模型中有许多相关变量时，它们的系数难以确定，并且表现出高方差。

降维方法：在包括模式识别、数据压缩、机器学习和数据库导航在内的广泛信息处理领域中，降维是一个重要的挑战。测量的数据向量是高维的，在许多情况下，数据位于一个低维流形附近。高维数据的主要挑战是它们是多元的；它们间接地测量了底层来源，这通常不能直接测量。降维也可以被视为推导出一组自由度的过程，这些自由度可以用来再现数据集的大部分变异性。

收缩方法 - 每日燃烧卡路里

为了比较人类的代谢率，基础代谢率（BMR）的概念在临床环境中至关重要，作为确定人类甲状腺状态的手段。哺乳动物的基础代谢率与体重成正比，与场代谢率具有相同的异速增长指数，以及许多生理和生化速率。Fitbit 作为一种设备，使用基础代谢率和一天中进行的活动来估计一天中燃烧的卡路里。

准备中

为了执行收缩方法，我们将使用从 Fitbit 收集的数据集和燃烧卡路里数据集。

第 1 步 - 收集和描述数据

应使用标题为fitbit_export_20160806.csv的 CSV 格式数据集。数据集是标准格式。有 30 行数据，10 个变量。数值变量如下：

• Calories Burned

• Steps

• Distance

• Floors

• Minutes Sedentary

• Minutes Lightly Active

• Minutes Fairly Active

• ExAng

• Minutes Very Active

• Activity Calories

非数值变量如下：

• Date

如何操作...

让我们深入了解。

第 2 步 - 探索数据

作为第一步，需要加载以下包：

    > install.packages("glmnet")
    > install.packages("dplyr")
    > install.packages("tidyr")
    > install.packages("ggplot2")
    > install.packages("caret")
    > install.packages("boot")
    > install.packages("RColorBrewer")
    > install.packages("Metrics")
    > library(dplyr)
    > library(tidyr)
    > library(ggplot2)
    > library(caret)
    > library(glmnet)
    > library(boot)
    > library(RColorBrewer)
    > library(Metrics)

备注

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）中进行了测试

让我们探索数据并了解变量之间的关系。我们将从导入名为fitbit_export_20160806.csv的 csv 数据文件开始。我们将把数据保存到fitbit_details框架中：

> fitbit_details <- read.csv("https://raw.githubusercontent.com/ellisp/ellisp.github.io/source/data/fitbit_export_20160806.csv", 
    + skip = 1, stringsAsFactors = FALSE) %>%
    + mutate(
    + Calories.Burned = as.numeric(gsub(",", "", Calories.Burned)),
    + Steps = as.numeric(gsub(",", "", Steps)),
    + Activity.Calories = as.numeric(gsub(",", "", Activity.Calories)),
    + Date = as.Date(Date, format = "%d/%m/%Y")
    + )

将fitbit_details数据框存储到fitbit数据框中：

> fitbit <- fitbit_details

打印fitbit数据框。head()函数返回fitbit数据框的第一部分。fitbit数据框作为输入参数传递：

 > head(fitbit)

结果如下：

将Activity.Calories和Date值设置为 NULL：

> fitbit$Activity.Calories <- NULL

> fitbit$Date <- NULL

将缩放系数设置为每千步卡路里。然后将结果设置为fitbit$Steps数据框：

> fitbit$Steps <- fitbit$Steps / 1000

打印fitbit$Steps数据框：

> fitbit$Steps

结果如下：

探索所有候选变量。计算相关系数的函数：

    > panel_correlations <- function(x, y, digits = 2, prefix = "", cex.cor, ...){
    # combining multiple plots into one overall graph
    + usr <- par("usr")
    + on.exit(par(usr))
    + par(usr = c(0, 1, 0, 1))
    # computing the absolute value
    + r <- abs(cor(x, y))
# Formatting object 
    + txt <- format(c(r, 0.123456789), digits = digits)[1]
    + txt <- paste0(prefix, txt)
    + if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    + text(0.5, 0.5, txt, cex = cex.cor * r)
    + }

生成散点矩阵。pairs()函数以矩阵形式生成散点图。"fitbit"是散点图的数据集。距离可以直接从"Steps"中几乎精确地计算出来：

> pairs(fitbit[ , -1], lower.panel = panel_correlations, main = "Pairwise Relationship - Fitbit's Measured Activities")

结果如下：

打印fitbit 数据框：

> ggplot(fitbit, aes(x = Distance / Steps)) + geom_rug() + geom_density() +ggtitle("Stride Length Reverse- Engineered from Fitbit Data", subtitle = "Not all strides identical, due to rounding or other jitter")

结果如下：

步骤 3 - 构建模型

使用"Steps"作为唯一解释变量和"Calories.Burned"作为响应变量构建普通最小二乘估计。"lm()"作为一个函数用于拟合线性模型。"Calories.Burned ~ Steps"是公式，而"fitbit"是数据框。结果存储在moderate数据框中：

> moderate <- lm(Calories.Burned ~ Steps, data = fitbit)

打印moderate数据框：

> moderate

结果如下：

将moderate数据框的值四舍五入：

> round(coef(moderate))

结果如下：

使用模型中的残差绘制预测卡路里。plot()函数是一个通用的绘图 R 对象的函数。moderate数据框作为函数值传递。bty参数确定围绕绘图绘制的框的类型：

> plot(moderate, which = 1, bty = "l", main = "Predicted Calories compared with Residuals")

结果如下：

检查残差的偏自相关函数。"pacf()"用于计算偏自相关。"resid()"作为一个函数计算因变量观测数据之间的差异。"moderate"作为一个数据框传递给"resid()"函数，以计算因变量观测数据之间的差异：

> pacf(resid(moderate), main = "Partial Autocorrelation of residuals from single variable regression")

grid()函数向绘制的数据添加网格：

> grid()

结果如下：

步骤 4 - 改进模型

根据所有七个解释变量预测每日卡路里。使用拟合模型在多个不同 alpha 值下对样本进行拟合，使用拟合模型从原始样本中预测未在重新抽样的样本中的外袋点。这是通过选择适当的 alpha 值来在岭回归和 lasso 估计的极端之间创建平衡。

通过标准化创建矩阵X。as.matrix()函数将fitbit[ , -1]（即除日期列之外的部分）转换为矩阵：

 > X <- as.matrix(fitbit[ , -1])

打印X数据框。"head()"函数返回X数据框的第一部分。"X"数据框作为输入参数传递：

> head(X)

结果如下：

通过标准化创建向量Y：

> Y <- fitbit$Calories.Burned

打印Y数据框：

> Y

结果如下：

> set.seed(123)

生成常规序列：

 > alphas <- seq(from = 0, to  = 1, length.out = 10)
 > res <- matrix(0, nrow = length(alphas), ncol = 6)

创建每个 CV 运行的五次重复：

    > for(i in 1:length(alphas)){
    + for(j in 2:6){
    # k-fold cross-validation for glmnet
    + cvmod <- cv.glmnet(X, Y, alpha = alphas[i])
    + res[i, c(1, j)] <- c(alphas[i], sqrt(min(cvmod$cvm)))
    + }
    + }

创建要使用的数据集。"data.frame()"函数用于根据紧密耦合的变量集创建数据框。这些变量共享矩阵的性质：

> res <- data.frame(res)

打印res数据框：

> res

结果如下：

创建average_rmse向量：

> res$average_rmse <- apply(res[ , 2:6], 1, mean)

打印res$average_rmse向量：

> res$average_rmse

结果如下：

将res$average_rmse按升序排列。结果存储在res数据框中：

> res <- res[order(res$average_rmse), ]

打印res数据框：

> res

结果如下：

    > names(res)[1] <- "alpha"
    > res %>%
    + select(-average_rmse) %>%
    + gather(trial, rmse, -alpha) %>%
    + ggplot(aes(x = alpha, y = rmse)) +
    + geom_point() +
    + geom_smooth(se = FALSE) +
    + labs(y = "Root Mean Square Error") +
    + ggtitle("Cross Validation best RMSE for differing values of alpha")

结果如下：

> bestalpha <- res[1, 1]

打印bestalpha数据框：

> bestalpha

使用弹性网络比较普通最小二乘等价物与八个系数（七个解释变量加一个截距）的估计值。

确定 alpha 的最佳值处的 lambda。通过调用cv.glmnet()函数计算glmnet的 k 折交叉验证：

> crossvalidated <- cv.glmnet(X, Y, alpha = bestalpha)

创建模型。glmnet()通过惩罚最大似然估计拟合广义线性模型。正则化路径在正则化参数 lambda 的值网格上计算，对于 lasso 或elasticnet惩罚。X是输入矩阵，而Y是响应变量。alpha是elasticnet混合参数，范围是 0 ≤ α ≤ 1：

> moderate1 <- glmnet(X, Y, alpha = bestalpha)

建立普通最小二乘估计，以fitbit作为唯一的解释变量，以Calories.Burned作为响应变量。使用lm()函数拟合线性模型。Calories.Burned ~ Steps是公式，而fitbit是数据框。结果存储在OLSmodel数据框中：

> OLSmodel <- lm(Calories.Burned ~ ., data = fitbit)

打印OLSmodel数据框：

> OLSmodel

结果如下：

比较普通最小二乘等价物与八个系数（七个解释变量加一个截距）的估计值。结果存储在coeffs数据框中：

 > coeffs <- data.frame(original = coef(OLSmodel), 
 + shrunk = as.vector(coef(moderate1, s = crossvalidated$lambda.min)),
 + very.shrunk = as.vector(coef(moderate1, s = crossvalidated$lambda.1se)))

打印coeffs数据框：

> coeffs

结果如下：

将moderate数据框的值四舍五入到三位有效数字：

> round(coeffs, 3)

结果如下：

创建模型。glmnet()通过惩罚最大似然估计拟合广义线性模型：

> moderate2 <- glmnet(X, Y, lambda = 0)

打印moderate2数据框：

> moderate2

结果如下：

将值四舍五入到三位有效数字：

> round(data.frame("elastic, lambda = 0" = as.vector(coef(moderate2)), "lm" = coef(OLSmodel), check.names = FALSE), 3)

结果如下：

创建模型。在消除距离列后，glmnet()通过惩罚最大似然估计拟合广义线性模型：

> moderate3 <- glmnet(X[ , -2], Y, lambda = 0)

打印moderate3数据框：

> moderate3

结果如下：

建立普通最小二乘估计Y ~ X[ , -2]是公式。结果存储在moderate4数据框中：

> moderate4 <- lm(Y ~ X[ , -2])

打印moderate4数据框：

> moderate4

结果如下：

将数值四舍五入到三位有效数字：

> round(data.frame("elastic, lambda = 0" = as.vector(coef(moderate3)), "lm" = coef(moderate4), check.names = FALSE), 3)

结果如下：

步骤 5 - 比较模型

通过使用 bootstrapping 比较不同模型的预测能力，其中建模方法应用于数据的 bootstrap 重采样。然后使用估计模型来预测完整的原始数据集。

传递给 boot 以进行弹性建模的函数：

    > modellingfucn1 <- function(data, i){
    + X <- as.matrix(data[i , -1])
    + Y <- data[i , 1]
    # k-fold cross-validation for glmnet
    + crossvalidated <- cv.glmnet(X, Y, alpha = 1, nfolds = 30)
    # Fitting a generalized linear model via penalized maximum likelihood
    + moderate1 <- glmnet(X, Y, alpha = 1)
    # Computing the root mean squared error
    + rmse(predict(moderate1, newx = as.matrix(data[ , -1]), s =     crossvalidated$lambda.min), data[ , 1])
    + }

生成应用于数据的统计量的 R bootstrap 副本。fitbit是数据集，statistic = modellingfucn1是函数，当应用于fitbit时，返回包含感兴趣统计量的向量。R = 99表示 bootstrap 副本的数量：

> elastic_boot <- boot(fitbit, statistic = modellingfucn1, R = 99)

打印elastic_boot数据框：

 > elastic_boot

结果如下：

传递给 boot 以进行 OLS 建模的函数：

    > modellingOLS <- function(data, i){
    + mod0 <- lm(Calories.Burned ~ Steps, data = data[i, ])
    + rmse(predict(moderate, newdata = data), data[ , 1])
    + }

生成应用于数据的统计量的 R bootstrap 副本。fitbit是数据集，statistic = modellingOLS是函数，当应用于fitbit时，返回包含感兴趣统计量的向量。R = 99表示 bootstrap 副本的数量：

> lmOLS_boot <- boot(fitbit, statistic = modellingOLS, R = 99)

打印lmOLS_boot数据框：

> lmOLS_boot

结果如下：

生成应用于数据的统计量的 R bootstrap 副本。fitbit是数据集，statistic = modellingfucn2是函数，当应用于fitbit时，返回包含感兴趣统计量的向量。R = 99表示 bootstrap 副本的数量：

> lm_boot <- boot(fitbit, statistic = modellingfucn2, R = 99)

打印lm_boot数据框：

> lm_boot

结果如下：

 > round(c("elastic modelling" = mean(elastic_boot$t), 
 + "OLS modelling" = mean(lm_boot$t),
 + "OLS modelling, only one explanatory variable" = mean(lmOLS_boot$t)), 1)

结果如下：

使用缩放变量重新拟合模型。

创建模型。glmnet()通过惩罚最大似然估计拟合广义线性模型。

 > ordering <- c(7,5,6,2,1,3,4)
 > par(mar = c(5.1, 4.1, 6.5, 1), bg = "grey90")
 > model_scaled <- glmnet(scale(X), Y, alpha = bestalpha)
 > the_palette <- brewer.pal(7, "Set1")
 > plot(model_scaled, xvar = "dev", label = TRUE, col = the_pallete, lwd = 2, main = "Increasing contribution of different explanatory variablesnas penalty for including them is relaxed")
 > legend("topleft", legend = colnames(X)[ordering], text.col = the_palette[ordering], lwd = 2, bty = "n", col = the_palette[ordering])

结果如下：

维度缩减方法 - Delta 的机队

航空公司战略规划过程中的一个部分是机队规划。机队是指航空公司运营的飞机总数，以及构成总机队的具体飞机类型。飞机采购的航空公司选择标准基于技术/性能特性、经济和财务影响、环境法规和限制、营销考虑以及政治现实。机队构成是航空公司公司的关键长期战略决策。每种飞机类型都有不同的技术性能特性，例如，携带有效载荷在最大飞行距离或范围内的能力。它影响财务状况、运营成本，尤其是服务特定路线的能力。

准备中

为了进行降维，我们将使用收集于 Delta 航空公司机队的数据集。

第一步 - 收集和描述数据

应使用标题为 delta.csv 的数据集。该数据集采用标准格式。共有 44 行数据，34 个变量。

如何做到这一点...

让我们深入了解细节。

步骤 2 - 探索数据

第一步是加载以下包：

    > install.packages("rgl")
    > install.packages("RColorBrewer")
    > install.packages("scales")
    > library(rgl)
    > library(RColorBrewer)
    > library(scales)

注意

版本信息：本页面的代码在 R 版本 3.3.2（2016-10-31）上进行了测试

让我们探索数据并了解变量之间的关系。我们将首先导入名为 delta.csv 的 csv 数据文件。我们将把数据保存到 delta 数据框中：

 > delta <- read.csv(file="d:/delta.csv", header=T, sep=",", row.names=1)

探索 delta 数据框的内部结构。str() 函数显示数据框的内部结构。将作为 R 对象传递给 str() 函数的详细信息：

> str(delta)

结果如下：

探索与飞机物理特性相关的中间定量变量：住宿、巡航速度、航程、引擎、翼展、尾高和 Length.Scatter 折线图矩阵。plot() 函数是一个用于绘制 R 对象的通用函数。将 delta[,16:22] 数据框作为函数值传递：

> plot(delta[,16:22], main = "Aircraft Physical Characteristics", col = "red")

结果如下：

所有这些变量之间都存在正相关关系，因为它们都与飞机的整体尺寸有关。

步骤 3 - 应用主成分分析

可视化高维数据集，如引擎数量。对数据进行主成分分析。princomp() 函数对 delta 数据矩阵执行主成分分析。结果是 principal_comp_analysis，它是一个 princomp 类的对象：

> principal_comp_analysis <- princomp(delta)

打印 principal_comp_analysis 数据框：

> principal_comp_analysis

结果如下：

绘制 principal_comp_analysis 数据：

> plot(principal_comp_analysis, main ="Principal Components Analysis of Raw Data", col ="blue")

结果如下：

可以证明，第一个主成分具有标准差，它解释了数据中超过 99.8% 的方差。

打印主成分分析加载项。loadings() 函数使用 principal_comp_analysis 主成分分析数据对象作为输入：

> loadings(principal_comp_analysis)

结果如下：

观察加载的第一个列，很明显，第一个主成分仅仅是航程，以英里为单位。数据集中每个变量的尺度都不同。

在常规尺度上绘制方差。barplot() 绘制垂直和水平条形图。sapply() 是一个包装函数，它返回与 delta.horiz=T 相同长度的列表，表示条形图将水平绘制，第一个在底部：

 > mar <- par()$mar
 > par(mar=mar+c(0,5,0,0))
 > barplot(sapply(delta, var), horiz=T, las=1, cex.names=0.8, main = "Regular Scaling of Variance", col = "Red", xlab = "Variance")

结果如下：

在对数尺度上绘制方差。barplot() 绘制垂直和水平条形：

> barplot(sapply(delta, var), horiz=T, las=1, cex.names=0.8, log='x', main = "Logarithmic  Scaling of Variance", col = "Blue", xlab = "Variance")

结果如下：

> par(mar=mar)

第 4 步 - 缩放数据

在某些情况下，缩放 delta 数据是有用的，因为变量跨越不同的范围。scale() 函数作为函数对 delta 矩阵的列进行中心化和/或缩放。结果存储在 delta2 数据框中：

> delta2 <- data.frame(scale(delta))

验证方差是否均匀：

> plot(sapply(delta2, var), main = "Variances Across Different Variables", ylab = "Variances")

结果如下：

现在方差在变量间是恒定的。

将主成分应用于缩放后的数据 delta2。princomp() 函数对 delta2 数据矩阵执行主成分分析。结果是 principal_comp_analysis，它是一个 princomp 类的对象：

> principal_comp_analysis <- princomp(delta2)

绘制 principal_comp_analysis 对象：

> plot(principal_comp_analysis, main ="Principal Components Analysis of Scaled Data", col ="red")

结果如下：

> plot(principal_comp_analysis, type='l', main ="Principal Components Analysis of Scaled Data")

结果如下：

使用 summary() 函数生成各种模型拟合函数结果的摘要：

> summary(principal_comp_analysis)

结果如下：

将主成分应用于缩放后的数据 delta2。prcomp() 函数对 delta2 数据矩阵执行主成分分析。结果是 principal_comp_analysis，它是一个 prcomp 类的对象：

> principal_comp_vectors <- prcomp(delta2)

创建 principal_comp_vectors 的数据框：

> comp <- data.frame(principal_comp_vectors$x[,1:4])

使用 k = 4 进行 k 均值聚类。kmeans() 函数对 comp 执行 k 均值聚类。nstart=25 表示要选择的随机集的数量。iter.max=1000 是允许的最大迭代次数：

> k_means <- kmeans(comp, 4, nstart=25, iter.max=1000)

创建一个包含九种连续颜色的向量：

> palette(alpha(brewer.pal(9,'Set1'), 0.5))

绘制 comp：

> plot(comp, col=k_means$clust, pch=16)

结果如下：

第 5 步 - 在 3D 图中可视化

在 3D 中绘制 comp$PC1、comp$PC2、comp$PC3：

> plot3d(comp$PC1, comp$PC2, comp$PC3, col=k_means$clust) 

结果如下：

在 3D 中绘制 comp$PC1、comp$PC3、comp$PC4：

> plot3d(comp$PC1, comp$PC3, comp$PC4, col=k_means$clust)

结果如下：

按照大小顺序检查簇：

> sort(table(k_means$clust))

结果如下：

> clust <- names(sort(table(k_means$clust)))

如第一簇中显示的名称：

> row.names(delta[k_means$clust==clust[1],])

结果如下：

如第二簇中显示的名称：

> row.names(delta[k_means$clust==clust[2],])

结果如下：

如第三簇中显示的名称：

> row.names(delta[k_means$clust==clust[3],])

结果如下：

如第四簇中显示的名称：

> row.names(delta[k_means$clust==clust[4],])

结果如下：

主成分分析 - 理解世界美食

食物是我们身份的强大象征。有许多类型的食物识别，如民族、宗教和阶级识别。在存在味觉外国人（如出国或在外国人访问家乡时）的情况下，民族食物偏好成为身份标志。

准备工作

为了进行主成分分析，我们将使用在 Epicurious 菜谱数据集上收集的数据集。

第 1 步 - 收集和描述数据

将使用标题为 epic_recipes.txt 的数据集。该数据集为标准格式。

如何操作...

让我们深入了解细节。

第 2 步 - 探索数据

第一步是加载以下包：

> install.packages("glmnet") 
    > library(ggplot2)
    > library(glmnet)

注意

版本信息：本页代码在 R 版本 3.3.2（2016-10-31）上进行了测试

让我们探索数据并了解变量之间的关系。我们将首先导入名为 epic_recipes.txt 的 TXT 数据文件。我们将数据保存到 datafile 数据框中：

> datafile <- file.path("d:","epic_recipes.txt")

从表格格式文件中读取文件并创建数据框。datafile 是文件名，作为输入传递：

> recipes_data <- read.table(datafile, fill=TRUE, col.names=1:max(count.fields(datafile)), na.strings=c("", "NA"), stringsAsFactors = FALSE)

第 3 步 - 准备数据

将数据拆分为子集。aggregate() 函数将 recipes_data[,-1] 拆分并计算汇总统计信息。recipes_data[,-1] 是一个分组元素列表，每个元素与数据框中的变量长度相同。结果存储在 agg 数据框中：

> agg <- aggregate(recipes_data[,-1], by=list(recipes_data[,1]), paste, collapse=",")

创建一个向量、数组或值列表：

> agg$combined <- apply(agg[,2:ncol(agg)], 1, paste, collapse=",")

替换所有模式出现。gsub() 函数在搜索 agg$combined 后将每个 ,NA 替换为 ""：

> agg$combined <- gsub(",NA","",agg$combined)

提取所有菜系的名称：

> cuisines <- as.data.frame(table(recipes_data[,1]))

打印菜系数据框：

> cuisines

结果如下：

提取成分的频率：

 > ingredients_freq <- lapply(lapply(strsplit(a$combined,","), table), as.data.frame) 
 > names(ingredients_freq) <- agg[,1]

标准化成分的频率：

 > proportion <- lapply(seq_along(ingredients_freq), function(i) {
 + colnames(ingredients_freq[[i]])[2] <- names(ingredients_freq)[i]
 + ingredients_freq[[i]][,2] <- ingredients_freq[[i]][,2]/cuisines[i,2] 
 + ingredients_freq[[i]]}
 + )

包含 26 个元素，每个元素对应一种菜系：

    > names(proportion) <- a[,1]
    > final <- Reduce(function(...) merge(..., all=TRUE, by="Var1"), proportion)
    > row.names(final) <- final[,1]
    > final <- final[,-1]
    > final[is.na(final)] <- 0
    > prop_matrix <- t(final)
    > s <- sort(apply(prop_matrix, 2, sd), decreasing=TRUE)

scale() 函数将 prop_matrix 矩阵的列进行居中和/或缩放。结果存储在 final_impdata 数据框中：

 > final_imp <- scale(subset(prop_matrix, select=names(which(s > 0.1))))

创建热图。final_imp 是作为输入传递的数据框。trace="none" 表示字符串，指示是否在行或列上绘制实线 "trace"，"both" 或 "none"。key=TRUE 值表示应显示颜色键：

> heatmap.2(final_imp, trace="none", margins = c(6,11), col=topo.colors(7), key=TRUE, key.title=NA, keysize=1.2, density.info="none")

结果如下：

第 4 步 - 应用主成分分析

对数据进行主成分分析。princomp() 函数对 final_imp 数据矩阵执行主成分分析。结果是 pca_computation，它是一个 princomp 类的对象：

> pca_computation <- princomp(final_imp) 

打印 pca_computation 数据框：

> pca_computation

结果如下：

生成双变量图。pca_computation 是一个 princomp 类的对象。pc.biplot=TRUE 表示它是一个主成分双变量图：

> biplot(pca_computation, pc.biplot=TRUE, col=c("black","red"), cex=c(0.9,0.8), xlim=c(-2.5,2.5), xlab="PC1, 39.7%", ylab="PC2, 24.5%")

结果如下：

第五章。非线性

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

• 广义加性模型 - 测量新西兰的家庭收入

• 光滑样条 - 理解汽车和速度

• 局部回归 - 理解干旱预警和影响

广义加性模型 - 测量新西兰的家庭收入

收入调查提供了人们和家庭收入水平的快照。它给出了来自大多数来源的中位数和平均每周收入。存在不同人口群体之间的收入比较。收入是间歇性收到的，而消费是随时间平滑的。因此，可以合理地预期，消费与当前生活水平比当前收入更直接相关，至少在短期参考期内是这样。

准备工作

为了执行收缩方法，我们将使用 2013 年新西兰人口普查收集的数据集。

步骤 1 - 收集和描述数据

nzcensus 包包含超过 60 个新西兰人口统计值。这些值在网格块、面积单位、地区当局和地区议会级别累积。

如何做到这一点...

让我们深入了解。

步骤 2 - 探索数据

第一步是加载以下包：

 > devtools::install_github("ellisp/nzelect/pkg2")
> library(leaflet)
> library(nzcensus)
> library(Metrics)
> library(ggplot2)
> library(scales)
> library(boot)
> library(dplyr)
> library(Hmisc)
> library(mgcv)
> library(caret)
> library(grid)
> library(stringr)
> library(ggrepel)
> library(glmnet)
> library(maps)

从数据集中删除查塔姆群岛。AreaUnits2013 是一个 esriGeometryPolygon 几何类型对象。它定义了来自 2013 年人口普查模式的面积单位：

    > tmp <- AreaUnits2013[AreaUnits2013$WGS84Longitude> 0 & !is.na(AreaUnits2013$MedianIncome2013), ]

创建颜色调色板函数：

    > palette <- colorQuantile("RdBu", NULL, n = 10)

为弹出窗口创建标签。paste0() 函数在转换为字符后将向量连接起来：

    > labels <- paste0(tmp$AU_NAM, " $", format(tmp$MedianIncome2013, big.mark = ","))

绘制地图：

> leaflet() %>%
+ addProviderTiles("CartoDB.Positron") %>%
+ addCircles(lng = tmp$WGS84Longitude, lat = tmp$WGS84Latitude,
+ color = pal(-tmp$MedianIncome2013),
+ popup = labs,
+ radius = 500) %>%
+ addLegend(
+ pal = pal,
+ values = -tmp$MedianIncome2013,
+ title = "Quantile of median<br>household income",
+ position = "topleft",
+ bins = 5) 

结果如下：

步骤 3 - 为模型设置数据

将数据整理成方便的形状。消除区域的代码和名称，以及冗余的坐标系：

    > au <- AreaUnits2013 %>%     +  select(-AU2014, -AU_NAM, -NZTM2000Easting, -NZTM2000Northing) %>%     +  select(-PropWorked40_49hours2013, -Prop35to39_2013, -PropFemale2013)     > row.names(au) <- AreaUnits2013$AU_NAM

替换所有重复模式的实例。gsub() 函数搜索模式 "_2013"、"2013" 和 "Prop"，然后将它们替换为 names(au)：

 names(au) <- gsub("_2013", "", names(au))
> names(au) <- gsub("2013", "", names(au))
> names(au) <- gsub("Prop", "", names(au))

获取一个逻辑向量，指示一组案例是否完整：

    > au <- au[complete.cases(au), ]

提供一个通用名称：

    > data_use <- au

探索 data_use 数据框的维度。dim() 函数返回 data_use 框架的维度。将 data_use 数据框作为输入参数传递。结果清楚地表明有 1785 行数据和 69 列：

    > dim(data_use)

结果如下：

    > data_use <- data_use[the_data$WGS84Longitude > 100, ]

从字符向量创建语法有效的名称并设置它们。names() 函数在从返回的字符向量创建语法有效的名称的同时设置 data_use 对象的名称：

    > names(data_use) <- make.names(names(data_use))

显示从 data_use 数据框创建的名称：

    > names(data_use)

结果如下：

步骤 4 - 构建模型

估计非参数模型的强度。spearman2()计算 Spearman 的 rho 秩相关系数的平方，以及其推广，其中x可以非单调地与y相关。这是通过计算*(rank(x), rank(x)²)与y*之间的 Spearman 多重 rho 平方来完成的：

    > reg_data <- spearman2(MedianIncome ~ ., data = data_use)

按降序排列数据：

    > reg_data[order(-reg_data[ ,6])[1:15], ]

结果如下：

将灵活的样条分配给前 15 个变量。terms()函数从多个 R 数据对象中提取terms对象：

> reg_formula <- terms(MedianIncome ~
s(FullTimeEmployed, k = 6) +
s(InternetHH, k = 6) +
s(NoQualification, k = 5) +
s(UnemploymentBenefit, k = 5) +
s(Smoker, k = 5) +
s(Partnered, k = 5) +
s(Managers, k = 4) +
s(Bachelor, k = 4) +
s(SelfEmployed, k = 4) +
s(NoMotorVehicle, k = 4) +
s(Unemployed, k = 3) +
s(Labourers, k = 3) +
s(Worked50_59hours, k = 3) +
s(Separated, k = 3) +
s(Maori, k = 3) +
s(WGS84Longitude, WGS84Latitude) +
.,
data = data_use)

拟合广义加性模型。reg_formula是公式，而data_use是数据集。

> gam_model <- gam(reg_formula, data = data_use) 

绘制gam_model。

    > par(bty = "l", mar = c(5,4, 2, 1))     > par(mar = rep(2, 4))     > plot(gam_model, residuals = TRUE, pages = 1, shade = TRUE, seWithMean = TRUE, ylab = "")

结果如下：

    > rmses_gam_boot <- boot(data = data_use, statistic = fit_gam, R = 99)

打印rmses_gam_boot数据框：

    > rmses_gam_boot

结果如下：

计算rmses_gam_boot$t的均值：

    > gam_rmse <- mean(rmses_gam_boot$t)

打印gam_rmse数据框：

    > gam_rmse

结果如下：

平滑样条 - 理解汽车和速度

为了确定用于拟合模型的统计参数，可以使用多种方法。在每种情况下，拟合都涉及从数据中估计少量参数。除了估计参数外，两个重要阶段是确定合适的模型和验证模型。这些平滑方法可以以多种方式使用：帮助理解并生成平滑图，从平滑数据形状中识别合适的参数模型，或者专注于感兴趣的效应，以消除无用的复杂效应。

如何做到这一点...

让我们深入了解细节。

步骤 1 - 探索数据

第一步是加载以下包：

> install.packages("graphics")
> install.packages("splines")
> library(graphics)
> library(splines)

创建矩阵。cbind()函数将数字序列组合成一个矩阵。然后将结果传递给matrix()函数，创建一个两行的矩阵。结果存储在矩阵中：

    > matrx = matrix(cbind(1,.99, .99,1),nrow=2)

步骤 2 - 创建模型

Cholesky 分解创建正定矩阵A，它可以分解为A=LL^T，其中L是下三角矩阵，对角线元素为正。chol()函数计算实对称正定方阵的 Cholesky 分解。结果存储在cholsky中：

> cholsky = t(chol(matrx))
> nvars = dim(cholsky)[1]

密度分布的观测数：

    > numobs = 1000     
> set.seed(1)

使用正态分布计算矩阵。rnorm()计算正态分布，numobs为要使用的观测数。然后将结果用于matrix()函数计算矩阵，nrow=nvars为两行，ncol=numobs为 1,000 列。结果存储在random_normal中：

    > random_normal = matrix(rnorm(nvars*numobs,10,1), nrow=nvars, ncol=numobs)

执行矩阵乘法。cholsky与矩阵random_normal相乘：

    > X = cholsky %*% random_normal

转置矩阵X：

    > newX = t(X)

创建矩阵的数据框。as.data.frame() 函数创建数据框 raw，它是一组紧密耦合的变量，这些变量共享矩阵 newX 的许多属性：

    > raw = as.data.frame(newX)

打印原始数据框。head() 函数返回原始数据框的第一部分。原始数据框作为输入参数传递：

    > head(raw)

结果如下：

创建 random_normal 的转置数据框。t() 函数创建 random_normal 矩阵的转置矩阵，然后将其转换为紧密耦合的变量集合。这些变量共享矩阵的许多属性：

    > raw_original = as.data.frame(t(random_normal))

结合响应和 predictor1 的名称。c() 函数将响应和 predictor1 作为参数结合，形成一个向量：

    > names(raw) = c("response","predictor1")

raw$predictor1 的指数增长到 3 次方：

    > raw$predictor1_3 = raw$predictor1³

打印 raw$predictor1_3 数据框。head() 函数返回 raw$predictor1_3 数据框的第一部分。raw$predictor1_3 数据框作为输入参数传递：

    > head(raw$predictor1_3)

结果如下：

raw$predictor1 的指数增长到 2 次方：

    > raw$predictor1_2 = raw$predictor1²

打印 raw$predictor1_2 数据框。head() 函数返回 raw$predictor1_2 数据框的第一部分。raw$predictor1_2 数据框作为输入参数传递：

    > head(raw$predictor1_2)

结果如下：

使用 raw$response ~ raw$predictor1_3 作为公式的普通最小二乘估计。lm() 函数用于拟合线性模型。raw$response ~ raw$predictor1_3 是公式。结果存储在拟合数据框中：

    > fit = lm(raw$response ~ raw$predictor1_3)

打印拟合数据框：

    > fit

结果如下：

绘制普通最小二乘估计公式。plot() 函数是用于绘制 R 对象的通用函数。raw$response ~ raw$predictor1_3 公式作为函数值传递：

    > plot(raw$response ~ raw$predictor1_3, pch=16, cex=.4, xlab="Predictor", ylab="Response", col ="red", main="Simulated Data with Slight Curve")

结果如下：

在当前图表上添加直线函数：

    > abline(fit)

结果如下：

在 x 轴上拟合汽车和速度的值：

    > x_axis <- with(cars, speed)

在 y 轴上拟合汽车和速度的值：

    > y_axis <- with(cars, dist)

设置平滑曲线评估的点数：

    > eval_length = 50

第 3 步 - 拟合平滑曲线模型

在两个变量之间拟合平滑曲线是一种非参数方法，因为传统回归方法的线性假设已经被放宽。它被称为局部回归，因为在点 x 处的拟合是加权向 x 附近的数据：

loess.smooth()函数在散点图上绘制并添加计算出的平滑曲线。x_axis，y_axis是提供给绘图 x 和 y 坐标的参数。例如evaluation = eval.length，eval_length = 50表示平滑曲线评估的点数。span=.75是平滑度参数。degree=1是局部多项式的次数：

> fit_loess <- loess.smooth(x_axis, y_axis, evaluation = eval_length, family="gaussian", span=.75, degree=1) 

打印fit_loess数据框：

    > fit_loess

结果如下：

在一个或多个数值预测的基础上，使用局部拟合拟合多项式表面。loess()函数拟合多项式表面。y_axis ~ x_axis表示公式。span=.75是平滑度参数。degree=1是局部多项式的次数：

    > fit_loess_2 <- loess(y_axis ~ x_axis, family="gaussian", span=.75, degree=1)

打印fit_loess_2数据框：

    > fit_loess_2

结果如下：

生成y轴最小值和最大值的常规序列。Seq()函数接受length.out=eval_length，例如eval_length = 50，表示从x轴的最小值和最大值生成的序列的期望长度：

    > new_x_axis = seq(min(x_axis),max(x_axis), length.out=eval_length)

打印new_x_axis数据框：

    > new_x_axis

结果如下：

在fit.loess模型上设置 95%的置信水平：

> conf_int = cbind( 
 + predict(fit_loess_2, data.frame(x=new_x_axis)), 
 + predict(fit_loess_2, data.frame(x=new_x_axis))+ 
 + predict(fit_loess_2, data.frame(x=new_x_axis), se=TRUE)$se.fit*qnorm(1-.05/2), 
 + predict(fit_loess_2, data.frame(x=new_x_axis))- 
+ predict(fit_loess_2, data.frame(x=new_x_axis), se=TRUE)$se.fit*qnorm(1-.05/2) 
 + )

使用y_axis ~ x_axis作为公式构建普通最小二乘估计。使用lm()函数拟合线性模型。y_axis ~ x_axis是公式。结果存储在fit_lm数据框中：

    > fit_lm = lm(y_axis ~ x_axis)

打印fit_lm数据框：

    > fit_lm

结果如下：

构建多项式函数。y_axis ~ poly(x_axis,3)是具有三个自由度的多项式函数。使用lm()函数拟合线性模型。y_axis ~ poly(x_axis,3)是公式。结果存储在fit_poly数据框中：

    > fit_poly = lm(y_axis ~ poly(x_axis,3) )

打印fit_poly数据框：

    > fit_poly

结果如下：

构建自然样条函数。y_axis ~ ns(x_axis, 3)是具有 3 个自由度的自然样条函数。使用lm()函数拟合线性模型。y_axis ~ ns(x_axis, 3)是公式。结果存储在fit_nat_spline数据框中：

    > fit_nat_spline = lm(y_axis ~ ns(x_axis, 3) )

打印fit_nat_spline数据框：

    > fit_nat_spline

结果如下：

样条曲线的平滑处理：

    > fit_smth_spline <- smooth.spline(y_axis ~ x_axis, nknots=15)

打印fit_smth_spline数据框：

    > fit_smth_spline

结果如下：

步骤 4 - 绘制结果

绘制模型：

    > plot(x_axis, y_axis, xlim=c(min(x_axis),max(x_axis)), ylim=c(min(y_axis),max(y_axis)), pch=16, cex=.5, ylab = "Stopping Distance (feet)", xlab= "Speed (MPH)", main="Comparison of Models", sub="Splines")

结果如下：

向图中添加额外的模型。绘制带有置信区间的 LOESS：

    > matplot(new_x_axis, conf_int, lty = c(1,2,2), col=c(1,2,2), type = "l", add=T)

结果如下：

绘制普通最小二乘估计。predict() 函数根据线性模型预测值。fit_lm 是 lm 类的对象：

    > lines(new_x_axis, predict(fit_lm, data.frame(x=new_x_axis)), col="red", lty=3)

结果如下：

绘制多项式函数估计：

    > lines(new_x_axis, predict(fit_poly, data.frame(x=new_x_axis)), col="blue", lty=4)

结果如下：

绘制自然样条函数：

    > lines(new_x_axis, predict(fit_nat_spline, data.frame(x=new_x_axis)), col="green", lty=5)

结果如下：

绘制平滑样条：

    > lines(fit_smth_spline, col="dark grey", lty=6)

结果如下：

绘制核曲线。ksmooth() 函数：

    > lines(ksmooth(x_axis, y_axis, "normal", bandwidth = 5), col = 'purple', lty=7)

结果如下：

局部回归 - 理解干旱预警和影响

干旱是一种自然灾害，其特征是低于预期的或低于正常水平的降雨。当这种状况在超过正常时间周期的更长时期内持续时，不足以满足人类活动的需求，对环境有害。干旱是一种暂时现象。干旱的三个主要特征是强度、持续时间和空间覆盖范围。干旱预警系统可以帮助识别气候变化，了解供水趋势，并为即将到来的紧急情况做好准备。干旱预警可以帮助决策者采取适当的措施来应对即将到来的挑战。然后，他们可以衡量影响的严重程度，并了解特定地点、特定人群或经济部门脆弱性的潜在原因，以降低风险。

准备就绪

让我们从配方开始。

第 1 步 - 收集和描述数据

dataRetrieval 包是一组函数，用于帮助检索美国地质调查局（USGS）和美国环境保护署（EPA）。

如何操作...

让我们深入了解细节。

第 2 步 - 收集和探索数据

第一步是加载以下包：

> library(dataRetrieval)
> library(dplyr)

获取站点编号。站点编号通常是一个八位数的数字，它被表示为一个字符串或向量：

> siteNumber <- c("01538000") 

获取参数代码：

    > parameterCd <- "00060"

使用站点编号和参数代码从 NWIS 网络服务导入数据。结果存储在 Q_daily 数据框中：

    > Q_daily <- readNWISdv(siteNumber, parameterCd)

打印 Q_daily 数据框。tail() 函数返回 Q_daily 数据框的最后部分。Q_daily 数据框作为输入参数传递：

    > tail(Q_daily)

结果如下：

探索 Q_daily 数据框的内部结构。str() 函数显示数据框的内部结构。Q_daily 作为 R 对象传递给 str() 函数：

    > str(Q_daily)

结果如下：

重命名列 -- renameNWISColumns() 函数重命名从 NWIS 获取的列。Q_daily 是从 NWIS 网站获取的每日或单位值数据集：

    > Q_daily <- renameNWISColumns(Q_daily)

打印重命名的 Q_daily 数据框。tail() 函数返回 Q_daily 数据框的最后部分。Q_daily 数据框作为输入参数传递：

    > tail(Q_daily)

结果如下：

从 USGS 文件站点导入数据。readNWISsite() 函数使用 8 位数字的 siteNumber，它代表 USGS 站点编号。结果存储在 stationInfo 数据框中：

    > stationInfo <- readNWISsite(siteNumber)

步骤 3 - 计算移动平均

检查缺失天数：

> if(as.numeric(diff(range(Q_daily$Date))) != (nrow(Q_daily)+1)){
+ fullDates <- seq(from=min(Q_daily$Date),
+ to = max(Q_daily$Date), by="1 day")
+ fullDates <- data.frame(Date = fullDates,
+ agency_cd = Q_daily$agency_cd[1],
+ site_no = Q_daily$site_no[1],
+ stringsAsFactors = FALSE)
+ Q_daily <- full_join(Q_daily, fullDates,
+ by=c("Date","agency_cd","site_no")) %>%
+ arrange(Date)
+ }

计算 30 天的移动平均。filter() 函数对时间序列应用线性滤波。sides=1，仅对过去值应用滤波系数：

> moving_avg <- function(x,n=30){stats::filter(x,rep(1/n,n), sides=1)}     > 
Q_daily <- Q_daily %>% mutate(rollMean = as.numeric(moving_avg(Flow)), day.of.year = as.numeric(strftime(Date, format = "%j")))

打印 Q_daily 数据框。tail() 函数返回 Q_daily 数据框的最后部分。Q_daily 数据框作为输入参数传递：

    > tail(Q_daily)

结果如下：

步骤 4 - 计算百分位数

计算历史百分位数。使用相应的概率计算各种分位数。然后，使用 summarize() 函数将数据框折叠成单行。最后，使用 group_by() 函数，将结果（以表格形式）转换为并分组到表格中：

> Q_summary >- Q_daily %>%
+ group_by(day.of.year) %>%
+ summarize(p75 = quantile(rollMean, probs = .75, na.rm = TRUE),
+ p25 = quantile(rollMean, probs = .25, na.rm = TRUE),
+ p10 = quantile(rollMean, probs = 0.1, na.rm = TRUE),
+ p05 = quantile(rollMean, probs = 0.05, na.rm = TRUE),
+ p00 = quantile(rollMean, probs = 0, na.rm = TRUE))

从系统中获取当前年份：

> current_year <- as.numeric(strftime(Sys.Date(), format = "%Y"))
> summary.0 <- Q_summary %>% mutate(Date = as.Date(day.of.year - 1,
origin = paste0(current_year-2,"-01-01")), day.of.year = day.of.year - 365)
> summary.1 <- Q_summary %>% mutate(Date = as.Date(day.of.year - 1,
origin = paste0(current_year-1,"-01-01")))
> summary.2 <- Q_summary %>% mutate(Date = as.Date(day.of.year - 1,
origin = paste0(current_year,"-01-01")), day.of.year = day.of.year + 365)

合并每个数据框：

    > Q_summary <- bind_rows(summary.0, summary.1, summary.2) 

打印 Q_summary 数据框：

    > Q_summary

结果如下：

    > smooth.span <- 0.3

根据线性模型预测值并拟合多项式曲面。loess() 函数拟合多项式曲面。p75~day.of.year 表示公式，而 span = smooth.span 例如 smooth.span= 0.3 控制平滑程度：

    > Q_summary$sm.75 <- predict(loess(p75~day.of.year, data = Q_summary, span = smooth.span))

打印 Q_summary$sm.75 数据框：

    > head(Q_summary$sm.75)

结果如下：

    > Q_summary$sm.25 <- predict(loess(p25~day.of.year, data = Q_summary, span = smooth.span))

打印 Q_summary$sm.25 数据框：

    > head(summaryQ$sm.25)

结果如下：

    > Q_summary$sm.10 <- predict(loess(p10~day.of.year, data = Q_summary, span = smooth.span))

打印 Q_summary$sm.10 数据框：

    > head(summaryQ$sm.10)

结果如下：

    > Q_summary$sm.05 <- predict(loess(p05~day.of.year, data = Q_summary, span = smooth.span))

打印 Q_summary$sm.05 数据框：

    > head(summaryQ$sm.05)

结果如下：

    > Q_summary$sm.00 <- predict(loess(p00~day.of.year, data = Q_summary, span = smooth.span))

打印 Q_summary$sm.05 数据框：

    > head(summaryQ$sm.00)

结果如下：

    > Q_summary <- select(Q_summary, Date, day.of.year, sm.75, sm.25, sm.10, sm.05, sm.00) %>% filter(Date >= as.Date(paste0(current_year-1,"-01-01")))

打印 Q_summary 数据框：

    > Q_summary

结果如下：

    > latest.years <- Q_daily %>% filter(Date >= as.Date(paste0(current_year-1,"-01-01"))) %>% mutate(day.of.year = 1:nrow(.))

步骤 5 - 绘制结果

绘制数据：

> title.text <- paste0(stationInfo$station_nm,"n", "Provisional Data - Subject to changen", "Record Start = ", min(Q_daily$Date), "  Number of years = ", as.integer (as.numeric(difftime(time1 = max(Q_daily$Date), time2 = min(Q_daily$Date), units = "weeks"))/52.25), "nDate of plot = ",Sys.Date(), "  Drainage Area = ",stationInfo$drain_area_va, "mi²")     > mid.month.days <- c(15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349)     > month.letters <- c("J","F","M","A","M","J","J","A","S","O","N","D")     > start.month.days <- c(1, 32, 61, 92, 121, 152, 182, 214, 245, 274, 305, 335)     > label.text <- c("Normal","DroughtWatch","DroughtWarning","Drought Emergency")     > year1_summary <- data.frame(Q_summary[2:366,])     > head(year1_summary) 

结果如下：

    > year2_summary <- data.frame(Q_summary[367:733,])     
> head(year2_summary)

结果如下：

> simple.plot <- ggplot(data = Q_summary, aes(x = day.of.year)) + 
+ geom_ribbon(aes(ymin = sm.25, ymax = sm.75, fill = "Normal")) + 
    + geom_ribbon(aes(ymin = sm.10, ymax = sm.25, fill =       "Drought Watch")) +
    + geom_ribbon(aes(ymin = sm.05, ymax = sm.10, fill = "Drought Warning")) +
+ geom_ribbon(aes(ymin = sm.00, ymax = sm.05, fill = "Drought Emergency")) + 
+ scale_y_log10(limits = c(1,1000)) + 
+ geom_line(data = latest.years, aes(x=day.of.year, y=rollMean, color = "30-Day Mean"),size=2) + 
+ geom_vline(xintercept = 365) 
    > simple.plot

结果如下：

第六章. 监督学习

在本章中，我们将涵盖以下内容：

• 决策树学习 - 胸痛患者的健康指导

• 决策树学习 - 房地产价值基于收入的分布

• 决策树学习 - 预测股票运动的方向

• 朴素贝叶斯 - 预测股票运动的方向

• 随机森林 - 货币交易策略

• 支持向量机 - 货币交易策略

• 随机梯度下降 - 成人收入

简介

决策树学习：决策树是分类和预测问题中非常流行的工具。决策树是一种递归地将实例空间或变量集进行划分的分类器。决策树以树结构表示，其中每个节点可以分类为叶节点或决策节点。叶节点包含目标属性的值，而决策节点指定对单个属性值要实施的规则。每个决策节点根据输入属性值的某个离散函数将实例空间划分为两个或更多子空间。每个测试考虑一个属性，因此实例空间根据属性值进行划分。在数值属性的情况下，条件指的是一个范围。在决策节点上实施规则后，子树是一个结果。每个叶节点都包含一个概率向量，表示目标属性具有某个值的概率。通过沿着路径的测试结果，从树的根节点导航到叶节点来对实例进行分类。

使用决策树挖掘数据的关键要求如下：

• 属性值描述：对象可以用一组固定的属性或属性来表示

• 预定义类别：要分配给示例的类别必须是监督数据

• 充足数据：使用多个训练案例

朴素贝叶斯：朴素贝叶斯是一种监督学习方法。它是一个线性分类器。它基于贝叶斯定理，该定理表明一个类别的特定特征的存在与任何其他特征的存在无关。它是一个健壮且高效的算法。贝叶斯分类器可以预测类成员概率，例如给定元组属于特定类的概率。贝叶斯信念网络是联合条件概率分布。它允许在变量子集之间定义类条件独立性。它提供了一个因果关系的图形模型，可以在其上进行学习。

随机森林：随机森林是决策树的集合，提供对数据结构的预测。它们是利用多个决策树在合理随机化、集成学习中的力量来产生预测模型的一种工具。它们为每个记录提供变量排名、缺失值、分割和报告，以确保深入理解数据。在每棵树构建完成后，所有数据都会通过树。对于每一对案例，计算邻近区域。如果两个案例占据相同的终端节点，它们的邻近区域增加一。运行结束后，通过树的数量进行归一化。邻近区域用于替换缺失数据、定位异常值和揭示数据的低维理解。训练数据，即袋外数据，用于估计分类错误和计算变量的重要性。

随机森林在大数据库上运行非常高效，产生准确的结果。它们处理多个变量而不删除，给出变量对解决分类问题重要性的估计。它们在森林构建过程中生成内部无偏估计的泛化误差。随机森林是估计缺失数据的有效方法，并且在大量数据缺失时保持准确性。

支持向量机：机器学习算法使用正确的特征集来解决学习问题。SVMs 利用一个（非线性）映射函数φ，将输入空间中的数据转换为特征空间中的数据，以便使问题线性可分。然后 SVM 发现最优的分离超平面，然后通过φ-1 将其映射回输入空间。在所有可能超平面中，我们选择距离最近数据点（边缘）距离尽可能大的那个超平面。

决策树学习 - 胸痛患者的健康指导文件

健康指导文件声明了关于个人在各种医疗条件下未来医疗保健的指示。它指导个人在紧急情况下或需要时做出正确的决定。该文件帮助个人了解其医疗保健决策的性质和后果，了解指导的性质和影响，自由自愿地做出这些决定，并以某种方式传达这些决定。

准备工作

为了执行决策树分类，我们将使用从心脏病患者数据集中收集的数据集。

第 1 步 - 收集和描述数据

将使用标题为Heart.csv的数据集，该数据集以 CSV 格式提供。数据集是标准格式。有 303 行数据。有 15 个变量。数值变量如下：

• Age

• Sex

• RestBP

• Chol

• Fbs

• RestECG

• MaxHR

• ExAng

• Oldpeak

• Slope

• Ca

非数值变量如下：

• ChestPain

• Thal

• AHD

如何做...

让我们深入了解。

第 2 步 - 探索数据

以下包需要在第一步执行时加载：

> install.packages("tree")
> install.packages("caret")
> install.packages("e1071")
> library(tree)
> library(caret)

注意

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）上进行了测试。

让我们探索数据并了解变量之间的关系。我们将首先导入名为 Heart.csv 的 CSV 数据文件。我们将数据保存到 AHD_data 数据框中：

    > AHD_data <- read.csv("d:/Heart.csv", header = TRUE)

探索 AHD_data 数据框的内部结构。str() 函数显示数据框的内部结构。AHD_data 作为 R 对象传递给 str() 函数：

> str(AHD_data) 

结果如下：

打印 AHD_data 数据框。head() 函数返回 AHD_data 数据框的前部分。AHD_data 数据框作为输入参数传递：

    > head(AHD_data)

结果如下：

探索 AHD_data 数据框的维度。dim() 函数返回 AHD_data 数据框的维度。将 AHD_data 数据框作为输入参数传递。结果清楚地表明有 303 行数据和 15 列：

    >dim(AHD_data)

结果如下：

第 3 步 - 准备数据

需要准备数据以执行模型构建和测试。数据分为两部分--一部分用于构建模型，另一部分用于测试模型，这将准备。

使用 createDataPartition() 函数创建数据的分割。将 AHD_data 作为参数传递给函数。进行随机抽样。表示用于训练的数据百分比的 p。在这里，p 的值为 0.5，这意味着 50% 的数据用于训练。List = 'FALSE' 避免以列表的形式返回数据。结果存储在数据框 split 中：

    > split <- createDataPartition(y=AHD_data$AHD, p = 0.5, list=FALSE)

调用 split 数据框显示用于训练目的的训练集数据：

    > split

结果如下：

将创建训练数据。使用 split 数据框创建训练数据。train 数据框用于存储训练数据的值：

    > train <- AHD_data[split,]

打印训练数据框：

    > train

结果如下：

将创建测试数据。使用 split 数据框创建测试数据。split 数据框前的 - 符号表示所有那些未被考虑用于训练目的的数据行。测试数据框用于存储测试数据的值：

    > test <- AHD_data[-split,]

打印测试数据框：

    > test

结果如下：

第 4 步 - 训练模型

模型现在将被准备并在训练数据集上训练。当数据集被分成组时使用决策树，与调查数值响应及其与一组描述符变量的关系相比。在 R 中使用 tree() 函数实现分类树。

使用 tree() 函数实现分类树。通过二分递归分割来生长树。训练数据集上的 AHD 字段用于形成分类树。结果数据框存储在 trees 数据框中：

    > trees <- tree(AHD ~., train)

将显示数据框的图形版本。plot() 函数是 R 对象绘图的通用函数。将数据框 trees 作为函数值传递：

    > plot(trees)

结果如下：

通过运行交叉验证实验来查找偏差或错误分类的数量。将使用 cv.tree() 函数。将 trees 数据框对象传递。FUN=prune.misclass 通过递归剪掉最不重要的分割来获取提供的 data frame trees 的嵌套子树序列。结果存储在 cv.trees 数据框中：

    > cv.trees <- cv.tree(trees, FUN=prune.misclass)

打印数据框 cv.trees 的结果：

    > cv.trees

$dev 字段给出了每个 K 的偏差。

结果如下：

使用 plot() 函数数据框，显示 cv.trees。$dev 值位于 y 轴（右侧）。$k 值位于顶部。$size 值位于 x 轴。

如清晰可见，当 $size = 1，$k = 30.000000，$dev = 1。我们使用以下方式绘制数据框：

    > plot(cv.trees)

结果如下：

步骤 5 - 改进模型

让我们通过分割偏差最低的树来改进模型。调用 prune.misclass() 函数来分割树。prune.misclass 通过递归剪掉最不重要的分割来获取提供的 data frame trees 的嵌套子树序列。结果存储在 prune.trees 数据框中。best=4 表示要返回的成本-复杂度序列中特定子树的大小（例如，终端节点的数量）：

    > prune.trees <- prune.misclass(trees, best=4)

使用 plot() 函数数据框，显示 prune.trees：

    > plot(prune.trees)

结果如下：

向前面的修剪树添加文本：

    > text(prune.trees, pretty=0)

结果如下：

为了根据线性模型对象预测值，我们将使用 predict() 函数。将 prune.trees 作为对象传递。将 test 数据对象传递作为查找预测变量的对象。结果将存储在 tree.pred 数据框中：

    > tree.pred <- predict(prune.trees, test, type='class')

显示变量 test.pred 的值：

    > tree.pred

结果如下：

总结模型的成果。confusionMatrix() 计算观察到的和预测的类别的交叉表。tree.pred 作为预测类别的因子传递：

    > confusionMatrix(tree.pred, test$AHD)

结果如下：

决策树学习 - 基于收入的房地产价值分布

收入一直是房地产作为一种资产类别提供的具有吸引力的长期总回报的一个基本组成部分。投资房地产产生的年度收入回报比股票高出 2.5 倍以上，仅落后于债券 50 个基点。房地产通常为租户支付的租金提供稳定的收入来源。

准备工作

为了执行决策树分类，我们将使用从房地产数据集中收集的数据集。

步骤 1 - 收集和描述数据

将使用标题为 RealEstate.txt 的数据集。此数据集以 TXT 格式提供，标题为 RealEstate.txt。数据集是标准格式。有 20,640 行数据。9 个数值变量如下：

• MedianHouseValue

• MedianIncome

• MedianHouseAge

• TotalRooms

• TotalBedrooms

• Population

• Households

• Latitude

• Longitude

如何做到这一点...

让我们深入了解细节。

步骤 2 - 探索数据

需要在第一步中加载以下包：

    > install.packages("tree")

注意

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）中进行了测试。

让我们探索数据并了解变量之间的关系。我们将从导入名为 RealEstate.txt 的 TXT 数据文件开始。我们将数据保存到 realEstate 数据框中：

    > realEstate <- read.table("d:/RealEstate.txt", header=TRUE)

探索 realEstate 数据框的维度。dim() 函数返回 realEstate 框的维度。realEstate 数据框作为输入参数传递。结果清楚地表明有 20,640 行数据和 9 列：

    > dim(realEstate)

结果如下：

探索 realEstate 数据框的内部结构。str() 函数显示数据框的内部结构。realEstate 作为 R 对象传递给 str() 函数：

    > str(realEstate)

结果如下：

打印 realEstate 数据框。head() 函数返回 realEstate 数据框的前部分。realEstate 数据框作为输入参数传递：

    > head(realEstate)

结果如下：

打印 realEstate 数据框的摘要。summary() 函数是一个多功能函数。summary() 是一个通用函数，它提供与单个对象或数据框相关的数据摘要。realEstate 数据框作为 R 对象传递给 summary() 函数：

    > summary(realEstate)

结果如下：

步骤 3 - 训练模型

模型现在将在数据集上准备。决策树是分类和预测的工具。它们代表人类可以理解并用于如数据库等知识系统的规则。它们通过从树的根开始并移动到叶节点来对实例进行分类。节点指定对单个属性的测试，叶节点指示目标属性的值，边分割出一个属性。

使用tree()函数实现分类树。通过二元递归分区来生长树。这些模型是计算密集型技术，因为它们根据响应变量与一个或多个预测变量的关系递归地将响应变量分割成子集。

公式表达式基于变量纬度和经度的总和。总和的结果存储在MedianHouseValue的对数值中。data=realEstate表示优先解释公式、权重和子集的数据框。

结果数据框存储在数据框treeModel中：

> treeModel <- tree(log(MedianHouseValue) ~ Longitude + Latitude, data=realEstate) 

将显示treeModel的摘要。摘要显示了所使用的公式，以及树中的终端节点或叶子的数量。还显示了残差的统计分布。

使用summary()函数显示treeModel的统计摘要。它是一个泛型，用于生成各种拟合函数的结果摘要。希望进行摘要的数据框是treeModel，它作为输入参数传递。

在这里，偏差表示均方误差：

    > summary(treeModel)

结果如下：

将显示treeModel数据框的图形版本。plot()函数是用于绘制 R 对象的泛型函数。treeModel数据框作为函数值传递：

> plot(treeModel) 

结果如下：

在显示treeModel数据框的图形版本后，需要插入文本以显示每个节点和叶子的值。使用text()函数在给定的坐标处插入标签向量中给出的字符串：

    > text(treeModel, cex=.75)

结果如下：

第 4 步 - 比较预测

将预测与反映全球价格趋势的数据集进行比较。我们希望总结MedianHouseValue的频率分布，以便于报告或比较。最直接的方法是使用分位数。分位数是分布中的点，与该分布中值的排名顺序相关。分位数将分割MedianHouseValue分布，使得观测值在分位数下方的比例是给定的。

quantile() 函数产生与给定概率相对应的样本分位数。realEstate$MedianHouseValue 是想要样本分位数的数值向量。quantile() 函数返回长度为的 priceDeciles 向量：

    > priceDeciles <- quantile(realEstate$MedianHouseValue, 0:10/10)

显示 priceDeciles 数据框的值：

    > priceDeciles

结果如下：

接下来，将显示 priceDeciles 的摘要。使用 summary() 函数显示 priceDeciles 的统计摘要。希望摘要的数据框是 priceDeciles，它作为输入参数传递：

    > summary(priceDeciles)

结果如下：

将 priceDeciles 向量划分为不同的范围。cut() 函数根据它们所属的区间来划分区间范围。realEstate 数据框中的数值向量 MedianHouseValue 需要通过切割转换为因子：

    > cutPrices <- cut(realEstate$MedianHouseValue, priceDeciles, include.lowest=TRUE)

打印 cutPrices 数据框。head() 函数返回 cutPrices 数据框的前部分。cutPrices 数据框作为输入参数传递：

    > head(cutPrices)

结果如下：

将显示 cutPrices 的摘要。使用 summary() 函数显示 treeModel 的统计摘要。希望摘要的数据框是 cutPrices，它作为输入参数传递：

    > summary(cutPrices)

结果如下：

绘制 cutPrices 的值。plot() 函数是 R 对象绘图的通用函数。realEstate 数据集中的经度变量代表图中点的 x 坐标。realEstate 数据集中的纬度变量代表图中点的 y 坐标。col=grey(10:2/11) 代表绘图颜色。pch=20 代表在绘图点时使用的符号大小。xlab="Longitude" 代表 x 轴的标题，而 ylab="Latitude" 代表 y 轴的标题：

> plot(realEstate$Longitude, realEstate$Latitude, col=grey(10:2/11)[cutPrices], pch=20, xlab="Longitude",ylab="Latitude") 

结果如下：

将显示 Longitude 的摘要。使用 summary() 函数显示统计摘要：

    > summary(realEstate$Longitude)

结果如下：

打印 Longitude 数据框。head() 函数返回 Longitude 数据框的前部分：

    > head(realEstate$Longitude)

结果如下：

将显示 Latitude 的摘要。使用 summary() 函数显示统计摘要：

    > summary(realEstate$Latitude)

结果如下：

打印 纬度 数据框。head() 函数返回 纬度 数据框的前部分：

    > head(realEstate$Latitude)

结果如下：

使用 partition.tree() 函数对涉及两个或更多变量的树进行分区。treeModel 作为树对象传递。ordvars=c("经度","纬度") 表示用于绘图的变量顺序。经度代表 x 轴，而 纬度 代表 y 轴。add=TRUE 表示添加到现有图形：

    > partition.tree(treeModel, ordvars=c("Longitude","Latitude"), add=TRUE)

结果如下：

步骤 5 - 改进模型

树中的叶子节点数量控制着树的灵活性。叶子节点的数量表示它们将树分割成多少个单元格。每个节点必须包含一定数量的点，并且添加节点必须至少减少一定的错误。min.dev 的默认值是 0.01。

接下来，我们将 min.dev 的值降低到 0.001。

使用 tree() 函数实现分类树。公式表达式基于变量 纬度 和 经度 的总和。总和的结果存储在 MedianHouseValue 的对数值中。data=realEstate 表示在其中的数据框中优先解释公式、权重和子集。min.dev 的值表示必须至少是根节点偏差的 0.001 倍才能进行节点分割。

结果数据框存储在 treeModel2 数据框中：

    > treeModel2 <- tree(log(MedianHouseValue) ~ Longitude + Latitude, data=realEstate, mindev=0.001)

将显示 treeModel2 的摘要。摘要显示使用的公式，以及树中的终端节点或叶子节点的数量。还显示了残差的统计分布。

使用 summary() 函数显示 treeModel2 的统计摘要。希望进行摘要的数据框是 treeModel2，它作为输入参数传递。

偏差在这里意味着均方误差：

    > summary(treeModel2)

结果如下：

与 treeModel 的摘要相比，treeModel2 中的叶子节点值从 12 增加到 68。对于 treeModel 和 treeModel2，偏差值分别从 0.1666 变为 0.1052。

将显示 treeModel2 数据框的图形版本。plot() 函数是用于绘图 R 对象的通用函数。将 treeModel2 数据框作为函数值传递：

    > plot(treeModel2)

结果如下：

在显示 treeModel2 数据框的图形版本后，需要插入文本以显示每个节点和叶子节点的值。使用 text() 函数在给定的坐标处插入向量标签中给出的字符串：

    > text(treeModel2, cex=.65)

结果如下：

在公式扩展中包含所有变量。

使用tree()函数实现分类树。公式表达式基于所有变量。

结果数据框存储在treeModel3数据框中：

    > treeModel3 <- tree(log(MedianHouseValue) ~ ., data=realEstate)

将显示treeModel3的摘要。摘要显示了使用的公式以及树中的终端节点或叶子节点的数量。还显示了残差的统计分布。

使用summary()函数显示treeModel3的统计摘要。希望进行摘要的数据框是treeModel3，它作为输入参数传递。

偏差在这里表示均方误差：

    > summary(treeModel3)

结果如下：

公式明确指出realEstate数据集中的所有变量。

将显示treeModel3数据框的图形版本。plot()函数是用于绘制 R 对象的通用函数。treeModel3数据框作为函数值传递：

    > plot(treeModel3)

结果如下：

在显示treeModel3数据框的图形版本后，需要插入文本以显示每个节点和叶子节点的值。使用text()函数在给定的坐标处插入向量标签中的字符串：

    > text(treeModel3, cex=.75)

结果如下：

决策树学习 - 预测股票运动方向

股票交易是统计学家试图解决的最具挑战性的问题之一。有多个技术指标，例如趋势方向、动量或市场中的动量不足、盈利潜力的波动性和用于监测市场流行度的成交量等。这些指标可以用来创建策略以创造高概率的交易机会。可以花费数天/周/月来发现技术指标之间的关系。可以使用像决策树这样的高效且节省时间的工具。决策树的主要优势是它是一个强大且易于解释的算法，为良好的起点提供了帮助。

准备工作

为了执行决策树分类，我们将使用从股票市场数据集中收集的数据集。

第 1 步 - 收集和描述数据

要使用的数据集是美国银行 2012 年 1 月 1 日至 2014 年 1 月 1 日的每日收盘价。此数据集在yahoo.com/上免费提供，我们将从那里下载数据。

如何做到这一点...

让我们深入了解。

第 2 步 - 探索数据

需要在第一步加载以下包：

> install.packages("quantmod")
> install.packages("rpart")
> install.packages("rpart.plot")

注意

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）上进行了测试。

上述每个库都需要安装：

> library("quantmod")
> library("rpart")
> library("rpart.plot")

让我们下载数据。我们将首先标记所需数据的时间段的开始和结束日期。

使用as.Date()函数将字符表示和Date类的对象转换为日历日期。

数据集的起始日期存储在startDate中，它代表日历日期的字符向量表示。表示的格式为YYYY-MM-DD：

    > startDate = as.Date("2012-01-01")

数据集的结束日期存储在endDate中，它代表日历日期的字符向量表示。表示的格式为YYYY-MM-DD：

    > endDate = as.Date("2014-01-01")

使用getSymbols()函数加载数据。该函数从多个来源加载数据，无论是本地还是远程。数据被检索并存储在指定的env中。env的默认值是.GlobalEnv。BAC是字符向量，指定要加载的符号名称。src = yahoo指定数据来源方法：

    > getSymbols("BAC", env = .GlobalEnv,  src = "yahoo", from = startDate, to = endDate)

步骤 3 - 计算指标

相对强弱指数（Relative Strength Index）已计算。它是最近上升价格变动与绝对价格变动的比率。使用RSI()函数来计算相对强弱指数。BAC符号用作价格序列。n = 3代表移动平均的周期数。结果存储在relativeStrengthIndex3数据框中：

> relativeStrengthIndex3 <- RSI(Op(BAC), n= 3) 

显示relativeStrengthIndex3的值：

    > relativeStrengthIndex3

结果如下：

计算移动平均。指数移动平均用于技术分析和作为技术指标。在简单移动平均中，序列中的每个值具有相等的权重。时间序列之外的价值不包括在平均中。然而，指数移动平均是一个累积计算，包括所有数据。过去的数据具有递减的价值，而最近的数据值具有更大的贡献。

EMA()使用BAC符号，并用作价格序列。n = 5代表平均的时间周期。结果存储在exponentialMovingAverage5数据框中：

    > exponentialMovingAverage5 <- EMA(Op(BAC),n=5)

显示exponentialMovingAverage5的值：

    > exponentialMovingAverage5

结果如下：

探索exponentialMovingAverage5数据框的维度。dim()函数返回exponentialMovingAverage5框架的维度。将exponentialMovingAverage5数据框作为输入参数传递。结果清楚地表明有 502 行数据和 1 列：

    > dim(exponentialMovingAverage5)

结果如下：

探索exponentialMovingAverage5数据框的内部结构。str()函数显示数据框的内部结构。将exponentialMovingAverage5作为 R 对象传递给str()函数：

    > str(exponentialMovingAverage5)

结果如下：

计算价格和计算出的exponentialMovingAverage5（例如，五年指数移动平均值）之间的差异。结果存储在exponentialMovingAverageDiff数据框中：

    > exponentialMovingAverageDiff <- Op(BAC)-exponentialMovingAverage5

比较 BAC 系列快速移动平均与 BAC 系列慢速移动平均。BAC作为价格矩阵传递。fast = 12表示快速移动平均的周期数，slow = 26表示慢速移动平均的周期数，signal = 9表示移动平均的信号：

    > MACD <- MACD(Op(BAC),fast = 12, slow = 26, signal = 9)

显示 MACD 值：

    > MACD

结果如下：

打印 MACD 数据框。head()函数返回MACD数据框的第一部分。MACD数据框作为输入参数传递：

    > head(MACD)

结果如下：

捕获信号线作为指标。结果存储在MACDsignal数据框中：

    > MACDsignal <- MACD[,2]

显示MACDsignal值：

    > MACDsignal

结果如下：

确定收盘价相对于高低范围的中间位置。为了确定每天收盘价相对于高低范围的位置，使用随机振荡器。SMI()函数用于动量指标。

BAC是包含高低收盘价矩阵。n = 13表示周期数。slow=25表示双平滑的周期数。fast=2表示初始平滑的周期数。signal=9表示信号线的周期数。结果存储在stochasticOscillator数据框中：

    > stochasticOscillator <- SMI(Op(BAC),n=13,slow=25,fast=2,signal=9)

显示stochasticOscillator值：

    > stochasticOscillator

结果如下：

捕获振荡器作为指标。结果存储在stochasticOscillatorSignal数据框中：

    > stochasticOscillatorSignal <- stochasticOscillator[,1]

显示stochasticOscillatorSignal值：

    > stochasticOscillatorSignal

结果如下：

第 4 步 - 准备变量以构建数据集

计算收盘价和开盘价之间的差异。Cl代表收盘价，Op代表开盘价。结果存储在PriceChange数据框中：

    > PriceChange <- Cl(BAC) - Op(BAC)

显示PriceChange值：

    > PriceChange

结果如下：

创建一个二元分类变量。ifelse()函数使用一个测试表达式来返回值，该值本身是一个向量，其长度与测试表达式相同。如果test表达式的对应值为TRUE，则返回x中的元素；如果test表达式的对应值为FALSE，则返回y中的元素。

在这里，PriceChange>0 是测试函数，将在逻辑模式下进行测试。UP 和 DOWN 执行逻辑测试。结果随后存储在 binaryClassification 数据框中：

    > binaryClassification <- ifelse(PriceChange>0,"UP","DOWN")

显示 binaryClassification 值：

    > binaryClassification

结果如下：

探索 binaryClassification 数据框的内部结构。str() 函数显示数据框的内部结构。binaryClassification 作为 R 对象传递给 str() 函数：

    > str(binaryClassification)

结果如下：

创建要使用的数据集。data.frame() 函数用于根据紧密耦合的变量集创建数据框。这些变量具有矩阵的性质。传递给 data.frame() 的参数变量有 relativeStrengthIndex3、exponentialMovingAverageDiff、MACDsignal、stochasticOscillator 和 binaryClassification。

结果随后存储在 DataSet 数据框中：

> AAPLDataSetNew >-
data.frame(weekDays,exponentialMovingAverageDiffRound,
binaryClassification) 

显示 DataSet 值：

    > DataSet

结果如下：

打印 DataSet 数据框。head() 函数返回 DataSet 数据框的第一部分。DataSet 数据框作为输入参数传递：

    > head(DataSet)

结果如下：

探索 DataSet 数据框的内部结构。str() 函数显示数据框的内部结构。DataSet 作为 R 对象传递给 str() 函数：

    > str(DataSet)

结果如下：

命名列。c() 函数用于将参数组合成向量。

传递给 c() 函数的参数变量有 relativeStrengthIndex3、exponentialMovingAverageDiff、MACDsignal、stochasticOscillator 和 binaryClassification：

    > colnames(DataSet) <- c("relativeStrengthIndex3", "exponentialMovingAverageDiff", "MACDsignal", "stochasticOscillator", "binaryClassification")

显示 colnames(DataSet) 值：

    > colnames(DataSet)

结果如下：

删除要计算指标的数据：

    > DataSet <- DataSet[-c(1:33),]

显示 DataSet 值：

    > DataSet

结果如下：

打印 DataSet 数据框。head() 函数返回 DataSet 数据框的第一部分。DataSet 数据框作为输入参数传递：

    > head(DataSet)

结果如下：

探索 DataSet 数据框的内部结构。str() 函数显示数据框的内部结构。DataSet 作为 R 对象传递给 str() 函数：

    > str(DataSet)

结果如下：

探索DataSet数据框的维度。dim()函数返回DataSet框的维度。将DataSet数据框作为输入参数传递。结果显示，共有 469 行数据和 5 列：

    > dim(DataSet)

结果如下：

构建训练数据集。DataSet数据框中的三分之二元素将用作训练数据集，而DataSet数据框中的一分之一元素将用作测试数据集。

训练数据集将存储在TrainingDataSet中：

    > TrainingDataSet <- DataSet[1:312,]

显示TrainingDataSet的值：

    > TrainingDataSet

结果如下：

探索TrainingDataSet数据框的内部结构。str()函数显示数据框的内部结构。将TrainingDataSet作为 R 对象传递给str()函数：

    > str(TrainingDataSet)

结果如下：

训练数据集将存储在TestDataSet中：

    > TestDataSet <- DataSet[313:469,]

显示TestDataSet的值：

    > TestDataSet

结果如下：

探索TestDataSet数据框的内部结构。str()函数显示数据框的内部结构。将TestDataSet作为 R 对象传递给str()函数：

    > str(TestDataSet)

结果如下：

步骤 5 - 构建模型

通过指定指标构建树模型。将使用rpart()函数。它将拟合模型。binaryClassification是结果，使用relativeStrengthIndex3、exponentialMovingAverageDiff、MACDsignal和stochasticOscillator的总和作为预测因子。data=TrainingDataSet表示数据框。cp=.001表示复杂性参数。该参数的主要作用是通过剪枝来节省计算时间。结果随后存储在DecisionTree数据框中：

    > DecisionTree <- rpart(binaryClassification~relativeStrengthIndex3+exponentialMovingAverageDiff+MACDsignal+stochasticOscillator,data=TrainingDataSet, cp=.001)

绘制树模型。将使用prp()函数绘制DecisionTree数据框。type=2将交替节点垂直移动：

    > prp(DecisionTree,type=2)

结果如下：

显示DecisionTree数据框的cp表。使用printcp()函数。将DecisionTree作为输入传递：

    > printcp(DecisionTree)

结果如下：

绘制树的几何平均。使用plotcp()函数。它提供了DecisionTree数据框交叉验证结果的视觉表示：

    > plotcp(DecisionTree,upper="splits")

结果如下：

步骤 6 - 改进模型

在剪枝后改进模型。使用prune()函数。DecisionTree是作为输入传递的数据框。cp=0.041428已被采用，因为这是最低的交叉验证错误值（x 错误）：

    > PrunedDecisionTree <- prune(DecisionTree,cp=0.041428)

绘制tree模型。将使用prp()函数绘制DecisionTree数据框。type=4将交替节点垂直移动：

    > prp(PrunedDecisionTree, type=4)

结果如下：

测试模型：

> table(predict(PrunedDecisionTree,TestDataSet), TestDataSet[,5],dnn=list('predicted','actual')) 

结果如下：

简单贝叶斯 - 预测股票运动的方向

股票交易是统计学家试图解决的最具挑战性的问题之一。市场中有多个技术指标，例如趋势方向、动量或市场动量的缺乏、波动性以衡量盈利潜力，以及用于监控市场流行度的成交量等，仅举几例。这些指标可以用来创建策略以捕捉高概率的交易机会。可能需要花费数日/数周/数月来发现技术指标之间的关系。可以使用像决策树这样的高效且节省时间的工具。决策树的主要优势在于它是一个强大且易于解释的算法，这为良好的起点提供了帮助。

准备工作

为了执行简单贝叶斯，我们将使用从股票市场数据集中收集的数据集。

第 1 步 - 收集和描述数据

要使用的数据集是 2012 年 1 月 1 日至 2014 年 1 月 1 日苹果公司每日收盘价。此数据集在www.yahoo.com/上免费提供，我们将从那里下载数据。

如何做到这一点...

让我们深入了解细节。

第 2 步 - 探索数据

以下包需要在执行第一步时加载：

    > install.packages("quantmod")
    > install.packages("lubridate")
    > install.packages("e1071")

注意

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）上进行了测试

以下每个库都需要安装：

    > library("quantmod")
    > library("lubridate")
    > library("e1071")

让我们下载数据。我们首先标记所需数据的时间段的开始和结束日期。

使用as.Date()函数将字符表示和Date类的对象转换为日历日期。

数据集的开始日期存储在startDate中，它表示日历日期的字符向量表示。表示的格式是YYYY-MM-DD：

    > startDate = as.Date("2012-01-01")

数据集的结束日期存储在endDate中，它表示日历日期的字符向量表示。表示的格式是 YYYY-MM-DD：

    > endDate = as.Date("2014-01-01")

使用getSymbols()函数加载数据。该函数从多个来源加载数据，无论是本地还是远程。数据被检索并保存在指定的env中。对于env，默认值是.GlobalEnv。AAPL是字符向量，指定要加载的符号名称。src = yahoo指定了数据来源方法：

    > getSymbols("AAPL", env = .GlobalEnv, src = "yahoo", from = startDate,  to = endDate)

探索数据可用的星期几。使用 wday() 函数。该函数以十进制格式返回星期几。AAPL 代表数据框。label = TRUE 将星期几显示为字符串，例如，星期日。结果随后存储在 weekDays 数据框中：

    > weekDays <- wday(AAPL, label=TRUE)

打印 weekDays 数据框。head() 函数返回 weekDays 数据框的前部分。将 weekDays 数据框作为输入参数传递：

    > head(weekDays)

结果如下：

第 3 步 - 准备构建数据集的变量

计算收盘价和开盘价之间的差异。Cl 代表收盘价，Op 代表开盘价。结果存储在 changeInPrices 数据框中：

    > changeInPrices <- Cl(AAPL) - Op(AAPL)

打印 changeInPrices 数据框。head() 函数返回 changeInPrices 数据框的前部分。将 changeInPrices 数据框作为输入参数传递：

    > head(changeInPrices)

结果如下：

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计，以生成 changeInPrices 数据框的结果摘要：

    > summary(changeInPrices)

结果如下：

探索 changeInPrices 数据框的维度。dim() 函数返回 changeInPrices 框的维度。将 changeInPrices 数据框作为输入参数传递。结果清楚地表明有 502 行数据和 1 列：

    > dim(changeInPrices)

结果如下：

创建一个二元分类变量。ifelse() 函数使用测试表达式返回值，该值本身是一个向量，其长度与测试表达式相同。如果测试表达式的对应值为 TRUE，则从 x 中返回向量中的一个元素，如果测试表达式的对应值为 FALSE，则从 y 中返回。

在这里，changeInPrices>0 是一个测试函数，用于测试逻辑模式。UP 和 DOWN 执行逻辑测试。结果随后存储在 binaryClassification 数据框中：

    > binaryClassification <- ifelse(changeInPrices>0,"UP","DOWN")

显示 binaryClassification 值：

    > binaryClassification

结果如下：

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计，以生成 binaryClassification 数据框的结果摘要：

    > summary(binaryClassification)

结果如下：

创建要使用的数据集。使用 data.frame() 函数根据紧密耦合的变量集创建数据框。这些变量具有矩阵的性质。

将作为data.frame()参数传递的变量是weekDays和binaryClassification。结果随后存储在DataSet数据框中：

    > AAPLDataSet <- data.frame(weekDays,binaryClassification)

显示AAPLDataSet值：

    > AAPLDataSet

结果如下：

打印AAPLDataSet数据框。head()函数返回AAPLDataSet数据框的前部分。将AAPLDataSet数据框作为输入参数传递：

    > head(AAPLDataSet)

结果如下：

探索AAPLDataSet数据框的维度。dim()函数返回AAPLDataSet数据框的维度。将AAPLDataSet数据框作为输入参数传递。结果明确指出有 502 行数据和 2 列：

    > dim(AAPLDataSet)

结果如下：

第 4 步 - 构建模型

通过指定指标构建朴素贝叶斯分类器。将使用naiveBayes()函数。该函数使用贝叶斯规则来计算给定一组独立预测变量的后验概率。该函数假设度量预测变量服从高斯分布。"NaiveBayesclassifier"是函数的输出结果，其中独立变量是AAPLDataSet[,1]，因变量是AAPLDataSet[,2]：

    > NaiveBayesclassifier <- naiveBayes(AAPLDataSet[,1], AAPLDataSet[,2])

显示NaiveBayesclassifier结果：

    > NaiveBayesclassifier

结果如下：

结果覆盖整个数据集，并显示价格增加或减少的概率。其本质上是看跌的。

第 5 步 - 创建新的、改进模型的数据

制定一个复杂的策略，展望超过一天。对模型计算 5 年的移动平均。EMA()使用 AAPL 符号作为价格序列。"n = 5"代表平均的时间段。结果随后存储在exponentialMovingAverage5数据框中：

    > exponentialMovingAverage5 <- EMA(Op(AAPL),n = 5)

显示exponentialMovingAverage5值：

    > exponentialMovingAverage5

结果如下：

探索价格变化的摘要。使用summary()函数。该函数提供一系列描述性统计量，以生成exponentialMovingAverage5数据框的结果摘要：

    > summary(exponentialMovingAverage5)

结果如下：

对模型计算 10 年的移动平均。

EMA()使用 AAPL 符号作为价格序列。"n = 10"代表平均的时间段。结果随后存储在exponentialMovingAverage10数据框中：

    > exponentialMovingAverage10 <- EMA(Op(AAPL),n = 10)

显示exponentialMovingAverage10值：

    > exponentialMovingAverage10

结果如下：

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 exponentialMovingAverage10 数据框的结果摘要：

    > summary(exponentialMovingAverage10)

结果如下：

探索 exponentialMovingAverage10 数据框的维度。dim() 函数返回 exponentialMovingAverage10 框的维度。将 exponentialMovingAverage10 数据框作为输入参数传递。结果清楚地表明有 502 行数据和 1 列：

    > dim(exponentialMovingAverage10)

结果如下：

计算 exponentialMovingAverage5 和 exponentialMovingAverage10 之间的差异：

    > exponentialMovingAverageDiff <- exponentialMovingAverage5 - exponentialMovingAverage10

显示 exponentialMovingAverageDiff 值：

    > exponentialMovingAverageDiff

结果如下：

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 exponentialMovingAverageDiff 数据框的结果摘要：

    > summary(exponentialMovingAverageDiff)

结果如下：

将 exponentialMovingAverageDiff 数据框四舍五入到两位有效数字：

    > exponentialMovingAverageDiffRound <- round(exponentialMovingAverageDiff, 2)

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 exponentialMovingAverageDiffRound 数据框的结果摘要：

    > summary(exponentialMovingAverageDiffRound)

结果如下：

步骤 6 - 改进模型

创建用于的数据集。使用 data.frame() 函数根据一组紧密耦合的变量创建数据框。这些变量具有矩阵的性质。传递给 data.frame() 的参数变量是 weekDays、exponentialMovingAverageDiffRound 和 binaryClassification。结果存储在 AAPLDataSetNew 数据框中：

> AAPLDataSetNew <- data.frame(weekDays,exponentialMovingAverageDiffRound, binaryClassification) 

显示 AAPLDataSetNew 值：

> AAPLDataSetNew 

结果如下：

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 AAPLDataSetNew 数据框的结果摘要：

    > summary(AAPLDataSetNew)

结果如下：

    > AAPLDataSetNew <- AAPLDataSetNew[-c(1:10),]

结果如下：

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 AAPLDataSetNew 数据框的结果摘要：

> summary(AAPLDataSetNew) 

结果如下：

探索 AAPLDataSetNew 数据框的维度。dim() 函数返回 AAPLDataSetNew 框的维度。将 AAPLDataSetNew 数据框作为输入参数传递。结果明确指出有 492 行数据和 3 列：

    > dim(AAPLDataSetNew)

结果如下：

构建训练数据集。AAPLDataSetNew 数据框中的三分之二元素将用作训练数据集，而 AAPLDataSetNew 数据框中的一分之一元素将用作测试数据集。

训练数据集将存储在 trainingDataSet 数据框中：

> trainingDataSet <- AAPLDataSetNew[1:328,] 

探索 trainingDataSet 数据框的维度。dim() 函数返回 trainingDataSet 数据框的维度。将 trainingDataSet 数据框作为输入参数传递。结果明确指出有 328 行数据和 3 列：

    > dim(trainingDataSet)

结果如下：

探索价格变化的摘要。使用 trainingDataSet() 函数。该函数提供一系列描述性统计量，以生成 trainingDataSet 数据框的结果摘要：

    > summary(trainingDataSet)

结果如下：

训练数据集将存储在 TestDataSet 数据框中：

    > TestDataSet <- AAPLDataSetNew[329:492,]

探索 TestDataSet 数据框的维度。dim() 函数返回 TestDataSet 框的维度。将 TestDataSet 数据框作为输入参数传递。结果明确指出有 164 行数据和 3 列：

    > dim(TestDataSet)

结果如下：

    > summary(TestDataSet)

结果如下：

通过指定指标构建朴素贝叶斯分类器。将使用 naiveBayes() 函数。它使用贝叶斯规则计算给定一组类别变量和独立预测变量后的后验概率。该函数假设度量预测变量的高斯分布。

exponentialMovingAverageDiffRoundModel 是函数的输出结果，其中自变量是 trainingDataSet[,1:2]，因变量是 trainingDataSet[,3]：

> exponentialMovingAverageDiffRoundModel <-
naiveBayes(trainingDataSet[,1:2],trainingDataSet[,3])

显示 exponentialMovingAverageDiffRoundModel 结果：

    > exponentialMovingAverageDiffRoundModel

结果如下：

测试结果：

    > table(predict(exponentialMovingAverageDiffRoundModel,TestDataSet),
TestDataSet[,3],dnn=list('Predicted','Actual')) 

结果如下：

随机森林 - 货币交易策略

在进行技术分析后，可以科学地实现预测外汇市场未来价格趋势的目标。外汇交易者根据市场趋势、成交量、范围、支撑和阻力水平、图表模式和指标等多种技术分析制定策略，并使用不同时间框架的图表进行多时间框架分析。基于过去市场行动的统计数据，如过去价格和过去成交量，创建技术分析策略以评估资产。分析的主要目标不是衡量资产的基本价值，而是计算市场的历史表现所指示的未来市场表现。

准备工作

为了执行随机森林，我们将使用从美元和英镑数据集收集的数据集。

第一步 - 收集和描述数据

将使用标题为 PoundDollar.csv 的数据集。数据集是标准格式。有 5,257 行数据和 6 个变量。数值变量如下：

• 日期

• 开盘价

• 最高价

• 最低价

• 收盘价

• 成交量

如何操作...

让我们深入了解细节。

第二步 - 探索数据

作为第一步要执行，以下包需要加载：

> install.packages("quantmod")
> install.packages("randomForest")
> install.packages("Hmisc")

备注

版本信息：本页代码在 R 版本 3.3.0（2016-05-03）中进行了测试。

以下每个库都需要安装：

> library("quantmod")
> library("randomForest")
> library("Hmisc")

让我们探索数据并了解变量之间的关系。我们将首先导入名为 PoundDollar.csv 的 CSV 数据文件。我们将把数据保存到 PoundDollar 数据框中：

    > PoundDollar <- read.csv("d:/PoundDollar.csv")

打印 PoundDollar 数据框。head() 函数返回 PoundDollar 数据框的前一部分。PoundDollar 数据框作为输入参数传递：

    > head(PoundDollar)

结果如下：

打印 PoundDollar 数据框的摘要。summary() 函数是一个多功能函数。summary() 是一个通用函数，它提供了与单个对象或数据框相关的数据的摘要。PoundDollar 数据框作为 R 对象传递给 summary() 函数：

    > summary(PoundDollar)

结果如下：

探索 PoundDollar 数据框的维度。dim() 函数返回 PoundDollar 框的维度。PoundDollar 数据框作为输入参数传递。结果清楚地表明有 5,257 行数据和 7 列：

    > dim(PoundDollar)

结果如下：

第三步 - 准备变量以构建数据集

表示日历日期和时间。as.POSIXlt() 函数将对象操作为表示日期和时间。PoundDollar 作为参数传递。format="%m/%d/%y %H:%M 表示日期时间格式。结果存储在 DateAndTime 数据框中：

    > DateAndTime <- as.POSIXlt(PoundDollar[,2],format="%m/%d/%y %H:%M")

捕获 最高价、最低价 和 收盘价 值：

    > HighLowClose <- PoundDollar[,4:6]

PoundDollar数据框捕获了第四、第五和第六列中的High、Low和Close值。打印HighLowClose数据框。head()函数返回HighLowClose数据框的第一部分。HighLowClose数据框被作为输入参数传递：

    > head(HighLowClose)

结果如下：

打印HighLowClose数据框的摘要。summary()函数是一个多功能函数。summary()是一个泛型函数，它提供了与单个对象或数据框相关的数据的摘要。HighLowClose数据框被作为 R 对象传递给summary()函数：

    > summary(HighLowClose)

结果如下：

探索HighLowClose数据框的内部结构。str()函数显示数据框的内部结构。将HighLowClose作为 R 对象传递给str()函数：

    > str(HighLowClose)

结果如下：

创建要使用的数据集。使用data.frame()函数根据紧密耦合的变量集创建数据框。这些变量具有矩阵的性质。将HighLowClose作为参数传递给data.frame()。然后将结果存储在HighLowClosets数据框中。row.names=DateAndTime表示一个整数字符串，指定用作行名的列。结果存储在HighLowClose数据框中：

> HighLowClosets <- data.frame(HighLowClose, row.names=DateAndTime) 

描述数据集。describe()函数提供项目分析。HighLowClosets作为输入参数传递：

    > describe(HighLowClosets)

结果如下：

创建时间序列对象。使用as.xts()函数。它将任意类别的数据对象转换为xts类，而不丢失原始格式的任何属性。HighLowClosets被作为输入对象传递：

    > HighLowClosexts <- as.xts(HighLowClosets)

计算布林带。布林带是一种范围指标，它从移动平均数计算标准差。布林带遵循的逻辑是，货币对的价格最有可能趋向于其平均值，因此当它偏离太多，例如两个标准差之外时，它将回溯到其移动平均数。使用BBands()函数来计算布林带。HighLowClosexts被作为对象传递，该对象被转换为包含高低收盘价的矩阵。n=20表示移动平均数的周期数。SMA 命名要调用的函数。sd=2表示两个标准差：

    > BollingerBands <- BBands(HighLowClosexts,n=20,SMA,sd=2)

描述数据集。describe()函数提供项目分析。BollingerBands作为输入参数传递：

    > describe(BollingerBands)

结果如下：

构建上限带：

    > Upper <- BollingerBands$up - HighLowClosexts$Close

打印上界数据框的摘要。summary() 函数是一个多功能函数。summary() 是一个通用函数，它提供了与单个对象或数据框相关的数据的摘要。Upper 数据框作为 R 对象传递给 summary() 函数：

    > summary(Upper)

结果如下：

构建下界带：

    > Lower <- BollingerBands$dn - HighLowClosexts$Close

打印下界数据框的摘要。summary() 函数是一个多功能函数。summary() 是一个通用函数，它提供了与单个对象或数据框相关的数据的摘要。下界数据框作为 R 对象传递给 summary() 函数：

    > summary(Upper)

结果如下：

构建中间带：

    > Middle <- BollingerBands$mavg - HighLowClosexts$Close

打印中间数据框的摘要。summary() 函数是一个多功能函数。summary() 是一个通用函数，它提供了与单个对象或数据框相关的数据的摘要。Middle 数据框作为 R 对象传递给 summary() 函数：

    > summary(Middle)

结果如下：

计算百分比变化。使用 Delt() 函数计算给定序列从一个时期到另一个时期的百分比变化。k=1 表示在各个时期的变化。结果存储在 PercentageChngpctB 数据框中：

    > PercentageChngpctB <- Delt(BollingerBands$pctB,k=1)

描述数据集。describe() 函数提供项目分析。PercentageChngpctB 作为输入参数传递：

    > describe(PercentageChngpctB)

结果如下：

计算上界数据框的百分比变化。k=1 表示在各个时期的变化：

    > PercentageChngUp <- Delt(Upper,k=1)

描述数据集。describe() 函数提供项目分析。PercentageChngUp 作为输入参数传递：

    > describe(PercentageChngUp)

结果如下：

计算下界数据框的百分比变化。k=1 表示在各个时期的变化：

    > PercentageChngLow <- Delt(Lower, k=1)

描述数据集。describe() 函数提供项目分析。PercentageChngLow 作为输入参数传递：

    > describe(PercentageChngLow)

结果如下：

计算中间数据框的百分比变化。k=1 表示在各个时期的变化：

    > PercentageChngMid <- Delt(Middle,k=1)

描述数据集。describe() 函数提供项目分析。PercentageChngMid 作为输入参数传递：

    > describe(PercentageChngMid)

结果如下：

计算变量 HighLowClosexts$Close 的百分比变化。k=1 表示在各个时期的变化：

    > Returns <- Delt(HighLowClosexts$Close, k=1)

第 4 步 - 构建模型

创建二元分类变量。ifelse() 函数使用测试表达式返回值，该值本身是一个向量，其长度与测试表达式相同。如果测试表达式的对应值为 TRUE，则从 x 中返回一个元素；如果测试表达式的对应值为 FALSE，则从 y 中返回一个元素。

在这里，Returns>0 是测试函数，需要在逻辑模式下进行测试。UP 和 DOWN 执行逻辑测试。结果随后存储在 binaryClassification 数据框中：

> binaryClassification <- ifelse(Returns>0,"Up","Down") 

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 binaryClassification 数据框的结果摘要：

    > summary(binaryClassification)

结果如下：

将类别回退一个：

    > ClassShifted <- binaryClassification[-1]

结合所有特征。使用 data.frame() 函数根据紧密耦合的变量集创建数据框。这些变量具有矩阵的性质。

传递给 data.frame() 的参数变量有 Upper、Lower、Middle、BollingerBands$pctB、PercentageChngpctB、PercentageChngUp、PercentageChngLow 和 PercentageChngMid。结果随后存储在 FeaturesCombined 数据框中：

    > FeaturesCombined <- data.frame(Upper, Lower, Middle, BollingerBands$pctB, PercentageChngpctB, PercentageChngUp, PercentageChngLow, PercentageChngMid)

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 FeaturesCombined 数据框的结果摘要：

    > summary(FeaturesCombined)

结果如下：

匹配类别：

    > FeaturesShifted <- FeaturesCombined[-5257,]

结合 FeaturesShifted 和 ClassShifted 数据框。传递给 data.frame() 的参数变量是 FeaturesShifted 和 ClassShifted。结果随后存储在 FeaturesClassData 数据框中：

    > FeaturesClassData <- data.frame(FeaturesShifted, ClassShifted)

探索价格变化的摘要。使用 summary() 函数。该函数提供一系列描述性统计量，以生成 FeaturesClassData 数据框的结果摘要：

    > summary(FeaturesClassData)

结果如下：

计算指标正在被移除：

    > FinalModelData <- FeaturesClassData[-c(1:20),]

命名列。使用 c() 函数将参数组合成向量：

    > colnames(FinalModelData) <- c("pctB","LowDiff","UpDiff","MidDiff","PercentageChngpctB","PercentageChngUp","PercentageChngLow","PercentageChngMid","binaryClassification")

探索 FinalModelData 数据框的内部结构。str() 函数显示数据框的内部结构。FinalModelData 作为 R 对象传递给 str() 函数：

    > str(FinalModelData)

结果如下：

设置初始随机变量：

    > set.seed(1)

使用类别（第 9 列）评估特征（第 1 至 9 列）以找到每棵树的最佳特征数量。"FinalModelData[,-9]" 表示预测变量数据框，"FinalModelData[,9]" 表示响应变量数据框。"ntreeTry=100" 表示在调整步骤中使用的树的数量。"stepFactor=1.5" 表示每次迭代的值，"mtry" 通过这个值膨胀（或缩水），"improve=0.01" 表示搜索必须继续的（相对）出袋误差的改善量。"trace=TRUE" 表示是否打印搜索的进度。"dobest=FALSE" 表示是否使用找到的最佳 "mtry" 运行森林：

    > FeatureNumber <- tuneRF(FinalModelData[,-9], FinalModelData[,9], ntreeTry=100, stepFactor=1.5, improve=0.01, trace=TRUE, plot=TRUE, dobest=FALSE)

使用所有特征进行分类预测，每棵树有两个特征。使用 "randomForest()" 函数。data=FinalModelData 表示包含模型中变量的数据框。"mtry=2" 表示在每次分割中随机采样的变量作为候选者的数量。"ntree=2000" 表示要生长的树的数量。"keep.forest=TRUE" 表示森林将保留在输出对象中。"importance=TRUE" 表示要评估预测变量的重要性：

    > RandomForest <- randomForest(binaryClassification~., data=FinalModelData, mtry=2,  ntree=2000, keep.forest=TRUE, importance=TRUE)

结果如下：

绘制随机森林：

    > varImpPlot(RandomForest, main = 'Random Forest: Measurement of Importance of Each Feature',pch=16,col='blue' )

结果如下：

支持向量机 - 货币交易策略

外汇市场是一个国际交易市场，各国货币可以自由买卖。一种货币的价格仅由市场参与者决定，由供求关系驱动。交易通过个别合约进行。标准合约规模（也称为一手）通常是 100,000 单位。这意味着，对于每份标准合约，控制的是 100,000 单位的基础货币。对于这个合约规模，每个点（最小的价格变动单位）价值 10 美元。根据交易者的交易策略，头寸可以维持非常短的时间，也可以维持更长的时间，甚至数年。有几个工具允许交易者理解和在市场上做出决策，这些工具基本上分为基本面分析或技术分析。基本面分析考虑了政治和经济信息的持续交换。技术分析基本上基于价格、时间和成交量——货币达到的最低和最高价格、时间段、交易次数等。技术分析还假设市场的重复性，它很可能在未来再次执行，就像它在过去已经执行的那样。它分析过去的报价，并根据统计和数学计算预测未来的价格。

准备中

为了执行支持向量机，我们将使用从美元和英镑数据集中收集的数据集。

步骤 1 - 收集和描述数据

将使用标题为 PoundDollar.csv 的数据集。数据集是标准格式。有 5,257 行数据，6 个变量。数值变量如下：

• 日期

• 开盘价

• 最高价

• 低

• 收盘价

• 成交量

如何操作...

让我们深入了解细节。

步骤 2 - 探索数据

作为第一步需要加载以下包：

> install.packages("quantmod")
> install.packages("e1071")
> install.packages("Hmisc")
> install.packages("ggplot2")

注意

版本信息：本页代码在 R 版本 3.3.0（2016-05-03）中进行了测试。

以下每个库都需要安装：

> library("quantmod")
> library("e1071")
> library("Hmisc")
> install.packages("ggplot2")

让我们探索数据并了解变量之间的关系。我们将从导入名为 PoundDollar.csv 的 CSV 数据文件开始。我们将数据保存到 PoundDollar 数据框中：

    > PoundDollar <- read.csv("d:/PoundDollar.csv")

打印 PoundDollar 数据框。head() 函数返回 PoundDollar 数据框的前一部分。PoundDollar 数据框作为输入参数传递：

    > head(PoundDollar)

结果如下：

探索 PoundDollar 数据框的内部结构。str() 函数显示数据框的内部结构。PoundDollar 作为 R 对象传递给 str() 函数：

    > str(PoundDollar)

结果如下：

步骤 3 - 计算指标

计算相对强弱指数（RSI）。它是最近向上价格变动与绝对价格变动的比率。使用 RSI() 函数计算相对强弱指数。PoundDollar 数据框用作价格序列。n = 3 表示移动平均的周期数。结果存储在 relativeStrengthIndex3 数据框中：

    > relativeStrengthIndex3 <- RSI(Op(PoundDollar), n= 3)

探索价格变化的摘要。使用 summary() 函数。该函数提供了一系列描述性统计量，以生成 relativeStrengthIndex3 数据框的结果摘要：

    > summary(relativeStrengthIndex3)

结果如下：

计算 PoundDollar 序列的 移动平均（MA）。SMA 计算过去一系列观察值的算术平均值。n=50 表示平均的周期数：

    > SeriesMeanAvg50 <- SMA(Op(PoundDollar), n=50)

打印 SeriesMeanAvg50 数据框的摘要。summary() 函数是一个多功能函数。summary() 是一个通用函数，它提供了与单个对象或数据框相关的数据的摘要。SeriesMeanAvg50 数据框作为 R 对象传递给 summary() 函数：

    > summary(SeriesMeanAvg50)

结果如下：

描述数据集。describe() 函数提供项目分析。SeriesMeanAvg50 作为输入参数传递：

    > describe(SeriesMeanAvg50)

结果如下：

测量趋势。找出开盘价与 50 期简单移动平均价之间的差异：

    > Trend <- Op(PoundDollar) - SeriesMeanAvg50

打印 SeriesMeanAvg50 数据框的摘要。Trend 数据框作为 R 对象传递给 summary() 函数：

    > summary(Trend)

结果如下：

计算收盘价和开盘价之间的价格差异。结果存储在数据框 PriceDiff 中：

    > PriceDiff <- Cl(PoundDollar) - Op(PoundDollar)

打印 PriceDiff 数据框的摘要。Trend 数据框作为 R 对象传递给 summary() 函数：

    > summary(PriceDiff)

结果如下：

步骤 4 - 准备变量以构建数据集

创建二元分类变量。ifelse() 函数使用测试表达式返回值，该值本身是一个向量，其长度与测试表达式相同。如果测试表达式的对应值为 TRUE，则返回 x 的元素；如果对应值为 FALSE，则返回 y 的元素。

这里，PriceChange>0 是测试函数，需要在逻辑模式下进行测试。UP 和 DOWN 执行逻辑测试。结果随后存储在 binaryClassification 数据框中：

    > binaryClassification <- ifelse(PriceDiff>0,"UP","DOWN")

打印 binaryClassification 数据框的摘要。Trend 数据框作为 R 对象传递给 summary() 函数：

    > summary(binaryClassification)

结果如下：

结合相对 StrengthIndex3、Trend 和 binaryClassification 数据框。传递给 data.frame() 的参数是 relativeStrengthIndex3、Trend 和 binaryClassification。结果存储在 DataSet 数据框中：

    > DataSet <- data.frame(relativeStrengthIndex3, Trend, binaryClassification)

打印 DataSet 数据框的摘要。Trend 数据框作为 R 对象传递给 summary() 函数：

    > summary(DataSet)

结果如下：

探索 DataSet 数据框的内部结构。str() 函数显示数据框的内部结构。DataSet 作为 R 对象传递给 str() 函数：

    > str(DataSet)

结果如下：

计算指标、创建数据集和删除点：

    > DataSet <- DataSet[-c(1:49),]

探索 DataSet 数据框的维度。dim() 函数返回 DataSet 框的维度。DataSet 数据框作为输入参数传递。结果清楚地表明有 5,208 行数据和 3 列：

> dim(DataSet) 

结果如下：

分离训练数据集：

    > TrainingDataSet <- DataSet[1:4528,]

探索 TrainingDataSet 数据框的维度。dim() 函数返回 TrainingDataSet 框的维度。TrainingDataSet 数据框作为输入参数传递。结果清楚地表明有 4,528 行数据和 3 列：

    > dim(TrainingDataSet)

结果如下：

打印TrainingDataSet数据框的摘要。TrainingDataSet数据框作为 R 对象传递给summary()函数：

    > summary(TrainingDataSet)

结果如下：

分离测试数据集：

    > TestDataSet <- DataSet[4529:6038,]

探索TestDataSet数据框的维度。dim()函数返回TestDataSet框的维度。将TestDataSet数据框作为输入参数传递。结果清楚地表明有 1,510 行数据和 3 列：

    > dim(TestDataSet)

打印TestDataSet数据框的摘要。TestDataSet数据框作为 R 对象传递给summary()函数：

    > summary(TestDataSet)

结果如下：

步骤 5 - 构建模型

使用svm()函数构建支持向量机。使用binaryClassification~relativeStrengthIndex3+Trend作为公式。data=TrainingDataSet用作包含模型变量的数据框。kernel="radial"表示在训练和预测中使用径向基核函数。cost=1表示违反约束的成本。gamma=1/2表示除线性核函数之外所有核函数所需的参数：

    > SVM <- svm(binaryClassification~relativeStrengthIndex3+Trend, data=TrainingDataSet, kernel="radial", cost=1, gamma=1/2)

打印SVM数据框的摘要。SVM数据框作为 R 对象传递给summary()函数：

    > summary(SVM)

结果如下：

为了根据模型对象预测值，我们将使用predict()函数。将SVM作为对象传递。将TrainingDataSet数据对象作为对象传递，在其中查找用于预测的变量：

    > TrainingPredictions <- predict(SVM, TrainingDataSet, type="class")

打印TrainingPredictions数据框的摘要。SVM数据框作为 R 对象传递给summary()函数：

    > summary(TrainingPredictions)

结果如下：

描述数据集。describe()函数提供项目分析。将TrainingPredictions作为输入参数传递：

    > describe(TrainingPredictions)

结果如下：

合并TrainingDataSet和TrainingPredictions数据框。传递给data.frame()函数的参数是TrainingDataSet和TrainingPredictions。结果存储在TrainingDatadata数据框中：

    > TrainingData <- data.frame (TrainingDataSet, TrainingPredictions)

打印TrainingData数据框的摘要。TrainingData数据框作为 R 对象传递给summary()函数：

    > summary(TrainingData)

结果如下：

打印TrainingData：

    > ggplot(TrainingData,aes(x=Trend,y=relativeStrengthIndex3))    +stat_density2d(geom="contour",aes(color=TrainingPredictions))    +labs(,x="Open - SMA50",y="RSI3",color="Training Predictions")

结果如下：

随机梯度下降 - 成人收入

随机梯度下降也称为增量梯度下降，是梯度下降优化方法的一个随机近似，该方法用于最小化一个表示为可微函数之和的目标函数。它通过迭代尝试找到最小值或最大值。在随机梯度下降中，*Q(w)*的真正梯度被一个单例的梯度近似：

当算法遍历训练集时，它会对每个训练示例执行上述更新。可以在训练集上多次遍历，直到算法收敛。如果这样做，则可以在每次遍历中打乱数据以防止循环。典型的实现可能使用自适应学习率，以便算法收敛。

准备工作

为了执行随机梯度下降，我们将使用从人口普查数据收集的数据集来预测收入。

第 1 步 - 收集和描述数据

将使用名为adult.txt的数据集。数据集是标准格式。有 32,561 行数据和 15 个变量。数值变量如下：

• 年龄

• fnlwgt

• 教育年限

• 资本收益

• 资本损失

• 每周工作小时数

非数值变量如下：

• 工作类别

• 教育

• 婚姻状况

• 职业

• 关系

• 种族

• 性别

• 国籍

• 收入范围

如何做到这一点...

让我们深入了解细节。

第 2 步 - 探索数据

以下每个库都需要安装：

> library("klar")
> library("caret")
> library ("stringr")

注意

版本信息：本页面的代码在 R 版本 3.3.0（2016-05-03）中进行了测试。

让我们探索数据并了解变量之间的关系。我们将从导入名为adult.txt的 TXT 数据文件开始。我们将数据保存到labels数据框中：

    > labels <- read.csv("d:/adult.txt")

探索allData数据框的内部结构。str()函数显示数据框的内部结构。将allData作为 R 对象传递给str()函数：

    > str(allData)

结果如下：

第 3 步 - 准备数据

从主文件中获取标签。使用as.factor()函数将allData[,15]向量编码为因子，以确保格式兼容性。然后，结果存储在labels数据框中：

    > labels <- as.factor(allData[,15])

在去除标签后获取数据的所有特征。结果存储在allFeatures数据框中：

    > allFeatures <- allData[,-c(15)]

打印allFeatures数据框。head()函数返回allFeatures数据框的前部分。将allFeatures数据框作为输入参数传递：

    > head(allFeatures)

结果如下：

标准化特征。均值和尺度转换为z分数，使得variance = 1。scale()函数的默认方法将数值矩阵的列中心化和/或缩放。continuousFeatures是数值矩阵。结果存储在continuousFeatures数据框中：

    > continuousFeatures <- scale(continuousFeatures)

打印continuousFeatures数据框。head()函数返回continuousFeatures数据框的前部分。continuousFeatures数据框作为输入参数传递：

    > head(continuousFeatures)

结果如下：

将标签转换为1或-1。使用rep()函数复制值。结果存储在labels.n数据框中：

    > labels.n = rep(0,length(labels))     
> labels.n[labels==" <=50K"] = -1     
> labels.n[labels==" >50K"] = 1     
> labels = labels.n     
> rm(labels.n)

分离训练数据集。createDataPartition()函数创建一组训练数据分区。y=labels表示结果向量。p=.8表示 80%的数据用于训练数据集：

    > trainingData <- createDataPartition(y=labels, p=.8, list=FALSE)

探索trainingData数据框的维度。dim()函数返回trainingData数据框的维度。trainingData数据框作为输入参数传递。结果清楚地表明有 26,049 行数据和单列：

    > dim(trainingData)

结果如下：

创建trainingData数据框的训练特征和训练标签：

    > trainingFeatures <- continuousFeatures[trainingData,]     
> trainingLabels <- labels[trainingData]

确定剩余 20%的数据用于测试和验证：

    > remainingLabels <- labels[-trainingData]     
> remainingFeatures <- continuousFeatures[-trainingData,]

创建trainingData数据框的测试特征和测试标签。在 20%的数据中，其中 50%用于测试目的，剩余的 50%用于验证目的。

createDataPartition()函数创建一组训练数据分区。y=remainingLabels表示结果向量。p=.5表示 50%的数据用于训练数据集。结果存储在testingData数据框中：

    > testingData <- createDataPartition(y=remainingLabels, p=.5, list=FALSE)     
> testingLabels <- remainingLabels[testingData]     
> testingFeatures <- remainingFeatures[testingData,]

创建testingData数据框的验证特征和测试标签：

    > validationLabels <- remainingLabels[-testingData]
    > validationFeatures <- remainingFeatures[-testingData,]

定义所需的准确度度量：

> getAccuracy >- function(a,b,features,labels){
+ estFxn = features %*% a + b;
+ predictedLabels = rep(0,length(labels));
+ predictedLabels [estFxn < 0] = -1 ;
+ predictedLabels [estFxn >= 0] = 1 ;
+ return(sum(predictedLabels == labels) / length(labels))
+ }

第 4 步 - 构建模型

设置初始参数：

> numEpochs = 100
> numStepsPerEpoch = 500
> nStepsPerPlot = 30
> evalidationSetSize = 50
> c1 = 0.01
> c2 = 50

组合一组参数。结果存储在lambda_vals数据框中：

    > lambda_vals = c(0.001, 0.01, 0.1, 1)     
> bestAccuracy = 0

探索lambda_vals数据框的内部结构。str()函数显示数据框的内部结构。lambda_vals作为 R 对象传递给str()函数：

    > str(lambda_vals)

结果如下：

从给定的一组值中创建每个 epoch 的矩阵。使用matrix()函数。nrow = (numStepsPerEpoch/nStepsPerPlot)*numEpochs+1表示矩阵的行数，而ncol = length(lambda_vals)表示矩阵的列数：

    > accMat <- matrix(NA, nrow = (numStepsPerEpoch/nStepsPerPlot)*numEpochs+1, ncol = length(lambda_vals))

从给定的一组值中创建用于验证集准确性的矩阵。matrix() 函数被使用。nrow = (numStepsPerEpoch/nStepsPerPlot)*numEpochs+1 表示矩阵的行数，而 ncol = length(lambda_vals) 表示矩阵的列数：

    > accMatv <- matrix(NA, nrow = (numStepsPerEpoch/nStepsPerPlot)*numEpochs+1, ncol = length(lambda_vals))

设置分类器模型：

for(i in 1:4){ 
lambda = lambda_vals[i] 
accMatRow = 1 
accMatCol = i 
a = rep(0,ncol(continuousFeatures)) 
b = 0 
stepIndex = 0 
       for (e in 1:numEpochs){

#createDataPartition() 函数创建一组训练数据分区。y= trainingLabels 表示结果向量。p = (1 - evalidationSetSize/length(trainingLabels)) 百分比的数据用于训练数据集。结果存储在 etrainingData 数据框中：

etrainingData <- createDataPartition(y=trainingLabels, p=(1 -   evalidationSetSize/length(trainingLabels)), list=FALSE) 
 etrainingFeatures <- trainingFeatures[etrainingData,] 
 etrainingLabels <- trainingLabels[etrainingData] 
 evalidationFeatures <- trainingFeatures[-etrainingData,] 
 evalidationLabels <- trainingLabels[-etrainingData] 
 steplength = 1 / (e*c1 + c2) 
 for (step in 1:numStepsPerEpoch){ 
 stepIndex = stepIndex+1 
 index = sample.int(nrow(etrainingFeatures),1) 
 xk = etrainingFeatures[index,] 
 yk = etrainingLabels[index] 
 costfxn = yk * (a %*% xk + b) 
 if(costfxn >= 1){ 
 a_dir = lambda * a 
 a = a - steplength * a_dir 
 } else { 
 a_dir = (lambda * a) - (yk * xk) 
 a = a - steplength * a_dir 
 b_dir = -yk 
 b = b - (steplength * b_dir) 
 } 

记录准确性。调用 getAccuracy()：

if (stepIndex %% nStepsPerPlot == 1){#30){ 
accMat[accMatRow,accMatCol] = getAccuracy(a,b,evalidationFeatures,evalidationLabels) 
accMatv[accMatRow,accMatCol] = getAccuracy(a,b,validationFeatures,validationLabels) 
accMatRow = accMatRow + 1 
} 
} 
} 
tempAccuracy = getAccuracy(a,b,validationFeatures,validationLabels) 
print(str_c("tempAcc = ", tempAccuracy," and bestAcc = ", bestAccuracy) ) 
if(tempAccuracy > bestAccuracy){ 
bestAccuracy = tempAccuracy 
best_a = a 
best_b = b 
best_lambdaIndex = i 
} 
   }

计算模型的准确性。使用先前定义的 getAccuracy()：

   > getAccuracy(best_a,best_b, testingFeatures, testingLabels)

步骤 5 - 绘制模型

绘制训练过程中模型的准确性。使用 c() 函数将参数组合成向量：

    > colors = c("red","blue","green","black")

设置用于图表的向量：

> xaxislabel = "Step"
> yaxislabels = c("Accuracy on Randomized Epoch Validation
Set","Accuracy on Validation Set")
>
> ylims=c(0,1)
> stepValues = seq(1,15000,length=500)

创建一个通用向量。调用 list()，它将 accMat 和 accMatv 数据框连接起来：

    > mats =  list(accMat,accMatv)

绘制图表：

> for(j in 1:length(mats)){
mat = mats[[j]]
for(i in 1:4){
if(i == 1){

# plot() 函数是一个用于绘制 R 对象的通用函数。将 stepValues 数据框作为函数值传递：

 plot(stepValues, mat[1:500,i], type = "l",xlim=c(0, 15000), ylim=ylims, 
 col=colors[i],xlab=xaxislabel,ylab=yaxislabels[j],main=title) 
 } else{ 
 lines(stepValues, mat[1:500,i], type = "l",xlim=c(0, 15000), ylim=ylims, 
 col=colors[i],xlab=xaxislabel,ylab=yaxislabels[j],main=title) 
 } 
 Sys.sleep(1) 
 } 
 legend(x=10000,y=.5,legend=c("lambda=.001","lambda=.01","lambda=.1","lambda=1"),fill=colors) 
 } 

生成的图表将如下所示：