在Rstudio中安装本节课所有需要的包,运行以下指令:
install.packages(c(
"dplyr", "tidyr", "data.table",
"ggplot2", "patchwork", "RColorBrewer",
"caret", "mlr3", "tidymodels",
"MASS", "e1071", "class", "nnet",
"rpart", "rpart.plot", "randomForest", "gbm", "xgboost",
"pROC", "ROCR", "yardstick",
"recipes", "dimRed", "RANN",
"mlbench", "ISLR",
"purrr", "broom", "knitr", "kableExtra"
))
| 章节 | 主题 | 核心内容 |
|---|---|---|
| 第一章 | 机器学习基础概念 | 定义、学习类型、问题类型、核心概念 |
| 第二章 | 数据预处理与特征工程 | 特征缩放、编码、缺失值、特征选择 |
| 第三章 | 模型评估与验证 | 评估指标、混淆矩阵、ROC曲线、交叉验证 |
| 第四章 | 线性模型 | 线性回归、逻辑回归、正则化、SVM |
| 第五章 | 树模型与集成学习 | 决策树、随机森林、GBDT、XGBoost |
| 第六章 | 支持向量机与核方法 | 最大间隔、核函数、核技巧、SVR |
| 第七章 | 最近邻与基于实例的方法 | KNN、距离度量、加权KNN |
| 第八章 | 贝叶斯与概率图模型 | 朴素贝叶斯、贝叶斯网络、高斯过程 |
机器学习(Machine Learning) 是人工智能的一个分支,其核心思想是让计算机从数据中”学习”规律,而不是由人类显式编程所有规则。
经典定义:
美国计算机科学家Tom Mitchell给出的定义:“如果一个计算机程序在某类任务T上的性能P,随着经验E的积累而不断改善,则称该程序对任务T进行了学习。”
这个定义包含三个核心要素:
传统编程 vs 机器学习:
| 对比维度 | 传统编程 | 机器学习 |
|---|---|---|
| 输入 | 数据 + 规则 | 数据 + 答案 |
| 输出 | 答案 | 规则(模型) |
| 核心逻辑 | 人类编写规则 | 算法自动学习规则 |
| 适用场景 | 规则明确、可穷举 | 规则复杂、难以显式表达 |
机器学习本质上是一个函数逼近问题:
\[\hat{f} \approx f\]
其中: - \(f\):真实的数据生成函数(未知) - \(\hat{f}\):学习到的模型(已知)
我们的目标是从假设空间中找到一个最优的 \(\hat{f}\),使得它在未见数据上的表现尽可能好。
# 机器学习本质示例:拟合一个未知函数
# 假设真实函数 f(x) = sin(x) + 噪声
# 生成模拟数据
set.seed(123)
n <- 100
x <- seq(0, 2 * pi, length.out = n)
y_true <- sin(x) # 真实函数
y <- y_true + rnorm(n, mean = 0, sd = 0.2) # 观测值(带噪声)
# 创建数据框
data <- data.frame(x = x, y = y, y_true = y_true)
# 用多项式回归"学习"这个函数
# 尝试不同阶数的多项式
model_3 <- lm(y ~ poly(x, degree = 3), data = data)
model_9 <- lm(y ~ poly(x, degree = 9), data = data)
# 预测
data$pred_3 <- predict(model_3)
data$pred_9 <- predict(model_9)
# 可视化:展示学习过程
ggplot2::ggplot(data, ggplot2::aes(x = x)) +
ggplot2::geom_point(ggplot2::aes(y = y), alpha = 0.5, color = "gray50") +
ggplot2::geom_line(ggplot2::aes(y = y_true), color = "black", linewidth = 1, linetype = "dashed") +
ggplot2::geom_line(ggplot2::aes(y = pred_3, color = "3阶多项式"), linewidth = 1) +
ggplot2::geom_line(ggplot2::aes(y = pred_9, color = "9阶多项式"), linewidth = 1) +
ggplot2::scale_color_manual(values = c("3阶多项式" = "blue", "9阶多项式" = "red")) +
ggplot2::labs(
title = "机器学习本质:从数据中学习函数",
subtitle = "黑色虚线为真实函数,灰点为观测数据",
x = "x", y = "y", color = "模型"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
核心意义:机器学习让计算机能够从数据中自动发现规律,解决那些规则复杂、难以显式编程的问题。
举例:基于患者的临床指标预测糖尿病风险
为什么这种数据适合用机器学习?
糖尿病的发生与年龄、BMI、血糖、血压、家族史等多种因素相关,这些因素之间存在复杂的交互作用和非线性关系,难以用简单的规则描述。机器学习可以从大量历史病例数据中自动学习这些复杂关系,建立预测模型。
研究目的是什么?
可能得到什么样的结果?
模型可能发现:年龄大于45岁、BMI超过28、空腹血糖高于6.1 mmol/L的患者,糖尿病风险显著升高。这些规律可能不是简单的线性组合,而是复杂的非线性关系。
需要注意什么?
机器学习模型是”黑箱”,其预测结果需要结合医学知识进行解释。同时,模型在训练数据上的表现好,不代表在新患者上也能准确预测,需要验证模型的泛化能力。
机器学习并非万能,以下情况适合使用机器学习:
| 适用条件 | 说明 | 示例 |
|---|---|---|
| 规则复杂 | 问题涉及大量变量和复杂关系 | 疾病诊断、股票预测 |
| 数据充足 | 有足够的历史数据用于学习 | 医疗影像、交易记录 |
| 模式存在 | 数据中确实存在可学习的规律 | 垃圾邮件识别 |
| 容忍误差 | 允许一定程度的预测错误 | 推荐系统 |
| 重复任务 | 需要对大量相似任务做决策 | 信用评分 |
| 不适用条件 | 说明 | 替代方案 |
|---|---|---|
| 规则简单明确 | 可以用简单的if-else解决 | 传统编程 |
| 数据不足 | 样本量太少无法学习 | 专家系统 |
| 需要完美准确 | 零容忍错误场景 | 规则引擎 |
| 需要完全解释 | 必须解释每个决策原因 | 统计模型 |
| 实时性要求极高 | 毫秒级响应要求 | 优化算法 |
核心意义:判断何时使用机器学习是项目成功的第一步,避免在不适合的场景浪费资源。
举例:医院是否需要用机器学习预测患者住院时长
为什么需要先判断适用性?
医院信息系统积累了大量患者数据,管理者想用机器学习预测住院时长以优化床位管理。但需要先评估:住院时长受病情、治疗方案、并发症等多种因素影响,规则复杂;历史数据充足(数万条记录);确实存在可学习的规律(同类疾病住院时长有分布特征);允许一定预测误差(用于规划而非精确安排)。这些条件都满足,适合使用机器学习。
研究目的是什么?
可能得到什么样的结果?
模型可能发现:年龄大于65岁、有并发症、手术类型为某类的患者,平均住院时长显著延长。预测准确率可能达到70-80%,足以支持规划决策。
需要注意什么?
如果数据质量差(缺失值多、记录不准),或预测精度要求极高(如ICU床位必须精确到天),则不适合使用机器学习,应考虑其他方案。
机器学习根据学习方式和数据特点,主要分为四种范式:
| 学习类型 | 是否有标签 | 数据特点 | 典型任务 |
|---|---|---|---|
| 监督学习 | 有 | 带标签的训练数据 | 分类、回归 |
| 无监督学习 | 无 | 无标签数据 | 聚类、降维 |
| 半监督学习 | 部分 | 少量标签+大量无标签 | 分类 |
| 强化学习 | 延迟奖励 | 环境交互反馈 | 决策控制 |
定义:从带标签的训练数据中学习一个映射函数,将输入映射到输出。
数学表达:
给定训练集 \(D = \{(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\}\),学习函数 \(f: X \rightarrow Y\),使得对新输入 \(x_{new}\),能预测其输出 \(y_{new}\)。
核心特点: 1. 训练数据包含输入和对应的正确答案(标签) 2. 学习目标是建立输入到输出的映射关系 3. 可以用测试数据评估模型性能
# 监督学习示例:鸢尾花分类
# 加载鸢尾花数据集
data(iris)
# 查看数据结构
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
# 可视化:不同种类的鸢尾花在特征空间中的分布
ggplot2::ggplot(iris, ggplot2::aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
ggplot2::geom_point(size = 3, alpha = 0.7) +
ggplot2::labs(
title = "监督学习示例:鸢尾花分类",
subtitle = "已知类别标签,学习分类边界",
x = "萼片长度", y = "花瓣长度", color = "种类"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
定义:从无标签数据中发现隐藏的结构或模式。
数学表达:
给定数据集 \(D = \{x_1, x_2, ..., x_n\}\)(无标签),发现数据的内在结构,如聚类或低维表示。
核心特点: 1. 训练数据没有标签 2. 学习目标是发现数据的内在结构 3. 评估标准不如监督学习明确
# 无监督学习示例:K-means聚类
# 使用鸢尾花数据,但假装不知道真实类别
# 只取特征列(去掉标签)
iris_features <- iris[, 1:4]
# K-means聚类(假设不知道有3类)
set.seed(42)
kmeans_result <- kmeans(iris_features, centers = 3, nstart = 10)
# 将聚类结果添加到数据框
iris_clustered <- iris
iris_clustered$cluster <- as.factor(kmeans_result$cluster)
# 可视化:聚类结果 vs 真实类别
p1 <- ggplot2::ggplot(iris_clustered, ggplot2::aes(x = Sepal.Length, y = Petal.Length, color = cluster)) +
ggplot2::geom_point(size = 3, alpha = 0.7) +
ggplot2::labs(title = "无监督学习:K-means聚类结果", x = "萼片长度", y = "花瓣长度", color = "聚类") +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
p2 <- ggplot2::ggplot(iris, ggplot2::aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
ggplot2::geom_point(size = 3, alpha = 0.7) +
ggplot2::labs(title = "真实类别标签", x = "萼片长度", y = "花瓣长度", color = "种类") +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
p1 + p2 + patchwork::plot_annotation(
title = "无监督学习 vs 真实标签对比",
subtitle = "无监督学习在没有标签的情况下发现数据结构"
)
定义:结合少量有标签数据和大量无标签数据进行学习。
适用场景: 1. 标签获取成本高(需要专家标注) 2. 无标签数据容易获取 3. 医学影像诊断(标注需要专业医生)
核心思想: 1. 利用有标签数据学习初步模型 2. 利用无标签数据增强模型(如伪标签、一致性正则化)
# 半监督学习示意图
# 假设只有10%的数据有标签
set.seed(123)
n_labeled <- 15 # 只有15个样本有标签
# 随机选择有标签的样本
labeled_idx <- sample(1:nrow(iris), n_labeled)
iris_semi <- iris
iris_semi$label_status <- ifelse(1:nrow(iris) %in% labeled_idx, "有标签", "无标签")
iris_semi$known_species <- ifelse(1:nrow(iris) %in% labeled_idx, as.character(iris$Species), "未知")
# 可视化半监督学习场景
ggplot2::ggplot(iris_semi, ggplot2::aes(x = Sepal.Length, y = Petal.Length)) +
ggplot2::geom_point(ggplot2::aes(color = known_species, shape = label_status), size = 3, alpha = 0.7) +
ggplot2::scale_shape_manual(values = c("有标签" = 16, "无标签" = 1)) +
ggplot2::labs(
title = "半监督学习场景示意",
subtitle = paste0("只有", n_labeled, "个样本有标签(实心点),其余无标签(空心点)"),
x = "萼片长度", y = "花瓣长度", color = "已知类别", shape = "标签状态"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
定义:智能体通过与环境交互,根据奖励/惩罚信号学习最优策略。
核心要素: 1. 智能体(Agent):学习者 2. 环境(Environment):智能体交互的对象 3. 动作(Action):智能体的行为 4. 状态(State):环境的当前情况 5. 奖励(Reward):环境对动作的反馈
数学框架:马尔可夫决策过程(MDP)
\[\text{目标:最大化累积奖励} \quad G_t = \sum_{k=0}^{\infty} \gamma^k R_{t+k+1}\]
其中 \(\gamma\) 是折扣因子,\(0 \leq \gamma < 1\)。
# 强化学习概念示意图
# 创建一个简单的网格世界示例
grid_world <- data.frame(
x = rep(1:4, each = 4),
y = rep(1:4, 4),
cell_type = c(
"起点", "普通", "普通", "普通",
"普通", "障碍", "普通", "普通",
"普通", "普通", "普通", "陷阱",
"普通", "普通", "普通", "终点"
),
reward = c(
0, -1, -1, -1,
-1, NA, -1, -1,
-1, -1, -1, -100,
-1, -1, -1, 100
)
)
# 添加标签
grid_world$label <- ifelse(is.na(grid_world$reward), "障碍",
ifelse(grid_world$reward == 100, "终点(+100)",
ifelse(grid_world$reward == -100, "陷阱(-100)",
ifelse(grid_world$reward == 0, "起点", ""))))
ggplot2::ggplot(grid_world, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = cell_type), color = "white", linewidth = 1) +
ggplot2::geom_text(ggplot2::aes(label = label), size = 3) +
ggplot2::scale_fill_manual(values = c(
"起点" = "green", "终点" = "gold", "障碍" = "gray50",
"陷阱" = "red", "普通" = "lightblue"
)) +
ggplot2::labs(
title = "强化学习示例:网格世界",
subtitle = "智能体从起点出发,学习避开陷阱、到达终点的最优路径",
x = "", y = "", fill = "格子类型"
) +
ggplot2::theme_minimal() +
ggplot2::theme(
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "bottom"
) +
ggplot2::coord_equal()
机器学习任务根据输出类型的不同,主要分为三大类:
| 问题类型 | 输出类型 | 目标 | 典型算法 |
|---|---|---|---|
| 分类问题 | 离散类别 | 预测样本所属类别 | 逻辑回归、决策树、SVM |
| 回归问题 | 连续数值 | 预测连续值 | 线性回归、随机森林回归 |
| 聚类问题 | 无(发现结构) | 将相似样本分组 | K-means、层次聚类 |
定义:预测样本属于哪个离散类别。
数学表达:
给定输入 \(x\),预测其类别 \(y \in \{C_1, C_2, ..., C_K\}\)
分类类型: 1. 二分类:只有两个类别(如:患病/健康) 2. 多分类:多个类别(如:肿瘤类型I/II/III期)
# 分类问题示例:乳腺癌诊断
# 使用威斯康星乳腺癌数据集
data(BreastCancer, package = "mlbench")
# 数据预处理
bc_data <- BreastCancer %>%
dplyr::select(-Id) %>%
dplyr::mutate_all(~ifelse(. == "?", NA, .)) %>%
tidyr::drop_na() %>%
dplyr::mutate_at(vars(-Class), as.numeric)
# 查看数据
head(bc_data[, 1:5])
## Cl.thickness Cell.size Cell.shape Marg.adhesion Epith.c.size
## 1 5 1 1 1 2
## 2 5 4 4 5 7
## 3 3 1 1 1 2
## 4 6 8 8 1 3
## 5 4 1 1 3 2
## 6 8 10 10 8 7
# 可视化:两类样本在特征空间中的分布
ggplot2::ggplot(bc_data, ggplot2::aes(x = Cl.thickness, y = Cell.size, color = Class)) +
ggplot2::geom_jitter(size = 2, alpha = 0.5, width = 0.2, height = 0.2) +
ggplot2::labs(
title = "分类问题示例:乳腺癌诊断",
subtitle = "根据细胞特征预测良性/恶性",
x = "细胞厚度", y = "细胞大小", color = "诊断结果"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
定义:预测连续的数值输出。
数学表达:
给定输入 \(x\),预测连续输出 \(y \in \mathbb{R}\)
与分类的区别: - 分类输出是离散的类别标签 - 回归输出是连续的数值
# 回归问题示例:医疗费用预测
# 使用模拟数据演示
set.seed(456)
n_patients <- 200
# 生成模拟数据
medical_cost <- data.frame(
age = round(runif(n_patients, 18, 80)),
bmi = round(rnorm(n_patients, 25, 4), 1),
smoker = sample(c("Yes", "No"), n_patients, replace = TRUE, prob = c(0.3, 0.7))
) %>%
dplyr::mutate(
# 医疗费用与年龄、BMI、吸烟状态相关
cost = 1000 + 50 * age + 100 * bmi +
ifelse(smoker == "Yes", 5000, 0) +
rnorm(n_patients, 0, 1000)
)
# 查看数据
head(medical_cost)
## age bmi smoker cost
## 1 24 25.5 No 5467.349
## 2 31 28.5 Yes 9982.663
## 3 63 24.6 No 6578.476
## 4 71 25.3 Yes 11175.188
## 5 67 18.3 No 6179.478
## 6 39 29.5 Yes 12297.854
# 可视化:费用与年龄的关系
ggplot2::ggplot(medical_cost, ggplot2::aes(x = age, y = cost, color = smoker)) +
ggplot2::geom_point(size = 2, alpha = 0.6) +
ggplot2::geom_smooth(method = "lm", se = FALSE) +
ggplot2::labs(
title = "回归问题示例:医疗费用预测",
subtitle = "预测连续的费用数值(元)",
x = "年龄(岁)", y = "医疗费用(元)", color = "是否吸烟"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
定义:将相似的样本归为同一组,发现数据的内在结构。
特点: 1. 无监督学习(无标签) 2. 目标是发现数据的自然分组 3. “正确”的聚类结果可能不唯一
# 聚类问题示例:患者分型
# 使用模拟数据演示
set.seed(789)
n_patients <- 150
# 生成三种不同类型的患者群体
patient_types <- data.frame(
# 年轻低风险组
age = c(rnorm(50, 35, 5), rnorm(50, 55, 8), rnorm(50, 70, 6)),
risk_score = c(rnorm(50, 20, 5), rnorm(50, 50, 8), rnorm(50, 75, 7)),
true_group = rep(c("低风险组", "中风险组", "高风险组"), each = 50)
)
# 进行K-means聚类(假装不知道真实分组)
kmeans_patients <- kmeans(patient_types[, c("age", "risk_score")], centers = 3, nstart = 10)
patient_types$cluster <- as.factor(kmeans_patients$cluster)
# 可视化聚类结果
ggplot2::ggplot(patient_types, ggplot2::aes(x = age, y = risk_score)) +
ggplot2::geom_point(ggplot2::aes(color = cluster), size = 3, alpha = 0.7) +
ggplot2::labs(
title = "聚类问题示例:患者风险分型",
subtitle = "无监督地将患者分为不同风险群体",
x = "年龄(岁)", y = "风险评分", color = "聚类结果"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
| 问题类型 | 输出类型 | 目标 | 典型算法 | 评估指标 | 医学示例 |
|---|---|---|---|---|---|
| 分类问题 | 离散类别 | 预测样本所属类别 | 逻辑回归、决策树、SVM | 准确率、F1分数 | 疾病诊断、病理分型 |
| 回归问题 | 连续数值 | 预测连续值 | 线性回归、随机森林回归 | MSE、R² | 费用预测、生存时间预测 |
| 聚类问题 | 无(发现结构) | 将相似样本分组 | K-means、层次聚类 | 轮廓系数、SSE | 患者分型、基因表达聚类 |
核心意义:正确识别问题类型是选择合适算法的前提,不同问题类型对应不同的建模方法和评估指标。
举例:心血管疾病风险预测项目
为什么需要区分问题类型?
同一个医学问题可能被建模为不同类型:预测”是否患病”是分类问题,预测”患病概率”仍是分类(但输出概率),预测”生存年数”是回归问题,发现”患者亚群”是聚类问题。明确问题类型后才能选择正确的方法。
研究目的是什么?
可能得到什么样的结果?
分类模型可能输出”高风险/低风险”标签,准确率85%;回归模型可能输出”5年生存概率78%“,R²=0.65;聚类模型可能发现3个患者亚群,各有不同的预后特征。
需要注意什么?
有时问题类型可以相互转化:回归问题可以通过设置阈值转为分类问题;分类问题的概率输出可以用于排序。选择哪种类型取决于业务需求和数据特点。
定义:描述样本的可测量属性或特性。
特征类型:
| 特征类型 | 数据类型 | 示例 | 编码方式 |
|---|---|---|---|
| 数值特征 | 连续/离散 | 年龄、血压 | 直接使用 |
| 类别特征 | 名义/有序 | 性别、疾病分期 | 独热编码/标签编码 |
| 文本特征 | 字符串 | 病历描述 | TF-IDF/词嵌入 |
| 图像特征 | 像素矩阵 | X光片 | CNN提取 |
# 特征示例:心脏病预测数据集
# 创建模拟的心脏病预测数据
set.seed(111)
n <- 100
heart_data <- data.frame(
# 数值特征
age = round(runif(n, 30, 80)), # 年龄(连续)
resting_bp = round(rnorm(n, 130, 20)), # 静息血压(连续)
cholesterol = round(rnorm(n, 200, 40)), # 胆固醇(连续)
max_heart_rate = round(220 - runif(n, 30, 80) + rnorm(n, 0, 10)), # 最大心率
# 类别特征
sex = sample(c("Male", "Female"), n, replace = TRUE), # 性别(名义)
chest_pain = sample(c("Typical", "Atypical", "Non-anginal", "Asymptomatic"),
n, replace = TRUE), # 胸痛类型(名义)
fasting_blood_sugar = sample(c(">120", "<=120"), n, replace = TRUE), # 空腹血糖(二值)
# 标签
heart_disease = sample(c("Yes", "No"), n, replace = TRUE, prob = c(0.4, 0.6))
)
# 查看数据结构
head(heart_data)
## age resting_bp cholesterol max_heart_rate sex chest_pain
## 1 60 134 196 183 Male Asymptomatic
## 2 66 161 125 153 Male Non-anginal
## 3 49 148 173 147 Female Typical
## 4 56 137 208 157 Male Asymptomatic
## 5 49 134 96 162 Male Asymptomatic
## 6 51 113 196 154 Female Non-anginal
## fasting_blood_sugar heart_disease
## 1 >120 No
## 2 >120 No
## 3 >120 No
## 4 >120 Yes
## 5 >120 No
## 6 <=120 No
# 特征类型总结
| 特征名称 | 特征类型 | 数据类型 | 取值范围 |
|---|---|---|---|
| age | 数值-连续 | 整数 | 30-80岁 |
| resting_bp | 数值-连续 | 整数 | 约80-180 mmHg |
| cholesterol | 数值-连续 | 整数 | 约100-300 mg/dL |
| sex | 类别-名义 | 字符 | Male/Female |
| chest_pain | 类别-名义 | 字符 | 4种类型 |
| fasting_blood_sugar | 类别-二值 | 字符 | >120/<=120 |
定义:监督学习中样本的目标值或正确答案。
标签类型:
| 任务类型 | 标签类型 | 示例 |
|---|---|---|
| 二分类 | 二值 | 患病/健康 |
| 多分类 | 多个离散值 | 肿瘤I/II/III期 |
| 回归 | 连续数值 | 医疗费用、生存时间 |
# 标签分布可视化
ggplot2::ggplot(heart_data, ggplot2::aes(x = heart_disease, fill = heart_disease)) +
ggplot2::geom_bar(width = 0.6) +
ggplot2::geom_text(stat = "count", ggplot2::aes(label = after_stat(count)), vjust = -0.5) +
ggplot2::labs(
title = "标签分布:心脏病诊断结果",
x = "诊断结果", y = "样本数量", fill = "诊断结果"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "none")
在机器学习中,数据通常表示为:
\[X = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{bmatrix}, \quad y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}\]
# 特征矩阵与标签向量的表示
# 提取数值特征作为特征矩阵
X <- heart_data %>%
dplyr::select(age, resting_bp, cholesterol, max_heart_rate) %>%
as.matrix()
# 提取标签向量
y <- heart_data$heart_disease
# 显示维度
cat("特征矩阵 X 的维度:", dim(X), "\n")
## 特征矩阵 X 的维度: 100 4
cat(" - 样本数量 n =", nrow(X), "\n")
## - 样本数量 n = 100
cat(" - 特征数量 p =", ncol(X), "\n")
## - 特征数量 p = 4
cat("\n标签向量 y 的长度:", length(y), "\n")
##
## 标签向量 y 的长度: 100
cat("标签类别:", unique(y), "\n")
## 标签类别: No Yes
# 特征矩阵前几行
cat("\n特征矩阵 X 前5行:\n")
##
## 特征矩阵 X 前5行:
print(head(X, 5))
## age resting_bp cholesterol max_heart_rate
## [1,] 60 134 196 183
## [2,] 66 161 125 153
## [3,] 49 148 173 147
## [4,] 56 137 208 157
## [5,] 49 134 96 162
核心意义:特征是机器学习模型的”原材料”,标签是监督学习的”答案”。特征的质量直接决定模型的上限。
举例:基于电子病历预测患者再入院风险
为什么特征选择很重要?
再入院风险受多种因素影响:年龄、诊断、用药、既往史、社会因素等。选择哪些作为特征,决定了模型能学到什么。如果遗漏重要特征(如既往再入院史),模型性能会受限;如果包含无关特征(如病历编号),可能引入噪声。
研究目的是什么?
可能得到什么样的结果?
模型可能发现:年龄大于65岁、有多次既往住院史、出院时用药超过5种的患者,再入院风险显著升高。这些特征组合可以生成风险评分。
需要注意什么?
特征工程是机器学习中最耗时的环节之一。需要注意:1. 特征的完整性和准确性;2. 特征之间的相关性;3. 特征的可解释性;4. 避免数据泄露(如使用了未来信息)。
在机器学习中,数据通常划分为三个独立的部分:
| 数据集 | 作用 | 比例建议 | 使用场景 |
|---|---|---|---|
| 训练集 | 训练模型参数 | 60-80% | 模型学习 |
| 验证集 | 调整超参数、模型选择 | 10-20% | 模型调优 |
| 测试集 | 评估最终性能 | 10-20% | 最终评估 |
核心问题:模型在训练数据上的表现好,不代表在新数据上也好。
目标:评估模型的泛化能力——在未见数据上的表现。
三种数据集的分工:
# 数据划分示例
# 使用iris数据集演示
set.seed(42)
# 获取样本数量
n <- nrow(iris)
# 随机打乱样本顺序
indices <- sample(1:n)
# 计算划分点
train_end <- round(0.7 * n)
val_end <- round(0.85 * n)
# 划分数据集
train_idx <- indices[1:train_end]
val_idx <- indices[(train_end + 1):val_end]
test_idx <- indices[(val_end + 1):n]
# 创建三个数据集
train_data <- iris[train_idx, ]
val_data <- iris[val_idx, ]
test_data <- iris[test_idx, ]
# 显示划分结果
cat("训练集样本数:", nrow(train_data), "\n")
## 训练集样本数: 105
cat("验证集样本数:", nrow(val_data), "\n")
## 验证集样本数: 23
cat("测试集样本数:", nrow(test_data), "\n")
## 测试集样本数: 22
# 可视化数据划分
split_viz <- data.frame(
dataset = rep(c("训练集", "验证集", "测试集"),
times = c(nrow(train_data), nrow(val_data), nrow(test_data))),
index = c(train_idx, val_idx, test_idx)
)
split_viz$dataset <- factor(split_viz$dataset, levels = c("训练集", "验证集", "测试集"))
ggplot2::ggplot(split_viz, ggplot2::aes(x = index, y = 1, color = dataset)) +
ggplot2::geom_point(size = 2, alpha = 0.7) +
ggplot2::labs(
title = "数据集划分可视化",
subtitle = "随机打乱后按比例划分",
x = "样本索引", y = "", color = "数据集"
) +
ggplot2::theme_minimal() +
ggplot2::theme(
axis.text.y = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
legend.position = "bottom"
)
# 使用caret包进行更规范的数据划分
library(caret)
# 创建训练集和测试集(70%-30%划分)
set.seed(42)
train_index <- caret::createDataPartition(iris$Species, p = 0.7, list = FALSE)
train_set <- iris[train_index, ]
test_set <- iris[-train_index, ]
# 在训练集内部再划分出验证集
set.seed(42)
val_index <- caret::createDataPartition(train_set$Species, p = 0.8, list = FALSE)
train_final <- train_set[val_index, ]
val_set <- train_set[-val_index, ]
# 显示最终划分
cat("训练集样本数:", nrow(train_final), "\n")
## 训练集样本数: 84
cat("验证集样本数:", nrow(val_set), "\n")
## 验证集样本数: 21
cat("测试集样本数:", nrow(test_set), "\n")
## 测试集样本数: 45
当类别不平衡时,需要使用分层抽样确保各数据集中类别比例一致。
# 分层抽样示例
# 创建不平衡数据集
set.seed(123)
imbalanced_data <- data.frame(
feature = rnorm(1000),
label = c(rep("Positive", 100), rep("Negative", 900)) # 1:9不平衡
)
# 普通随机划分 vs 分层划分
set.seed(42)
# 普通随机划分
random_split <- sample(1:nrow(imbalanced_data), 0.7 * nrow(imbalanced_data))
train_random <- imbalanced_data[random_split, ]
# 分层划分(使用caret)
stratified_index <- caret::createDataPartition(imbalanced_data$label, p = 0.7, list = FALSE)
train_stratified <- imbalanced_data[stratified_index, ]
# 比较两种划分方式的类别比例
cat("原始数据Positive比例:", sum(imbalanced_data$label == "Positive") / nrow(imbalanced_data) * 100, "%\n")
## 原始数据Positive比例: 10 %
cat("随机划分Positive比例:", sum(train_random$label == "Positive") / nrow(train_random) * 100, "%\n")
## 随机划分Positive比例: 10.14286 %
cat("分层划分Positive比例:", sum(train_stratified$label == "Positive") / nrow(train_stratified) * 100, "%\n")
## 分层划分Positive比例: 10 %
核心意义:合理的数据划分是评估模型泛化能力的基础,避免”用考试题训练”的作弊行为。
举例:构建糖尿病预测模型
为什么要划分三个数据集?
假设我们有1000名患者的数据。如果全部用于训练,模型可能”记住”了这些患者的特征,在新患者上表现很差。划分后:700人用于训练,150人用于调参,150人用于最终测试。测试集的150人模拟”新患者”,评估模型的真实泛化能力。
研究目的是什么?
可能得到什么样的结果?
训练集准确率85%,验证集准确率82%,测试集准确率80%。如果测试集准确率远低于训练集,说明模型过拟合,需要简化模型或增加数据。
需要注意什么?
过拟合(Overfitting):模型在训练数据上表现很好,但在新数据上表现差。模型”记住”了训练数据的噪声和细节,而非真正的规律。
欠拟合(Underfitting):模型在训练数据和新数据上都表现不好。模型太简单,无法捕捉数据中的规律。
# 过拟合与欠拟合示意图
# 使用多项式回归演示
set.seed(42)
n <- 30
# 生成真实数据(带噪声的二次函数)
x <- seq(-3, 3, length.out = n)
y_true <- 0.5 * x^2 - 1
y <- y_true + rnorm(n, sd = 0.5)
# 创建数据框
poly_data <- data.frame(x = x, y = y, y_true = y_true)
# 拟合不同复杂度的模型
# 欠拟合:线性模型(太简单)
model_underfit <- lm(y ~ x, data = poly_data)
# 适度拟合:二次多项式
model_good <- lm(y ~ poly(x, degree = 2, raw = TRUE), data = poly_data)
# 过拟合:高阶多项式(太复杂)
model_overfit <- lm(y ~ poly(x, degree = 15, raw = TRUE), data = poly_data)
# 预测
poly_data$pred_underfit <- predict(model_underfit)
poly_data$pred_good <- predict(model_good)
poly_data$pred_overfit <- predict(model_overfit)
# 可视化
p1 <- ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(size = 2, alpha = 0.7) +
ggplot2::geom_line(ggplot2::aes(y = pred_underfit), color = "blue", linewidth = 1) +
ggplot2::labs(title = "欠拟合(线性模型)", subtitle = "模型太简单,无法捕捉规律", x = "", y = "") +
ggplot2::theme_minimal() +
ggplot2::ylim(-3, 5)
p2 <- ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(size = 2, alpha = 0.7) +
ggplot2::geom_line(ggplot2::aes(y = pred_good), color = "green", linewidth = 1) +
ggplot2::labs(title = "适度拟合(二次多项式)", subtitle = "模型复杂度适中", x = "", y = "") +
ggplot2::theme_minimal() +
ggplot2::ylim(-3, 5)
p3 <- ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(size = 2, alpha = 0.7) +
ggplot2::geom_line(ggplot2::aes(y = pred_overfit), color = "red", linewidth = 1) +
ggplot2::labs(title = "过拟合(15阶多项式)", subtitle = "模型太复杂,拟合了噪声", x = "", y = "") +
ggplot2::theme_minimal() +
ggplot2::ylim(-3, 5)
p1 + p2 + p3 + patchwork::plot_annotation(
title = "欠拟合 vs 适度拟合 vs 过拟合",
subtitle = "蓝点为训练数据,曲线为模型预测"
)
| 特征 | 欠拟合 | 过拟合 |
|---|---|---|
| 训练误差 | 高 | 低 |
| 测试误差 | 高 | 高 |
| 原因 | 模型太简单 | 模型太复杂 |
| 表现 | 学不到规律 | 记住噪声 |
| 解决方法 | 增加模型复杂度 | 正则化、增加数据、简化模型 |
# 计算训练误差和测试误差
# 生成新数据作为测试集
set.seed(999)
x_test <- seq(-3, 3, length.out = 100)
y_test_true <- 0.5 * x_test^2 - 1
y_test <- y_test_true + rnorm(100, sd = 0.5)
test_data <- data.frame(x = x_test, y = y_test)
# 预测测试集
test_data$pred_underfit <- predict(model_underfit, newdata = test_data)
test_data$pred_good <- predict(model_good, newdata = test_data)
test_data$pred_overfit <- predict(model_overfit, newdata = test_data)
# 计算均方误差
calc_mse <- function(actual, predicted) {
mean((actual - predicted)^2)
}
cat("欠拟合(线性)- 训练MSE:", calc_mse(poly_data$y, poly_data$pred_underfit), "测试MSE:", calc_mse(test_data$y, test_data$pred_underfit), "\n")
## 欠拟合(线性)- 训练MSE: 2.525686 测试MSE: 2.151376
cat("适度拟合(二次)- 训练MSE:", calc_mse(poly_data$y, poly_data$pred_good), "测试MSE:", calc_mse(test_data$y, test_data$pred_good), "\n")
## 适度拟合(二次)- 训练MSE: 0.3484976 测试MSE: 0.2811028
cat("过拟合(15阶)- 训练MSE:", calc_mse(poly_data$y, poly_data$pred_overfit), "测试MSE:", calc_mse(test_data$y, test_data$pred_overfit), "\n")
## 过拟合(15阶)- 训练MSE: 0.2188181 测试MSE: 0.4459991
# 模型复杂度与误差的关系图
# 模拟不同复杂度下的训练误差和测试误差
complexity <- 1:15
train_errors <- c()
test_errors <- c()
for (deg in complexity) {
model <- lm(y ~ poly(x, degree = deg, raw = TRUE), data = poly_data)
train_errors <- c(train_errors, calc_mse(poly_data$y, predict(model)))
test_errors <- c(test_errors, calc_mse(test_data$y, predict(model, newdata = test_data)))
}
error_data <- data.frame(
complexity = rep(complexity, 2),
error = c(train_errors, test_errors),
type = rep(c("训练误差", "测试误差"), each = length(complexity))
)
ggplot2::ggplot(error_data, ggplot2::aes(x = complexity, y = error, color = type)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::geom_point(size = 2) +
ggplot2::geom_vline(xintercept = 2, linetype = "dashed", color = "gray50") +
ggplot2::annotate("text", x = 2.5, y = max(error_data$error) * 0.9,
label = "最优复杂度", hjust = 0) +
ggplot2::labs(
title = "模型复杂度与误差的关系",
subtitle = "训练误差持续下降,测试误差先降后升(过拟合)",
x = "模型复杂度(多项式阶数)", y = "均方误差(MSE)", color = ""
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
期望预测误差可以分解为三部分:
\[E[(y - \hat{f}(x))^2] = \underbrace{Bias[\hat{f}(x)]^2}_{偏差} + \underbrace{Var[\hat{f}(x)]}_{方差} + \underbrace{\sigma^2}_{不可约误差}\]
其中: - 偏差(Bias):模型预测的期望值与真实值的差异 \[Bias[\hat{f}(x)] = E[\hat{f}(x)] - f(x)\] - 方差(Variance):模型预测的波动程度 \[Var[\hat{f}(x)] = E[(\hat{f}(x) - E[\hat{f}(x)])^2]\] - 不可约误差(Irreducible Error):数据本身的噪声,无法消除
| 概念 | 含义 | 高偏差表现 | 高方差表现 |
|---|---|---|---|
| 偏差 | 预测的准确性 | 预测值偏离真实值 | - |
| 方差 | 预测的稳定性 | - | 不同训练集得到差异大的模型 |
靶心类比:
# 偏差-方差权衡示意图
# 使用靶心图演示
set.seed(42)
# 创建四种情况的模拟数据
n_points <- 50
# 低偏差低方差
low_bias_low_var <- data.frame(
x = rnorm(n_points, mean = 0, sd = 0.3),
y = rnorm(n_points, mean = 0, sd = 0.3),
type = "低偏差 低方差\n(理想)"
)
# 高偏差低方差
high_bias_low_var <- data.frame(
x = rnorm(n_points, mean = 2, sd = 0.3),
y = rnorm(n_points, mean = 2, sd = 0.3),
type = "高偏差 低方差\n(欠拟合)"
)
# 低偏差高方差
low_bias_high_var <- data.frame(
x = rnorm(n_points, mean = 0, sd = 1.5),
y = rnorm(n_points, mean = 0, sd = 1.5),
type = "低偏差 高方差\n(过拟合)"
)
# 高偏差高方差
high_bias_high_var <- data.frame(
x = rnorm(n_points, mean = 2, sd = 1.5),
y = rnorm(n_points, mean = 2, sd = 1.5),
type = "高偏差 高方差\n(最差)"
)
# 合并数据
target_data <- rbind(low_bias_low_var, high_bias_low_var,
low_bias_high_var, high_bias_high_var)
target_data$type <- factor(target_data$type, levels = unique(target_data$type))
# 绘制靶心图
ggplot2::ggplot(target_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(alpha = 0.6, size = 2) +
ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
ggplot2::annotate("point", x = 0, y = 0, shape = 21, size = 5,
fill = "red", color = "black", stroke = 1) +
ggplot2::facet_wrap(~ type) +
ggplot2::labs(
title = "偏差-方差权衡:靶心图",
subtitle = "红点为靶心(真实值),蓝点为预测值",
x = "", y = ""
) +
ggplot2::theme_minimal() +
ggplot2::coord_fixed(xlim = c(-4, 4), ylim = c(-4, 4))
# 模型复杂度与偏差/方差的关系
complexity <- 1:20
bias_squared <- 1 / complexity + 0.1 # 偏差随复杂度降低
variance <- complexity / 20 # 方差随复杂度增加
total_error <- bias_squared + variance + 0.1 # 总误差
bv_data <- data.frame(
complexity = rep(complexity, 3),
value = c(bias_squared, variance, total_error),
component = rep(c("偏差²", "方差", "总误差"), each = length(complexity))
)
ggplot2::ggplot(bv_data, ggplot2::aes(x = complexity, y = value, color = component)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::geom_point(size = 2) +
ggplot2::geom_vline(xintercept = which.min(total_error), linetype = "dashed", color = "gray50") +
ggplot2::annotate("text", x = which.min(total_error) + 0.5, y = max(bv_data$value) * 0.9,
label = "最优复杂度", hjust = 0) +
ggplot2::labs(
title = "偏差-方差权衡",
subtitle = "总误差 = 偏差² + 方差 + 不可约误差",
x = "模型复杂度", y = "误差", color = "组成部分"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
泛化能力(Generalization):模型在未见数据上的表现能力。
泛化误差(Generalization Error):模型在新数据上的期望误差。
\[\text{泛化误差} = E_{(x,y) \sim P} [L(y, \hat{f}(x))]\]
其中 \(P\) 是数据的真实分布,\(L\) 是损失函数。
训练误差 vs 泛化误差:
目标:最小化泛化误差,而非训练误差。
# 泛化能力示意图
# 使用不同大小的训练集演示
set.seed(42)
# 生成数据
n_total <- 1000
x_all <- runif(n_total, 0, 10)
y_all <- 2 * x_all + 1 + rnorm(n_total, sd = 2) # 真实关系: y = 2x + 1 + 噪声
all_data <- data.frame(x = x_all, y = y_all)
# 测试集(固定)
test_idx <- sample(1:n_total, 200)
test_data <- all_data[test_idx, ]
train_pool <- all_data[-test_idx, ]
# 不同训练集大小下的训练误差和测试误差
train_sizes <- c(20, 50, 100, 200, 400, 600, 800)
results <- data.frame()
for (n_train in train_sizes) {
# 多次实验取平均
train_errors <- c()
test_errors <- c()
for (i in 1:10) {
idx <- sample(1:nrow(train_pool), n_train)
train_data <- train_pool[idx, ]
model <- lm(y ~ x, data = train_data)
train_errors <- c(train_errors, mean((train_data$y - predict(model, train_data))^2))
test_errors <- c(test_errors, mean((test_data$y - predict(model, test_data))^2))
}
results <- rbind(results, data.frame(
train_size = n_train,
train_error = mean(train_errors),
test_error = mean(test_errors)
))
}
# 可视化
results_long <- tidyr::pivot_longer(results, cols = c(train_error, test_error),
names_to = "type", values_to = "error")
ggplot2::ggplot(results_long, ggplot2::aes(x = train_size, y = error, color = type)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::geom_point(size = 2) +
ggplot2::labs(
title = "训练数据量与泛化能力的关系",
subtitle = "随着训练数据增加,训练误差和测试误差趋于接近",
x = "训练样本数量", y = "均方误差(MSE)", color = ""
) +
ggplot2::scale_color_discrete(labels = c("训练误差", "测试误差(泛化误差估计)")) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
| 方法 | 原理 | 适用场景 |
|---|---|---|
| 增加训练数据 | 减少过拟合 | 数据不足时 |
| 正则化 | 限制模型复杂度 | 过拟合时 |
| 交叉验证 | 更准确估计泛化误差 | 模型选择时 |
| 特征选择 | 减少噪声特征 | 特征过多时 |
| 集成学习 | 降低方差 | 单模型不稳定时 |
# 学习曲线示例
# 使用caret包绘制学习曲线
set.seed(42)
# 加载数据
data(Sonar, package = "mlbench")
# 定义训练控制
train_control <- caret::trainControl(
method = "cv",
number = 5
)
# 计算不同训练集大小下的性能
train_sizes <- seq(0.1, 1, by = 0.1)
n <- nrow(Sonar)
results <- data.frame()
for (prop in train_sizes) {
n_train <- round(prop * n)
# 多次采样取平均
accs <- c()
for (i in 1:5) {
idx <- sample(1:n, n_train)
train_subset <- Sonar[idx, ]
model <- caret::train(
Class ~ .,
data = train_subset,
method = "rpart",
trControl = train_control,
tuneLength = 3
)
accs <- c(accs, max(model$results$Accuracy))
}
results <- rbind(results, data.frame(
train_size = n_train,
accuracy = mean(accs),
sd = sd(accs)
))
}
# 绘制学习曲线
ggplot2::ggplot(results, ggplot2::aes(x = train_size, y = accuracy)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_point(color = "steelblue", size = 2) +
ggplot2::geom_ribbon(ggplot2::aes(ymin = accuracy - sd, ymax = accuracy + sd),
alpha = 0.2, fill = "steelblue") +
ggplot2::labs(
title = "学习曲线示例",
subtitle = "模型性能随训练数据量增加而提升",
x = "训练样本数量", y = "交叉验证准确率"
) +
ggplot2::theme_minimal() +
ggplot2::ylim(0.5, 1)
# 模拟不同情况的学习曲线
train_sizes <- 1:100
# 欠拟合:性能低且平稳
underfit_curve <- 0.6 + 0.05 * (1 - exp(-train_sizes / 20))
# 过拟合:训练性能高,验证性能低
overfit_train <- 0.95 + 0.04 * (1 - exp(-train_sizes / 30))
overfit_val <- 0.5 + 0.3 * (1 - exp(-train_sizes / 50))
# 适度拟合:两条曲线接近
good_fit_train <- 0.85 + 0.1 * (1 - exp(-train_sizes / 20))
good_fit_val <- 0.8 + 0.12 * (1 - exp(-train_sizes / 25))
# 创建数据框
lc_data <- data.frame(
train_size = rep(train_sizes, 6),
accuracy = c(underfit_curve, underfit_curve,
overfit_train, overfit_val,
good_fit_train, good_fit_val),
type = rep(c("欠拟合", "过拟合", "适度拟合"), each = 200),
curve = rep(rep(c("训练集", "验证集"), each = 100), 3)
)
lc_data$type <- factor(lc_data$type, levels = c("欠拟合", "适度拟合", "过拟合"))
# 绘制
ggplot2::ggplot(lc_data, ggplot2::aes(x = train_size, y = accuracy, color = curve)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::facet_wrap(~ type) +
ggplot2::labs(
title = "不同情况下的学习曲线",
subtitle = "训练集曲线(红)vs 验证集曲线(蓝)",
x = "训练样本数量", y = "准确率", color = ""
) +
ggplot2::scale_color_manual(values = c("训练集" = "red", "验证集" = "blue")) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
方法:将数据分成K份,进行K次训练和验证,每次用其中1份验证、其余K-1份训练。
数学表达:
\[CV_{(K)} = \frac{1}{K} \sum_{i=1}^{K} E_i\]
其中 \(E_i\) 是第 \(i\) 折的验证误差。
# K折交叉验证示意图
# 使用5折交叉验证演示
set.seed(42)
n <- 50
data <- data.frame(
id = 1:n,
value = rnorm(n)
)
# 创建5折划分
k <- 5
folds <- sample(rep(1:k, length.out = n))
# 可视化交叉验证过程
cv_viz <- data.frame(
id = rep(1:n, k),
fold = rep(1:k, each = n),
role = rep(folds, k)
)
for (i in 1:k) {
cv_viz$role[(i-1)*n + 1:n] <- ifelse(folds == i, "验证集", "训练集")
}
cv_viz$fold <- factor(cv_viz$fold)
ggplot2::ggplot(cv_viz, ggplot2::aes(x = id, y = fold, fill = role)) +
ggplot2::geom_tile(color = "white", linewidth = 0.5) +
ggplot2::scale_fill_manual(values = c("训练集" = "lightblue", "验证集" = "salmon")) +
ggplot2::labs(
title = "5折交叉验证示意图",
subtitle = "每一行代表一次验证,红色为验证集,蓝色为训练集",
x = "样本编号", y = "折数", fill = "角色"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
# 使用caret包进行K折交叉验证
set.seed(42)
# 定义交叉验证控制
cv_control <- caret::trainControl(
method = "cv", # K折交叉验证
number = 5, # 5折
verboseIter = FALSE # 不显示每次迭代的进度
)
# 训练模型(使用5折交叉验证)
model_cv <- caret::train(
Species ~ .,
data = iris,
method = "rpart",
trControl = cv_control,
tuneLength = 3
)
# 查看交叉验证结果
print(model_cv)
## CART
##
## 150 samples
## 4 predictor
## 3 classes: 'setosa', 'versicolor', 'virginica'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 120, 120, 120, 120, 120
## Resampling results across tuning parameters:
##
## cp Accuracy Kappa
## 0.00 0.9400000 0.91
## 0.44 0.7600000 0.64
## 0.50 0.3333333 0.00
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.
方法:每次留一个样本做验证,其余全部做训练。适用于小样本数据。
特点: 1. 估计最准确(偏差最小) 2. 计算成本高(需要训练n次) 3. 方差可能较大
# LOOCV示例
# 使用caret包
set.seed(42)
loocv_control <- caret::trainControl(
method = "LOOCV" # 留一交叉验证
)
# 训练模型(使用LOOCV)
model_loocv <- caret::train(
Species ~ .,
data = iris,
method = "rpart",
trControl = loocv_control
)
# 查看结果
cat("LOOCV准确率:", model_loocv$results$Accuracy, "\n")
## LOOCV准确率: 0.9333333 0.3333333 0
方法:确保每折中各类别的比例与整体数据一致。
适用场景:类别不平衡的分类问题。
# 分层K折交叉验证
set.seed(42)
# 创建不平衡数据集
imbalanced_iris <- iris %>%
dplyr::filter(!(Species == "setosa" & row_number() > 10)) # setosa只保留10个
# 普通K折 vs 分层K折
# 使用caret的createFolds函数
normal_folds <- caret::createFolds(imbalanced_iris$Species, k = 5)
stratified_folds <- caret::createFolds(imbalanced_iris$Species, k = 5, returnTrain = TRUE)
# 比较每折中各类别的比例
compare_folds <- data.frame()
for (i in 1:5) {
test_idx <- normal_folds[[i]]
compare_folds <- rbind(compare_folds, data.frame(
折数 = i,
方法 = "普通K折",
setosa比例 = mean(imbalanced_iris$Species[test_idx] == "setosa")
))
}
# 分层K折的比例应该更一致
for (i in 1:5) {
test_idx <- normal_folds[[i]]
cat("折", i, "setosa比例:", mean(imbalanced_iris$Species[test_idx] == "setosa"), "\n")
}
## 折 1 setosa比例: 0.09090909
## 折 2 setosa比例: 0.09090909
## 折 3 setosa比例: 0.09090909
## 折 4 setosa比例: 0.09090909
## 折 5 setosa比例: 0.09090909
不同任务、不同场景需要关注不同的性能维度:
| 场景 | 关注重点 | 适合指标 |
|---|---|---|
| 疾病筛查 | 不漏诊 | 召回率 |
| 确诊检查 | 不误诊 | 精确率 |
| 综合评估 | 平衡 | F1分数 |
| 模型比较 | 整体性能 | AUC |
| 指标 | 公式 | 含义 |
|---|---|---|
| 准确率 | \(\frac{TP+TN}{TP+TN+FP+FN}\) | 预测正确的比例 |
| 精确率 | \(\frac{TP}{TP+FP}\) | 预测为正的样本中真正为正的比例 |
| 召回率 | \(\frac{TP}{TP+FN}\) | 真正为正的样本被正确预测的比例 |
| F1分数 | \(2 \times \frac{精确率 \times 召回率}{精确率 + 召回率}\) | 精确率和召回率的调和平均 |
| 指标 | 公式 | 含义 |
|---|---|---|
| MAE | \(\frac{1}{n}\sum|y_i - \hat{y}_i|\) | 平均绝对误差 |
| MSE | \(\frac{1}{n}\sum(y_i - \hat{y}_i)^2\) | 均方误差 |
| RMSE | \(\sqrt{MSE}\) | 均方根误差 |
| R² | \(1 - \frac{SS_{res}}{SS_{tot}}\) | 决定系数 |
# 评价指标计算示例
# 使用模拟数据演示
set.seed(42)
n <- 100
# 模拟真实值和预测值
actual <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
predicted_prob <- runif(n)
predicted <- ifelse(predicted_prob > 0.5, 1, 0)
# 计算混淆矩阵
conf_matrix <- table(实际 = actual, 预测 = predicted)
print(conf_matrix)
## 预测
## 实际 0 1
## 0 31 35
## 1 13 21
# 手动计算评价指标
TP <- conf_matrix[2, 2]
TN <- conf_matrix[1, 1]
FP <- conf_matrix[1, 2]
FN <- conf_matrix[2, 1]
accuracy <- (TP + TN) / sum(conf_matrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
f1 <- 2 * precision * recall / (precision + recall)
cat("准确率:", accuracy, "\n")
## 准确率: 0.52
cat("精确率:", precision, "\n")
## 精确率: 0.375
cat("召回率:", recall, "\n")
## 召回率: 0.6176471
cat("F1分数:", f1, "\n")
## F1分数: 0.4666667
# 使用caret包计算评价指标
set.seed(42)
# 训练模型
model <- caret::train(
Species ~ .,
data = iris,
method = "rpart",
trControl = caret::trainControl(method = "cv", number = 5)
)
# 查看详细结果
# 确保因子水平一致
obs <- factor(model$pred$obs, levels = levels(iris$Species))
pred <- factor(model$pred$pred, levels = levels(iris$Species))
conf_mat <- caret::confusionMatrix(pred, obs)
print(conf_mat)
## Confusion Matrix and Statistics
##
## Reference
## Prediction setosa versicolor virginica
## setosa 0 0 0
## versicolor 0 0 0
## virginica 0 0 0
##
## Overall Statistics
##
## Accuracy : NaN
## 95% CI : (NA, NA)
## No Information Rate : NA
## P-Value [Acc > NIR] : NA
##
## Kappa : NaN
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: setosa Class: versicolor Class: virginica
## Sensitivity NA NA NA
## Specificity NA NA NA
## Pos Pred Value NA NA NA
## Neg Pred Value NA NA NA
## Prevalence NaN NaN NaN
## Detection Rate NaN NaN NaN
## Detection Prevalence NaN NaN NaN
## Balanced Accuracy NA NA NA
核心意义:选择正确的评价指标是模型评估的基础,不同指标反映模型的不同侧面。
举例:癌症筛查模型的评价指标选择
场景:开发一个乳腺癌筛查模型,用于大规模人群筛查。
为什么召回率比精确率更重要?
在筛查场景中,漏诊(假阴性)的代价远高于误诊(假阳性)。漏诊可能导致患者错过最佳治疗时机,而误诊只是需要进一步检查确认。因此,应该优先优化召回率,确保尽可能多地识别出真正的患者。
如果精确率和召回率都很重要呢?
使用F1分数作为综合指标,或者根据实际成本设定不同的权重(如加权F分数)。
需要注意什么?
在类别不平衡的数据中,准确率可能具有误导性。例如,如果95%的样本是健康人,一个”全预测为健康”的模型准确率也有95%,但完全没有实用价值。
不同特征的量纲和数值范围可能差异巨大,这会导致:
| 方法 | 公式 | 结果范围 | 适用场景 |
|---|---|---|---|
| 标准化(Z-score) | \(\frac{x - \mu}{\sigma}\) | 均值0,标准差1 | 正态分布数据 |
| 归一化(Min-Max) | \(\frac{x - x_{min}}{x_{max} - x_{min}}\) | [0, 1] | 有明确边界的数据 |
| 鲁棒缩放 | \(\frac{x - Q_2}{Q_3 - Q_1}\) | 无固定范围 | 有异常值的数据 |
# 特征缩放示例
set.seed(42)
# 创建模拟医疗数据
medical_data <- data.frame(
age = round(rnorm(100, mean = 50, sd = 15)), # 年龄:20-80岁
blood_pressure = round(rnorm(100, mean = 120, sd = 20)), # 血压:80-180 mmHg
cholesterol = round(rnorm(100, mean = 200, sd = 50)), # 胆固醇:100-400 mg/dL
income = round(rnorm(100, mean = 50000, sd = 20000)) # 收入:1万-10万元
)
# 查看原始数据范围
summary(medical_data)
## age blood_pressure cholesterol income
## Min. : 5.00 Min. : 80.0 Min. : 65.0 Min. :16350
## 1st Qu.:40.75 1st Qu.:108.0 1st Qu.:164.8 1st Qu.:39346
## Median :51.00 Median :118.5 Median :199.0 Median :49086
## Mean :50.51 Mean :118.2 Mean :199.5 Mean :50659
## 3rd Qu.:60.00 3rd Qu.:129.2 3rd Qu.:232.8 3rd Qu.:63496
## Max. :84.00 Max. :174.0 Max. :323.0 Max. :98443
公式:
\[x_{scaled} = \frac{x - \mu}{\sigma}\]
特点: 1. 转换后均值为0,标准差为1 2. 保留异常值信息 3. 适用于近似正态分布的数据
# 标准化实现
scale_standard <- function(x) {
(x - mean(x)) / sd(x)
}
# 应用标准化
medical_data_std <- as.data.frame(lapply(medical_data, scale_standard))
# 查看标准化后的数据
summary(medical_data_std)
## age blood_pressure cholesterol income
## Min. :-2.90227 Min. :-2.11674 Min. :-2.643720 Min. :-1.95782
## 1st Qu.:-0.62242 1st Qu.:-0.56723 1st Qu.:-0.682459 1st Qu.:-0.64557
## Median : 0.03125 Median : 0.01383 Median :-0.009044 Median :-0.08975
## Mean : 0.00000 Mean : 0.00000 Mean : 0.000000 Mean : 0.00000
## 3rd Qu.: 0.60520 3rd Qu.: 0.60874 3rd Qu.: 0.654540 3rd Qu.: 0.73253
## Max. : 2.13573 Max. : 3.08519 Max. : 2.429013 Max. : 2.72680
# 使用caret包的preProcess函数
preproc_std <- caret::preProcess(medical_data, method = c("center", "scale"))
medical_data_std2 <- predict(preproc_std, medical_data)
# 可视化对比
viz_data <- data.frame(
特征 = rep(names(medical_data), each = nrow(medical_data)),
原始值 = unlist(medical_data),
标准化后 = unlist(medical_data_std)
)
viz_data$特征 <- factor(viz_data$特征, levels = names(medical_data))
p1 <- ggplot2::ggplot(viz_data, ggplot2::aes(x = 特征, y = 原始值, fill = 特征)) +
ggplot2::geom_boxplot() +
ggplot2::labs(title = "原始数据分布", y = "数值") +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "none", axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))
p2 <- ggplot2::ggplot(viz_data, ggplot2::aes(x = 特征, y = 标准化后, fill = 特征)) +
ggplot2::geom_boxplot() +
ggplot2::labs(title = "标准化后数据分布", y = "标准化数值") +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "none", axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))
p1 + p2 + patchwork::plot_annotation(
title = "特征缩放前后对比",
subtitle = "标准化消除了量纲差异"
)
公式:
\[x_{scaled} = \frac{x - x_{min}}{x_{max} - x_{min}}\]
特点: 1. 转换后数据在[0, 1]范围内 2. 对异常值敏感 3. 适用于有明确边界的数据
# 归一化实现
scale_minmax <- function(x) {
(x - min(x)) / (max(x) - min(x))
}
# 应用归一化
medical_data_norm <- as.data.frame(lapply(medical_data, scale_minmax))
# 查看归一化后的数据范围
summary(medical_data_norm)
## age blood_pressure cholesterol income
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.4525 1st Qu.:0.2979 1st Qu.:0.3866 1st Qu.:0.2801
## Median :0.5823 Median :0.4096 Median :0.5194 Median :0.3988
## Mean :0.5761 Mean :0.4069 Mean :0.5212 Mean :0.4179
## 3rd Qu.:0.6962 3rd Qu.:0.5239 3rd Qu.:0.6502 3rd Qu.:0.5743
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
# 使用caret包
preproc_norm <- caret::preProcess(medical_data, method = "range")
medical_data_norm2 <- predict(preproc_norm, medical_data)
公式:
\[x_{scaled} = \frac{x - Q_2}{IQR} = \frac{x - Q_2}{Q_3 - Q_1}\]
特点: 1. 使用中位数和四分位距,对异常值不敏感 2. 适用于有异常值的数据 3. 结果没有固定范围
# 鲁棒缩放实现
scale_robust <- function(x) {
(x - median(x)) / IQR(x)
}
# 添加异常值演示
medical_data_outlier <- medical_data
medical_data_outlier$cholesterol[1] <- 1000 # 添加异常值
# 对比标准化和鲁棒缩放对异常值的处理
std_with_outlier <- scale_standard(medical_data_outlier$cholesterol)
robust_with_outlier <- scale_robust(medical_data_outlier$cholesterol)
cat("标准化 - 均值:", mean(std_with_outlier), "标准差:", sd(std_with_outlier), "\n")
## 标准化 - 均值: -8.632418e-17 标准差: 1
cat("鲁棒缩放 - 均值:", mean(robust_with_outlier), "标准差:", sd(robust_with_outlier), "\n")
## 鲁棒缩放 - 均值: 0.1125088 标准差: 1.331805
核心意义:特征缩放是数据预处理的基础步骤,确保不同特征对模型的贡献均衡。
举例:基于患者生理指标预测疾病风险
为什么需要特征缩放?
假设模型使用年龄(20-80岁)、血压(80-180 mmHg)、血糖(3-30 mmol/L)三个特征。如果不缩放,血压的数值范围最大,在KNN等距离敏感算法中会主导距离计算,导致模型主要根据血压做决策,忽略其他特征。
如何选择缩放方法?
需要注意什么?
缩放参数(均值、标准差、最小值、最大值)必须只从训练集计算,然后应用到测试集。不能在整个数据集上计算缩放参数,否则会造成数据泄露。
| 方法 | 原理 | 适用场景 | 缺点 |
|---|---|---|---|
| 独热编码 | 每个类别创建一个二进制列 | 名义变量、类别少 | 维度爆炸 |
| 标签编码 | 每个类别映射为一个整数 | 有序变量 | 引入虚假顺序 |
| 目标编码 | 用目标变量的均值编码 | 高基数变量、树模型 | 容易过拟合 |
# 创建模拟医疗数据
set.seed(42)
patient_data <- data.frame(
patient_id = 1:20,
gender = sample(c("Male", "Female"), 20, replace = TRUE),
blood_type = sample(c("A", "B", "AB", "O"), 20, replace = TRUE),
disease_stage = sample(c("I", "II", "III", "IV"), 20, replace = TRUE),
hospital = sample(paste0("H", 1:10), 20, replace = TRUE), # 高基数变量
outcome = sample(c(0, 1), 20, replace = TRUE)
)
head(patient_data)
## patient_id gender blood_type disease_stage hospital outcome
## 1 1 Male AB III H5 1
## 2 2 Male AB III H10 0
## 3 3 Male A II H2 1
## 4 4 Male A II H6 0
## 5 5 Female AB IV H6 0
## 6 6 Female O IV H2 1
原理:为每个类别创建一个二进制列,存在为1,不存在为0。
# 独热编码实现
# 使用caret包的dummyVars函数
# 对性别进行独热编码
dummy_model <- caret::dummyVars(~ gender, data = patient_data)
gender_encoded <- predict(dummy_model, patient_data)
head(gender_encoded)
## genderFemale genderMale
## 1 0 1
## 2 0 1
## 3 0 1
## 4 0 1
## 5 1 0
## 6 1 0
# 对血型进行独热编码
dummy_model2 <- caret::dummyVars(~ blood_type, data = patient_data)
blood_type_encoded <- predict(dummy_model2, patient_data)
head(blood_type_encoded)
## blood_typeA blood_typeAB blood_typeB blood_typeO
## 1 0 1 0 0
## 2 0 1 0 0
## 3 1 0 0 0
## 4 1 0 0 0
## 5 0 1 0 0
## 6 0 0 0 1
# 使用model.matrix函数(R基础方法)
gender_onehot <- model.matrix(~ gender - 1, data = patient_data)
head(gender_onehot)
## genderFemale genderMale
## 1 0 1
## 2 0 1
## 3 0 1
## 4 0 1
## 5 1 0
## 6 1 0
原理:将每个类别映射为一个整数。
适用场景:有序分类变量(如疾病分期I < II < III < IV)。
# 标签编码实现
# 使用因子转换
# 疾病分期是有序变量
patient_data$disease_stage_ordered <- factor(
patient_data$disease_stage,
levels = c("I", "II", "III", "IV"),
ordered = TRUE
)
# 查看编码
patient_data$disease_stage_encoded <- as.numeric(patient_data$disease_stage_ordered)
head(patient_data[, c("disease_stage", "disease_stage_ordered", "disease_stage_encoded")])
## disease_stage disease_stage_ordered disease_stage_encoded
## 1 III III 3
## 2 III III 3
## 3 II II 2
## 4 II II 2
## 5 IV IV 4
## 6 IV IV 4
# 注意:对于名义变量,标签编码会引入虚假顺序
# 例如血型的标签编码
patient_data$blood_type_encoded <- as.numeric(factor(patient_data$blood_type))
# 这会引入虚假的顺序关系:A < B < AB < O,但实际上血型没有顺序
原理:用该类别对应的目标变量均值来编码。
公式:
\[x_{encoded} = \frac{\sum_{i \in C} y_i + \alpha \cdot \bar{y}}{n_C + \alpha}\]
其中 \(C\) 是当前类别,\(n_C\) 是该类别的样本数,\(\bar{y}\) 是全局目标均值,\(\alpha\) 是平滑参数。
# 目标编码实现
# 对医院变量进行目标编码
# 计算每个医院的目标变量均值
hospital_target_mean <- patient_data %>%
dplyr::group_by(hospital) %>%
dplyr::summarise(
target_mean = mean(outcome),
count = n(),
.groups = "drop"
)
# 全局目标均值
global_mean <- mean(patient_data$outcome)
# 平滑参数
alpha <- 5
# 计算平滑后的目标编码
hospital_target_mean <- hospital_target_mean %>%
dplyr::mutate(
target_encoded = (target_mean * count + global_mean * alpha) / (count + alpha)
)
# 合并到原数据
patient_data <- patient_data %>%
dplyr::left_join(
hospital_target_mean %>% dplyr::select(hospital, target_encoded),
by = "hospital"
)
print(hospital_target_mean)
## # A tibble: 9 × 4
## hospital target_mean count target_encoded
## <chr> <dbl> <int> <dbl>
## 1 H1 1 2 0.679
## 2 H10 0.667 3 0.594
## 3 H2 0.667 3 0.594
## 4 H3 0 1 0.458
## 5 H4 1 1 0.625
## 6 H5 0.333 3 0.469
## 7 H6 0 4 0.306
## 8 H7 1 2 0.679
## 9 H8 1 1 0.625
| 类型 | 定义 | 示例 | 处理策略 |
|---|---|---|---|
| MCAR | 完全随机缺失 | 数据录入随机遗漏 | 可删除或填补 |
| MAR | 随机缺失 | 男性更可能不报告体重 | 可用其他变量预测 |
| MNAR | 非随机缺失 | 重症患者未完成问卷 | 需要专门方法 |
# 创建带缺失值的数据
set.seed(42)
data_with_na <- data.frame(
patient_id = 1:100,
age = round(rnorm(100, 50, 15)),
bmi = round(rnorm(100, 25, 4), 1),
blood_pressure = round(rnorm(100, 120, 20)),
cholesterol = round(rnorm(100, 200, 50))
)
# 随机引入缺失值
data_with_na$bmi[sample(1:100, 15)] <- NA
data_with_na$blood_pressure[sample(1:100, 10)] <- NA
data_with_na$cholesterol[sample(1:100, 20)] <- NA
# 查看缺失情况
missing_summary <- data.frame(
变量 = names(data_with_na)[-1],
缺失数 = sapply(data_with_na[-1], function(x) sum(is.na(x))),
缺失比例 = sapply(data_with_na[-1], function(x) mean(is.na(x)) * 100)
)
print(missing_summary)
## 变量 缺失数 缺失比例
## age age 0 0
## bmi bmi 15 15
## blood_pressure blood_pressure 10 10
## cholesterol cholesterol 20 20
# 可视化缺失模式
visdat <- data_with_na[-1]
visdat$rown <- 1:nrow(visdat)
visdat_long <- tidyr::pivot_longer(visdat, cols = -rown, names_to = "variable", values_to = "value")
visdat_long$is_na <- is.na(visdat_long$value)
ggplot2::ggplot(visdat_long, ggplot2::aes(x = variable, y = rown, fill = is_na)) +
ggplot2::geom_tile() +
ggplot2::scale_fill_manual(values = c("FALSE" = "lightblue", "TRUE" = "red")) +
ggplot2::labs(
title = "缺失值模式可视化",
x = "变量", y = "样本", fill = "是否缺失"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
# 均值填补
data_mean_imputed <- data_with_na
data_mean_imputed$bmi[is.na(data_mean_imputed$bmi)] <- mean(data_mean_imputed$bmi, na.rm = TRUE)
data_mean_imputed$blood_pressure[is.na(data_mean_imputed$blood_pressure)] <- mean(data_mean_imputed$blood_pressure, na.rm = TRUE)
data_mean_imputed$cholesterol[is.na(data_mean_imputed$cholesterol)] <- mean(data_mean_imputed$cholesterol, na.rm = TRUE)
# 中位数填补(对异常值更稳健)
data_median_imputed <- data_with_na
data_median_imputed$bmi[is.na(data_median_imputed$bmi)] <- median(data_median_imputed$bmi, na.rm = TRUE)
data_median_imputed$blood_pressure[is.na(data_median_imputed$blood_pressure)] <- median(data_median_imputed$blood_pressure, na.rm = TRUE)
data_median_imputed$cholesterol[is.na(data_median_imputed$cholesterol)] <- median(data_median_imputed$cholesterol, na.rm = TRUE)
# 使用caret包的preProcess
preproc_impute <- caret::preProcess(data_with_na[-1], method = "medianImpute")
data_caret_imputed <- predict(preproc_impute, data_with_na[-1])
# 比较填补结果
impute_compare <- data.frame(
变量 = c("bmi", "blood_pressure", "cholesterol"),
原始均值 = c(mean(data_with_na$bmi, na.rm = TRUE),
mean(data_with_na$blood_pressure, na.rm = TRUE),
mean(data_with_na$cholesterol, na.rm = TRUE)),
均值填补后均值 = c(mean(data_mean_imputed$bmi),
mean(data_mean_imputed$blood_pressure),
mean(data_mean_imputed$cholesterol)),
中位数填补后均值 = c(mean(data_median_imputed$bmi),
mean(data_median_imputed$blood_pressure),
mean(data_median_imputed$cholesterol))
)
print(impute_compare)
## 变量 原始均值 均值填补后均值 中位数填补后均值
## 1 bmi 24.66353 24.66353 24.669
## 2 blood_pressure 121.02222 121.02222 121.020
## 3 cholesterol 202.73750 202.73750 201.690
原理:利用K个最相似样本的值进行填补。
# KNN填补
preproc_knn <- caret::preProcess(data_with_na[-1], method = "knnImpute")
data_knn_imputed <- predict(preproc_knn, data_with_na[-1])
# 注意:KNN填补需要先标准化数据
# caret会自动处理
# 查看填补结果
head(data_knn_imputed)
## age bmi blood_pressure cholesterol
## 1 1.3066901 1.4338767 -2.0204150 -0.06333831
## 2 -0.5427005 1.2663831 0.2944158 0.81587837
## 3 0.2863367 -0.8440364 1.0824434 -0.01706374
## 4 0.5414250 2.1596823 1.9689744 -0.46129954
## 5 0.3501088 -0.6597935 -1.4293944 -0.22529927
## 6 -0.1600679 0.2055902 -1.1831358 -0.13275015
原理:用其他变量作为特征,训练模型预测缺失值。
# 使用mice包进行多重插补(如果可用)
# 这里演示简单的模型预测填补
# 用其他变量预测bmi的缺失值
bmi_complete <- data_with_na[!is.na(data_with_na$bmi), ]
bmi_missing <- data_with_na[is.na(data_with_na$bmi), ]
# 训练线性回归模型
bmi_model <- lm(bmi ~ age + blood_pressure + cholesterol, data = bmi_complete)
# 预测缺失值
predicted_bmi <- predict(bmi_model, newdata = bmi_missing)
# 填补
data_model_imputed <- data_with_na
data_model_imputed$bmi[is.na(data_model_imputed$bmi)] <- predicted_bmi
# 可视化填补效果
ggplot2::ggplot() +
ggplot2::geom_histogram(ggplot2::aes(x = data_with_na$bmi, fill = "原始(含缺失)"),
alpha = 0.5, bins = 20, na.rm = TRUE) +
ggplot2::geom_histogram(ggplot2::aes(x = data_model_imputed$bmi, fill = "模型填补后"),
alpha = 0.5, bins = 20) +
ggplot2::labs(
title = "模型预测填补效果",
x = "BMI", y = "频数", fill = ""
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
| 方法 | 原理 | 优点 | 缺点 |
|---|---|---|---|
| 过滤法 | 基于统计指标筛选 | 快速、独立于模型 | 忽略特征间关系 |
| 包装法 | 用模型性能评估特征组合 | 考虑特征间关系 | 计算成本高 |
| 嵌入法 | 在模型训练中自动选择 | 高效、模型特定 | 依赖特定模型 |
# 使用Sonar数据集演示特征选择
data(Sonar, package = "mlbench")
# 查看数据维度
cat("数据维度:", dim(Sonar), "\n")
## 数据维度: 208 61
cat("特征数量:", ncol(Sonar) - 1, "\n")
## 特征数量: 60
# 划分训练集和测试集
set.seed(42)
train_index <- caret::createDataPartition(Sonar$Class, p = 0.7, list = FALSE)
train_data <- Sonar[train_index, ]
test_data <- Sonar[-train_index, ]
原理:基于统计指标(如相关性、方差、信息增益)评估每个特征的重要性。
# 方差阈值法:删除方差过小的特征
feature_vars <- apply(train_data[, -ncol(train_data)], 2, var)
low_var_features <- names(feature_vars[feature_vars < 0.01])
cat("低方差特征数量:", length(low_var_features), "\n")
## 低方差特征数量: 22
# 相关性过滤:删除与目标变量相关性低的特征
# 计算每个特征与目标变量的相关性
correlations <- sapply(names(train_data)[-ncol(train_data)], function(col) {
cor(as.numeric(train_data[[col]]), as.numeric(train_data$Class))
})
# 选择相关性最高的前20个特征
top_features <- names(sort(abs(correlations), decreasing = TRUE))[1:20]
# 可视化特征相关性
cor_df <- data.frame(
feature = names(correlations),
correlation = abs(correlations)
)
cor_df <- cor_df[order(cor_df$correlation, decreasing = TRUE), ]
cor_df$feature <- factor(cor_df$feature, levels = cor_df$feature)
ggplot2::ggplot(head(cor_df, 30), ggplot2::aes(x = feature, y = correlation)) +
ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
ggplot2::labs(
title = "特征与目标变量的相关性(前30)",
x = "特征", y = "相关系数绝对值"
) +
ggplot2::theme_minimal() +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))
原理:通过训练模型评估特征子集的性能,搜索最优特征组合。
# 递归特征消除(Recursive Feature Elimination, RFE)
set.seed(42)
# 定义控制参数
rfe_control <- caret::rfeControl(
functions = caret::rfFuncs, # 使用随机森林
method = "cv", # 交叉验证
number = 5 # 5折
)
# 执行RFE
rfe_result <- caret::rfe(
x = train_data[, -ncol(train_data)],
y = train_data$Class,
sizes = c(5, 10, 20, 30, 40, 50),
rfeControl = rfe_control
)
# 查看结果
print(rfe_result)
##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (5 fold)
##
## Resampling performance over subset size:
##
## Variables Accuracy Kappa AccuracySD KappaSD Selected
## 5 0.7875 0.5714 0.08620 0.17187
## 10 0.7870 0.5691 0.07636 0.15293
## 20 0.8208 0.6385 0.05967 0.12223
## 30 0.8489 0.6955 0.04050 0.08250
## 40 0.8553 0.7081 0.05310 0.10718
## 50 0.8351 0.6672 0.03162 0.06428
## 60 0.8694 0.7364 0.05197 0.10351 *
##
## The top 5 variables (out of 60):
## V11, V12, V9, V52, V36
# 最优特征数量
cat("\n最优特征数量:", rfe_result$bestSubset, "\n")
##
## 最优特征数量: 60
cat("选中的特征:\n")
## 选中的特征:
print(rfe_result$optVariables)
## [1] "V11" "V12" "V9" "V52" "V36" "V48" "V10" "V13" "V47" "V28" "V46" "V21"
## [13] "V49" "V27" "V37" "V51" "V20" "V16" "V17" "V18" "V2" "V45" "V14" "V5"
## [25] "V1" "V4" "V34" "V23" "V8" "V31" "V15" "V3" "V33" "V19" "V44" "V59"
## [37] "V42" "V6" "V43" "V24" "V26" "V25" "V54" "V35" "V53" "V39" "V30" "V32"
## [49] "V60" "V58" "V22" "V40" "V38" "V55" "V7" "V56" "V29" "V41" "V57" "V50"
原理:在模型训练过程中自动进行特征选择。
# 使用带正则化的模型(如LASSO)
# 使用glmnet包
# 准备数据
x_train <- as.matrix(train_data[, -ncol(train_data)])
y_train <- ifelse(train_data$Class == "M", 1, 0)
# 使用caret训练LASSO模型
set.seed(42)
lasso_model <- caret::train(
x = x_train,
y = train_data$Class,
method = "glmnet",
trControl = caret::trainControl(method = "cv", number = 5),
tuneGrid = data.frame(alpha = 1, lambda = seq(0.001, 0.1, length.out = 20))
)
# 查看模型
print(lasso_model)
## glmnet
##
## 146 samples
## 60 predictor
## 2 classes: 'M', 'R'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 117, 118, 116, 117, 116
## Resampling results across tuning parameters:
##
## lambda Accuracy Kappa
## 0.001000000 0.7386864 0.4749746
## 0.006210526 0.7796223 0.5568986
## 0.011421053 0.7660755 0.5296068
## 0.016631579 0.7724959 0.5426466
## 0.021842105 0.7582102 0.5130910
## 0.027052632 0.7786700 0.5554115
## 0.032263158 0.7786700 0.5554115
## 0.037473684 0.7791461 0.5550158
## 0.042684211 0.7860427 0.5685974
## 0.047894737 0.7998522 0.5952532
## 0.053105263 0.8003448 0.5965009
## 0.058315789 0.8005747 0.5962016
## 0.063526316 0.8008046 0.5973741
## 0.068736842 0.7939080 0.5825084
## 0.073947368 0.7867652 0.5688028
## 0.079157895 0.7870115 0.5695442
## 0.084368421 0.7870115 0.5695442
## 0.089578947 0.7732020 0.5410120
## 0.094789474 0.7798686 0.5549569
## 0.100000000 0.7658292 0.5260499
##
## Tuning parameter 'alpha' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were alpha = 1 and lambda = 0.06352632.
# 获取非零系数的特征
lasso_coefs <- coef(lasso_model$finalModel, lasso_model$bestTune$lambda)
selected_features <- rownames(lasso_coefs)[lasso_coefs[, 1] != 0][-1] # 排除截距
cat("\nLASSO选中的特征数量:", length(selected_features), "\n")
##
## LASSO选中的特征数量: 9
cat("选中的特征:\n")
## 选中的特征:
print(selected_features)
## [1] "V11" "V12" "V21" "V28" "V36" "V45" "V48" "V49" "V52"
原理:找到数据方差最大的方向,将数据投影到低维空间。
数学表达:
寻找投影方向 \(w\),使得投影后的方差最大:
\[\max_w \text{Var}(Xw) = \max_w w^T \Sigma w\]
其中 \(\Sigma\) 是协方差矩阵。
# PCA示例:使用iris数据集
data(iris)
# 只取数值特征
iris_numeric <- iris[, 1:4]
# 执行PCA
pca_result <- prcomp(iris_numeric, scale. = TRUE)
# 查看主成分
print(summary(pca_result))
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion 0.7296 0.9581 0.99482 1.00000
# 方差解释比例
var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2)
var_df <- data.frame(
PC = paste0("PC", 1:4),
方差解释比例 = var_explained,
累积解释比例 = cumsum(var_explained)
)
print(var_df)
## PC 方差解释比例 累积解释比例
## 1 PC1 0.729624454 0.7296245
## 2 PC2 0.228507618 0.9581321
## 3 PC3 0.036689219 0.9948213
## 4 PC4 0.005178709 1.0000000
# 碎石图
ggplot2::ggplot(var_df, ggplot2::aes(x = PC, y = 方差解释比例)) +
ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
ggplot2::geom_line(ggplot2::aes(group = 1), color = "red", linewidth = 1) +
ggplot2::geom_point(color = "red", size = 2) +
ggplot2::labs(
title = "PCA碎石图",
x = "主成分", y = "方差解释比例"
) +
ggplot2::theme_minimal()
# 可视化前两个主成分
pca_scores <- as.data.frame(pca_result$x)
pca_scores$Species <- iris$Species
ggplot2::ggplot(pca_scores, ggplot2::aes(x = PC1, y = PC2, color = Species)) +
ggplot2::geom_point(size = 3, alpha = 0.7) +
ggplot2::labs(
title = "PCA降维可视化",
subtitle = "前两个主成分解释了95.8%的方差",
x = paste0("PC1 (", round(var_explained[1] * 100, 1), "%)"),
y = paste0("PC2 (", round(var_explained[2] * 100, 1), "%)")
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
原理:找到能最好区分类别的投影方向,最大化类间方差、最小化类内方差。
数学表达:
\[\max_w \frac{w^T S_B w}{w^T S_W w}\]
其中 \(S_B\) 是类间散度矩阵,\(S_W\) 是类内散度矩阵。
# LDA示例
# 使用MASS包的lda函数
lda_result <- MASS::lda(Species ~ ., data = iris)
# 查看结果
print(lda_result)
## Call:
## lda(Species ~ ., data = iris)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.006 3.428 1.462 0.246
## versicolor 5.936 2.770 4.260 1.326
## virginica 6.588 2.974 5.552 2.026
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.8293776 -0.02410215
## Sepal.Width 1.5344731 -2.16452123
## Petal.Length -2.2012117 0.93192121
## Petal.Width -2.8104603 -2.83918785
##
## Proportion of trace:
## LD1 LD2
## 0.9912 0.0088
# 预测并可视化
lda_pred <- predict(lda_result, iris)
lda_scores <- as.data.frame(lda_pred$x)
lda_scores$Species <- iris$Species
ggplot2::ggplot(lda_scores, ggplot2::aes(x = LD1, y = LD2, color = Species)) +
ggplot2::geom_point(size = 3, alpha = 0.7) +
ggplot2::labs(
title = "LDA降维可视化",
subtitle = "LDA寻找最佳分类方向",
x = "LD1", y = "LD2"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
# 对比PCA和LDA的降维效果
p1 <- ggplot2::ggplot(pca_scores, ggplot2::aes(x = PC1, y = PC2, color = Species)) +
ggplot2::geom_point(size = 2, alpha = 0.7) +
ggplot2::labs(title = "PCA:无监督降维", x = "PC1", y = "PC2") +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "none")
p2 <- ggplot2::ggplot(lda_scores, ggplot2::aes(x = LD1, y = LD2, color = Species)) +
ggplot2::geom_point(size = 2, alpha = 0.7) +
ggplot2::labs(title = "LDA:有监督降维", x = "LD1", y = "LD2") +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "none")
p1 + p2 + patchwork::plot_annotation(
title = "PCA vs LDA 降维效果对比",
subtitle = "LDA利用类别信息,分离效果更好"
)
有时特征之间存在交互效应,单个特征的影响取决于另一个特征的值。
医学示例:药物A的疗效可能取决于患者是否同时服用药物B。
# 创建模拟数据演示特征交互
set.seed(42)
n <- 200
interaction_data <- data.frame(
age = round(runif(n, 30, 80)),
bmi = round(rnorm(n, 25, 4), 1),
smoking = sample(c(0, 1), n, replace = TRUE),
exercise = sample(c(0, 1), n, replace = TRUE)
)
# 生成目标变量(包含交互效应)
# 吸烟对风险的影响随年龄增加而增强
interaction_data$risk_score <-
0.5 * interaction_data$age +
2 * interaction_data$bmi +
20 * interaction_data$smoking +
0.3 * interaction_data$age * interaction_data$smoking + # 交互效应
rnorm(n, 0, 5)
# 查看数据
head(interaction_data)
## age bmi smoking exercise risk_score
## 1 76 29.8 1 1 145.54570
## 2 77 29.2 1 1 144.57387
## 3 44 21.0 1 1 97.18772
## 4 72 32.4 1 0 143.08005
## 5 62 22.3 1 1 110.59923
## 6 56 25.4 0 1 77.80938
# 手动创建交互特征
interaction_data$age_smoking <- interaction_data$age * interaction_data$smoking
interaction_data$bmi_exercise <- interaction_data$bmi * interaction_data$exercise
# 使用模型验证交互效应的重要性
model_no_interaction <- lm(risk_score ~ age + bmi + smoking + exercise,
data = interaction_data)
model_with_interaction <- lm(risk_score ~ age + bmi + smoking + exercise +
age:smoking + bmi:exercise,
data = interaction_data)
# 比较模型
anova_result <- anova(model_no_interaction, model_with_interaction)
print(anova_result)
## Analysis of Variance Table
##
## Model 1: risk_score ~ age + bmi + smoking + exercise
## Model 2: risk_score ~ age + bmi + smoking + exercise + age:smoking + bmi:exercise
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 195 5809.0
## 2 193 4988.8 2 820.14 15.864 4.181e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 查看交互项系数
summary(model_with_interaction)
##
## Call:
## lm(formula = risk_score ~ age + bmi + smoking + exercise + age:smoking +
## bmi:exercise, data = interaction_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.3919 -3.5348 -0.0625 3.7397 17.4802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.48097 3.80750 -0.652 0.515
## age 0.54573 0.03428 15.919 < 2e-16 ***
## bmi 1.98366 0.13167 15.065 < 2e-16 ***
## smoking 22.58550 2.88465 7.830 3.17e-13 ***
## exercise 5.33206 4.77227 1.117 0.265
## age:smoking 0.26817 0.04969 5.397 1.97e-07 ***
## bmi:exercise -0.23562 0.18994 -1.240 0.216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.084 on 193 degrees of freedom
## Multiple R-squared: 0.9562, Adjusted R-squared: 0.9548
## F-statistic: 702.1 on 6 and 193 DF, p-value: < 2.2e-16
原理:添加特征的幂次项,捕捉非线性关系。
# 多项式特征示例
set.seed(42)
n <- 100
# 生成非线性关系数据
poly_data <- data.frame(
x = runif(n, -5, 5)
)
poly_data$y <- 0.5 * poly_data$x^2 - 2 * poly_data$x + 3 + rnorm(n, 0, 2)
# 线性模型
model_linear <- lm(y ~ x, data = poly_data)
# 多项式模型(二次)
model_poly2 <- lm(y ~ poly(x, degree = 2, raw = TRUE), data = poly_data)
# 多项式模型(三次)
model_poly3 <- lm(y ~ poly(x, degree = 3, raw = TRUE), data = poly_data)
# 预测
poly_data$pred_linear <- predict(model_linear)
poly_data$pred_poly2 <- predict(model_poly2)
poly_data$pred_poly3 <- predict(model_poly3)
# 可视化
ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(alpha = 0.5) +
ggplot2::geom_line(ggplot2::aes(y = pred_linear, color = "线性"), linewidth = 1) +
ggplot2::geom_line(ggplot2::aes(y = pred_poly2, color = "二次多项式"), linewidth = 1) +
ggplot2::geom_line(ggplot2::aes(y = pred_poly3, color = "三次多项式"), linewidth = 1) +
ggplot2::labs(
title = "多项式特征捕捉非线性关系",
x = "x", y = "y", color = "模型"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
在医学数据中,类别不平衡非常常见,如罕见疾病的诊断。
问题:模型可能倾向于预测多数类,导致少数类的预测性能很差。
# 创建不平衡数据集
set.seed(42)
n <- 1000
imbalanced_data <- data.frame(
feature1 = c(rnorm(100, mean = 2, sd = 1), rnorm(900, mean = 0, sd = 1)),
feature2 = c(rnorm(100, mean = 2, sd = 1), rnorm(900, mean = 0, sd = 1)),
class = factor(c(rep("Positive", 100), rep("Negative", 900)))
)
# 查看类别分布
table(imbalanced_data$class)
##
## Negative Positive
## 900 100
# 可视化
ggplot2::ggplot(imbalanced_data, ggplot2::aes(x = feature1, y = feature2, color = class)) +
ggplot2::geom_point(size = 1.5, alpha = 0.6) +
ggplot2::labs(
title = "不平衡数据集可视化",
subtitle = "Positive: 100 (10%), Negative: 900 (90%)",
x = "Feature 1", y = "Feature 2"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
原理:增加少数类样本,使两类数量接近。
# 过采样:复制少数类样本
positive_samples <- imbalanced_data[imbalanced_data$class == "Positive", ]
negative_samples <- imbalanced_data[imbalanced_data$class == "Negative", ]
# 简单复制过采样
oversampled <- rbind(
positive_samples[sample(1:nrow(positive_samples), 800, replace = TRUE), ],
negative_samples
)
# 查看平衡后的分布
table(oversampled$class)
##
## Negative Positive
## 900 800
原理:减少多数类样本,使两类数量接近。
# 欠采样:随机抽取多数类样本
undersampled <- rbind(
positive_samples,
negative_samples[sample(1:nrow(negative_samples), 100), ]
)
# 查看平衡后的分布
table(undersampled$class)
##
## Negative Positive
## 100 100
原理:在少数类样本之间插值生成新样本。
# SMOTE原理示意
set.seed(42)
# 选择两个少数类样本
sample1 <- as.numeric(positive_samples[1, 1:2])
sample2 <- as.numeric(positive_samples[2, 1:2])
# 在两点之间随机插值生成新样本
alpha <- runif(1)
new_sample <- sample1 + alpha * (sample2 - sample1)
cat("原始样本1:", sample1, "\n")
## 原始样本1: 3.370958 4.325058
cat("原始样本2:", sample2, "\n")
## 原始样本2: 1.435302 2.524122
cat("生成的新样本:", new_sample, "\n")
## 生成的新样本: 1.600208 2.677551
# 可视化SMOTE
smote_viz <- data.frame(
x = c(sample1[1], sample2[1], new_sample[1]),
y = c(sample1[2], sample2[2], new_sample[2]),
type = c("原始样本1", "原始样本2", "SMOTE生成")
)
ggplot2::ggplot(smote_viz, ggplot2::aes(x = x, y = y, color = type)) +
ggplot2::geom_point(size = 3) +
ggplot2::geom_line(data = data.frame(x = c(sample1[1], sample2[1]),
y = c(sample1[2], sample2[2])),
ggplot2::aes(x = x, y = y), color = "gray50", linetype = "dashed") +
ggplot2::labs(
title = "SMOTE原理示意",
subtitle = "在少数类样本之间插值生成新样本",
x = "Feature 1", y = "Feature 2", color = ""
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
对于时间序列数据,需要创建能够捕捉时间依赖性的特征。
# 创建模拟时间序列数据
set.seed(42)
n_days <- 365
time_data <- data.frame(
date = seq(as.Date("2023-01-01"), by = "day", length.out = n_days),
patients = round(100 + 20 * sin(2 * pi * 1:n_days / 365) + # 年周期
5 * sin(2 * pi * 1:n_days / 7) + # 周周期
rnorm(n_days, 0, 10)) # 随机波动
)
# 查看数据
head(time_data)
## date patients
## 1 2023-01-01 118
## 2 2023-01-02 100
## 3 2023-01-03 107
## 4 2023-01-04 106
## 5 2023-01-05 101
## 6 2023-01-06 97
# 可视化
ggplot2::ggplot(time_data, ggplot2::aes(x = date, y = patients)) +
ggplot2::geom_line(color = "steelblue") +
ggplot2::labs(
title = "医院每日门诊量时间序列",
x = "日期", y = "门诊人数"
) +
ggplot2::theme_minimal()
原理:使用过去时间点的值作为当前时间点的特征。
# 创建滞后特征
time_data <- time_data %>%
dplyr::mutate(
lag_1 = dplyr::lag(patients, 1), # 前1天
lag_7 = dplyr::lag(patients, 7), # 前7天(上周同一天)
lag_30 = dplyr::lag(patients, 30) # 前30天
)
# 查看滞后特征
head(time_data[30:35, ], 10)
## date patients lag_1 lag_7 lag_30
## 30 2023-01-30 108 118 111 NA
## 31 2023-01-31 117 108 122 118
## 32 2023-02-01 115 117 125 100
## 33 2023-02-02 116 115 99 107
## 34 2023-02-03 101 116 102 106
## 35 2023-02-04 116 101 92 101
原理:计算过去一段时间窗口内的统计量。
# 创建滚动统计特征
time_data <- time_data %>%
dplyr::mutate(
# 7天滚动均值
rolling_mean_7 = zoo::rollmean(patients, k = 7, fill = NA, align = "right"),
# 7天滚动标准差
rolling_sd_7 = zoo::rollapply(patients, width = 7, FUN = sd, fill = NA, align = "right"),
# 7天滚动最大值
rolling_max_7 = zoo::rollapply(patients, width = 7, FUN = max, fill = NA, align = "right")
)
# 可视化滚动统计
ggplot2::ggplot(time_data[30:n_days, ], ggplot2::aes(x = date)) +
ggplot2::geom_line(ggplot2::aes(y = patients, color = "原始值"), alpha = 0.5) +
ggplot2::geom_line(ggplot2::aes(y = rolling_mean_7, color = "7天滚动均值"), linewidth = 1) +
ggplot2::labs(
title = "滚动统计特征",
x = "日期", y = "门诊人数", color = ""
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
# 提取时间分解特征
time_data <- time_data %>%
dplyr::mutate(
year = as.numeric(format(date, "%Y")),
month = as.numeric(format(date, "%m")),
day = as.numeric(format(date, "%d")),
weekday = as.numeric(format(date, "%w")), # 0-6, 周日为0
is_weekend = weekday %in% c(0, 6)
)
# 按星期几统计
weekday_stats <- time_data %>%
dplyr::group_by(weekday) %>%
dplyr::summarise(
mean_patients = mean(patients),
sd_patients = sd(patients),
.groups = "drop"
)
weekday_stats$weekday_name <- c("周日", "周一", "周二", "周三", "周四", "周五", "周六")
ggplot2::ggplot(weekday_stats, ggplot2::aes(x = weekday_name, y = mean_patients)) +
ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
ggplot2::geom_errorbar(ggplot2::aes(ymin = mean_patients - sd_patients,
ymax = mean_patients + sd_patients),
width = 0.2) +
ggplot2::labs(
title = "不同星期几的门诊量分布",
x = "", y = "平均门诊人数"
) +
ggplot2::theme_minimal()
数据分箱(Binning/Discretization):将连续变量转换为离散的区间(箱)。
用途: 1. 处理非线性关系 2. 减少噪声影响 3. 提高可解释性
# 创建模拟数据
set.seed(42)
binning_data <- data.frame(
age = round(runif(500, 18, 90)),
income = round(rnorm(500, mean = 50000, sd = 20000))
)
# 查看原始数据分布
p1 <- ggplot2::ggplot(binning_data, ggplot2::aes(x = age)) +
ggplot2::geom_histogram(bins = 30, fill = "steelblue", color = "white") +
ggplot2::labs(title = "年龄分布(原始)", x = "年龄", y = "频数") +
ggplot2::theme_minimal()
p2 <- ggplot2::ggplot(binning_data, ggplot2::aes(x = income)) +
ggplot2::geom_histogram(bins = 30, fill = "steelblue", color = "white") +
ggplot2::labs(title = "收入分布(原始)", x = "收入", y = "频数") +
ggplot2::theme_minimal()
p1 + p2
原理:将数据范围等分为若干区间。
# 等距分箱:将年龄分为5个等宽区间
binning_data$age_equal_width <- cut(
binning_data$age,
breaks = 5,
labels = c("青年", "中青年", "中年", "中老年", "老年")
)
# 查看分箱结果
table(binning_data$age_equal_width)
##
## 青年 中青年 中年 中老年 老年
## 114 92 110 82 102
# 可视化
ggplot2::ggplot(binning_data, ggplot2::aes(x = age_equal_width)) +
ggplot2::geom_bar(fill = "steelblue") +
ggplot2::labs(
title = "等距分箱结果",
x = "年龄区间", y = "频数"
) +
ggplot2::theme_minimal()
原理:每个箱包含大致相同数量的样本。
# 等频分箱:将年龄分为5个等频区间
binning_data$age_equal_freq <- cut(
binning_data$age,
breaks = quantile(binning_data$age, probs = seq(0, 1, 0.2)),
include.lowest = TRUE,
labels = c("Q1", "Q2", "Q3", "Q4", "Q5")
)
# 查看分箱结果
table(binning_data$age_equal_freq)
##
## Q1 Q2 Q3 Q4 Q5
## 102 104 97 101 96
# 可视化
ggplot2::ggplot(binning_data, ggplot2::aes(x = age_equal_freq)) +
ggplot2::geom_bar(fill = "steelblue") +
ggplot2::labs(
title = "等频分箱结果",
x = "年龄区间(按分位数)", y = "频数"
) +
ggplot2::theme_minimal()
原理:使用决策树找到最优分割点,使目标变量的区分度最大。
# 创建带目标变量的数据
set.seed(42)
binning_data$risk <- ifelse(
binning_data$age < 40,
sample(c(0, 1), sum(binning_data$age < 40), replace = TRUE, prob = c(0.8, 0.2)),
ifelse(
binning_data$age < 60,
sample(c(0, 1), sum(binning_data$age >= 40 & binning_data$age < 60), replace = TRUE, prob = c(0.5, 0.5)),
sample(c(0, 1), sum(binning_data$age >= 60), replace = TRUE, prob = c(0.3, 0.7))
)
)
# 使用决策树找到最优分割点
tree_model <- rpart::rpart(
risk ~ age,
data = binning_data,
method = "class",
control = rpart::rpart.control(maxdepth = 3, minbucket = 30)
)
# 查看决策树
print(tree_model)
## n= 500
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 500 242 1 (0.4840000 0.5160000)
## 2) age< 39.5 151 25 0 (0.8344371 0.1655629) *
## 3) age>=39.5 349 116 1 (0.3323782 0.6676218) *
# 可视化决策树
rpart.plot::rpart.plot(tree_model, type = 3, extra = 101)
# 根据决策树结果进行分箱
binning_data$age_tree_bin <- cut(
binning_data$age,
breaks = c(-Inf, 39.5, 59.5, Inf),
labels = c("低风险", "中风险", "高风险")
)
# 查看分箱后的风险分布
risk_by_bin <- binning_data %>%
dplyr::group_by(age_tree_bin) %>%
dplyr::summarise(
n = n(),
risk_rate = mean(risk),
.groups = "drop"
)
print(risk_by_bin)
## # A tibble: 3 × 3
## age_tree_bin n risk_rate
## <fct> <int> <dbl>
## 1 低风险 151 0.166
## 2 中风险 152 0.612
## 3 高风险 197 0.711
准确率(Accuracy):
\[Accuracy = \frac{TP + TN}{TP + TN + FP + FN}\]
精确率(Precision):
\[Precision = \frac{TP}{TP + FP}\]
召回率(Recall/Sensitivity):
\[Recall = \frac{TP}{TP + FN}\]
F1分数:
\[F1 = 2 \times \frac{Precision \times Recall}{Precision + Recall}\]
# 创建模拟预测结果
set.seed(42)
n <- 200
# 模拟真实标签和预测概率
actual <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
pred_prob <- runif(n)
predicted <- ifelse(pred_prob > 0.5, 1, 0)
# 计算混淆矩阵
conf_matrix <- table(实际 = actual, 预测 = predicted)
print(conf_matrix)
## 预测
## 实际 0 1
## 0 77 52
## 1 39 32
# 提取混淆矩阵元素
TN <- conf_matrix[1, 1]
FP <- conf_matrix[1, 2]
FN <- conf_matrix[2, 1]
TP <- conf_matrix[2, 2]
# 计算各指标
accuracy <- (TP + TN) / sum(conf_matrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
f1 <- 2 * precision * recall / (precision + recall)
# 结果汇总
metrics_df <- data.frame(
指标 = c("准确率", "精确率", "召回率", "F1分数"),
公式 = c("(TP+TN)/Total", "TP/(TP+FP)", "TP/(TP+FN)", "2*P*R/(P+R)"),
值 = c(accuracy, precision, recall, f1)
)
print(metrics_df)
## 指标 公式 值
## 1 准确率 (TP+TN)/Total 0.5450000
## 2 精确率 TP/(TP+FP) 0.3809524
## 3 召回率 TP/(TP+FN) 0.4507042
## 4 F1分数 2*P*R/(P+R) 0.4129032
# 使用caret包计算完整的混淆矩阵
confusion_result <- caret::confusionMatrix(
factor(predicted, levels = c(0, 1)),
factor(actual, levels = c(0, 1)),
positive = "1"
)
print(confusion_result)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 77 39
## 1 52 32
##
## Accuracy : 0.545
## 95% CI : (0.4733, 0.6154)
## No Information Rate : 0.645
## P-Value [Acc > NIR] : 0.9986
##
## Kappa : 0.0457
##
## Mcnemar's Test P-Value : 0.2084
##
## Sensitivity : 0.4507
## Specificity : 0.5969
## Pos Pred Value : 0.3810
## Neg Pred Value : 0.6638
## Prevalence : 0.3550
## Detection Rate : 0.1600
## Detection Prevalence : 0.4200
## Balanced Accuracy : 0.5238
##
## 'Positive' Class : 1
##
对于二分类问题,混淆矩阵包含四个基本元素:
# 可视化混淆矩阵
conf_df <- data.frame(
预测 = rep(c("预测为正", "预测为负"), each = 2),
实际 = rep(c("实际为正", "实际为负"), 2),
数值 = c(TP, FN, FP, TN)
)
conf_df$预测 <- factor(conf_df$预测, levels = c("预测为正", "预测为负"))
conf_df$实际 <- factor(conf_df$实际, levels = c("实际为正", "实际为负"))
ggplot2::ggplot(conf_df, ggplot2::aes(x = 预测, y = 实际, fill = 数值)) +
ggplot2::geom_tile(color = "white", linewidth = 1) +
ggplot2::geom_text(ggplot2::aes(label = 数值), size = 6) +
ggplot2::scale_fill_gradient(low = "lightblue", high = "steelblue") +
ggplot2::labs(
title = "混淆矩阵可视化",
x = "预测值", y = "实际值"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "none")
# 多分类混淆矩阵示例
set.seed(42)
multi_actual <- sample(c("A", "B", "C"), 200, replace = TRUE)
multi_pred <- sample(c("A", "B", "C"), 200, replace = TRUE,
prob = c(0.4, 0.35, 0.25))
# 计算多分类混淆矩阵
multi_conf <- table(实际 = multi_actual, 预测 = multi_pred)
print(multi_conf)
## 预测
## 实际 A B C
## A 28 31 16
## B 37 22 19
## C 24 11 12
# 可视化
multi_conf_df <- as.data.frame(multi_conf)
ggplot2::ggplot(multi_conf_df, ggplot2::aes(x = 预测, y = 实际, fill = Freq)) +
ggplot2::geom_tile(color = "white", linewidth = 1) +
ggplot2::geom_text(ggplot2::aes(label = Freq), size = 5) +
ggplot2::scale_fill_gradient(low = "white", high = "steelblue") +
ggplot2::labs(
title = "多分类混淆矩阵",
x = "预测类别", y = "实际类别"
) +
ggplot2::theme_minimal()
# 更多衍生指标
specificity <- TN / (TN + FP) # 特异度
npv <- TN / (TN + FN) # 阴性预测值
fpr <- FP / (TN + FP) # 假正率
fnr <- FN / (TP + FN) # 假负率
derived_metrics <- data.frame(
指标 = c("特异度(Specificity)", "阴性预测值(NPV)",
"假正率(FPR)", "假负率(FNR)"),
公式 = c("TN/(TN+FP)", "TN/(TN+FN)", "FP/(TN+FP)", "FN/(TP+FN)"),
值 = c(specificity, npv, fpr, fnr)
)
print(derived_metrics)
## 指标 公式 值
## 1 特异度(Specificity) TN/(TN+FP) 0.5968992
## 2 阴性预测值(NPV) TN/(TN+FN) 0.6637931
## 3 假正率(FPR) FP/(TN+FP) 0.4031008
## 4 假负率(FNR) FN/(TP+FN) 0.5492958
ROC曲线(Receiver Operating Characteristic Curve):展示不同阈值下真阳性率(TPR)和假阳性率(FPR)的关系。
AUC(Area Under Curve):ROC曲线下的面积,衡量模型整体分类能力。
# 使用pROC包绘制ROC曲线
set.seed(42)
n <- 200
# 模拟数据
actual <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
# 好的模型:预测概率与真实标签相关
pred_prob <- ifelse(actual == 1,
rbeta(n, 5, 2), # 正类:概率偏高
rbeta(n, 2, 5)) # 负类:概率偏低
# 计算ROC曲线
roc_obj <- pROC::roc(actual, pred_prob, quiet = TRUE)
# 绘制ROC曲线
ggplot2::ggplot(data.frame(
specificity = roc_obj$specificities,
sensitivity = roc_obj$sensitivities
), ggplot2::aes(x = 1 - specificity, y = sensitivity)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
ggplot2::annotate("text", x = 0.75, y = 0.25,
label = paste("AUC =", round(pROC::auc(roc_obj), 3)),
size = 5) +
ggplot2::labs(
title = "ROC曲线",
subtitle = "曲线越靠近左上角,模型性能越好",
x = "假阳性率(1 - 特异度)", y = "真阳性率(敏感度)"
) +
ggplot2::theme_minimal()
| AUC值 | 模型性能 |
|---|---|
| 0.5 | 随机猜测 |
| 0.5-0.7 | 较差 |
| 0.7-0.8 | 一般 |
| 0.8-0.9 | 良好 |
| 0.9-1.0 | 优秀 |
# 比较不同模型的ROC曲线
set.seed(42)
# 模拟三个不同质量的模型
# 好模型
pred_good <- ifelse(actual == 1, rbeta(n, 6, 1.5), rbeta(n, 1.5, 6))
# 中等模型
pred_medium <- ifelse(actual == 1, rbeta(n, 3, 2), rbeta(n, 2, 3))
# 差模型
pred_poor <- ifelse(actual == 1, rbeta(n, 2, 2), rbeta(n, 2, 2))
# 计算ROC
roc_good <- pROC::roc(actual, pred_good, quiet = TRUE)
roc_medium <- pROC::roc(actual, pred_medium, quiet = TRUE)
roc_poor <- pROC::roc(actual, pred_poor, quiet = TRUE)
# 创建数据框
roc_compare <- rbind(
data.frame(FPR = 1 - roc_good$specificities, TPR = roc_good$sensitivities,
Model = paste("好模型 (AUC =", round(pROC::auc(roc_good), 2), ")")),
data.frame(FPR = 1 - roc_medium$specificities, TPR = roc_medium$sensitivities,
Model = paste("中等模型 (AUC =", round(pROC::auc(roc_medium), 2), ")")),
data.frame(FPR = 1 - roc_poor$specificities, TPR = roc_poor$sensitivities,
Model = paste("差模型 (AUC =", round(pROC::auc(roc_poor), 2), ")"))
)
ggplot2::ggplot(roc_compare, ggplot2::aes(x = FPR, y = TPR, color = Model)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
ggplot2::labs(
title = "不同模型的ROC曲线对比",
x = "假阳性率", y = "真阳性率"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
PR曲线(Precision-Recall Curve):展示不同阈值下精确率和召回率的关系。
适用场景:类别不平衡时,PR曲线比ROC曲线更能反映模型性能。
# 绘制PR曲线
set.seed(42)
# 创建不平衡数据
actual_imb <- sample(c(0, 1), 500, replace = TRUE, prob = c(0.9, 0.1))
pred_imb <- ifelse(actual_imb == 1, rbeta(500, 5, 2), rbeta(500, 2, 5))
# 计算不同阈值下的精确率和召回率
thresholds <- seq(0, 1, by = 0.01)
pr_data <- data.frame()
for (thresh in thresholds) {
pred_class <- ifelse(pred_imb >= thresh, 1, 0)
TP <- sum(pred_class == 1 & actual_imb == 1)
FP <- sum(pred_class == 1 & actual_imb == 0)
FN <- sum(pred_class == 0 & actual_imb == 1)
precision <- ifelse(TP + FP == 0, 1, TP / (TP + FP))
recall <- ifelse(TP + FN == 0, 0, TP / (TP + FN))
pr_data <- rbind(pr_data, data.frame(
threshold = thresh,
precision = precision,
recall = recall
))
}
# 绘制PR曲线
ggplot2::ggplot(pr_data, ggplot2::aes(x = recall, y = precision)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_hline(yintercept = mean(actual_imb), linetype = "dashed",
color = "gray50") +
ggplot2::labs(
title = "PR曲线",
subtitle = "虚线表示随机分类器的基线(正类比例)",
x = "召回率(Recall)", y = "精确率(Precision)"
) +
ggplot2::theme_minimal()
# 对比ROC和PR曲线在不平衡数据上的表现
# 计算ROC
roc_imb <- pROC::roc(actual_imb, pred_imb, quiet = TRUE)
# 创建对比图
p1 <- ggplot2::ggplot(data.frame(
FPR = 1 - roc_imb$specificities,
TPR = roc_imb$sensitivities
), ggplot2::aes(x = FPR, y = TPR)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
ggplot2::labs(
title = paste("ROC曲线 (AUC =", round(pROC::auc(roc_imb), 2), ")"),
x = "假阳性率", y = "真阳性率"
) +
ggplot2::theme_minimal()
p2 <- ggplot2::ggplot(pr_data, ggplot2::aes(x = recall, y = precision)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_hline(yintercept = mean(actual_imb), linetype = "dashed", color = "gray50") +
ggplot2::labs(
title = "PR曲线",
x = "召回率", y = "精确率"
) +
ggplot2::theme_minimal()
p1 + p2 + patchwork::plot_annotation(
title = "类别不平衡数据:ROC vs PR曲线",
subtitle = paste("正类比例:", mean(actual_imb) * 100, "%")
)
| 指标 | 公式 | 特点 |
|---|---|---|
| MAE | \(\frac{1}{n}\sum|y_i - \hat{y}_i|\) | 对异常值稳健 |
| MSE | \(\frac{1}{n}\sum(y_i - \hat{y}_i)^2\) | 对大误差敏感 |
| RMSE | \(\sqrt{MSE}\) | 与原数据同量纲 |
| R² | \(1 - \frac{SS_{res}}{SS_{tot}}\) | 解释方差比例 |
# 回归评估指标示例
set.seed(42)
n <- 100
# 模拟数据
x <- runif(n, 0, 10)
y_true <- 2 * x + 3
y_pred <- y_true + rnorm(n, 0, 2) # 添加预测误差
# 计算各指标
residuals <- y_true - y_pred
MAE <- mean(abs(residuals))
MSE <- mean(residuals^2)
RMSE <- sqrt(MSE)
SS_res <- sum(residuals^2)
SS_tot <- sum((y_true - mean(y_true))^2)
R2 <- 1 - SS_res / SS_tot
# 调整R²(假设有1个自变量)
n <- length(y_true)
p <- 1
R2_adj <- 1 - (1 - R2) * (n - 1) / (n - p - 1)
reg_metrics <- data.frame(
指标 = c("MAE", "MSE", "RMSE", "R²", "调整R²"),
值 = c(MAE, MSE, RMSE, R2, R2_adj)
)
print(reg_metrics)
## 指标 值
## 1 MAE 1.4575998
## 2 MSE 3.4293694
## 3 RMSE 1.8518557
## 4 R² 0.9050144
## 5 调整R² 0.9040451
# 预测值 vs 实际值
plot_data <- data.frame(
实际值 = y_true,
预测值 = y_pred
)
ggplot2::ggplot(plot_data, ggplot2::aes(x = 实际值, y = 预测值)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
ggplot2::geom_smooth(method = "lm", se = FALSE, color = "steelblue") +
ggplot2::labs(
title = "预测值 vs 实际值",
subtitle = "虚线为理想预测线(y=x)",
x = "实际值", y = "预测值"
) +
ggplot2::theme_minimal()
# 残差图
resid_data <- data.frame(
预测值 = y_pred,
残差 = residuals
)
ggplot2::ggplot(resid_data, ggplot2::aes(x = 预测值, y = 残差)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
ggplot2::labs(
title = "残差图",
subtitle = "检查残差是否随机分布",
x = "预测值", y = "残差"
) +
ggplot2::theme_minimal()
对于多分类问题,需要将各类别的指标进行汇总:
| 方法 | 计算方式 | 特点 |
|---|---|---|
| 宏平均 | 各类别指标的算术平均 | 平等对待每个类别 |
| 微平均 | 汇总所有TP/FP/FN后计算 | 受多数类影响大 |
| 加权平均 | 按类别样本数加权平均 | 考虑类别不平衡 |
# 多分类评估示例
set.seed(42)
# 模拟三分类数据
n <- 300
multi_actual <- sample(c("A", "B", "C"), n, replace = TRUE, prob = c(0.5, 0.3, 0.2))
# 模拟预测(有一定准确率)
multi_pred <- sapply(multi_actual, function(x) {
if (x == "A") sample(c("A", "B", "C"), 1, prob = c(0.7, 0.2, 0.1))
else if (x == "B") sample(c("A", "B", "C"), 1, prob = c(0.2, 0.6, 0.2))
else sample(c("A", "B", "C"), 1, prob = c(0.1, 0.2, 0.7))
})
# 混淆矩阵
multi_conf <- table(实际 = multi_actual, 预测 = multi_pred)
print(multi_conf)
## 预测
## 实际 A B C
## A 112 21 18
## B 24 60 7
## C 4 10 44
# 计算每个类别的指标
calc_class_metrics <- function(conf_matrix, class_name) {
TP <- conf_matrix[class_name, class_name]
FP <- sum(conf_matrix[, class_name]) - TP
FN <- sum(conf_matrix[class_name, ]) - TP
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
f1 <- 2 * precision * recall / (precision + recall)
return(c(Precision = precision, Recall = recall, F1 = f1))
}
class_metrics <- data.frame()
for (cls in c("A", "B", "C")) {
metrics <- calc_class_metrics(multi_conf, cls)
class_metrics <- rbind(class_metrics, data.frame(
类别 = cls,
Precision = metrics["Precision"],
Recall = metrics["Recall"],
F1 = metrics["F1"]
))
}
# 计算宏平均
macro_avg <- data.frame(
类别 = "宏平均",
Precision = mean(class_metrics$Precision),
Recall = mean(class_metrics$Recall),
F1 = mean(class_metrics$F1)
)
# 计算加权平均
class_counts <- table(multi_actual)
weights <- class_counts / sum(class_counts)
weighted_avg <- data.frame(
类别 = "加权平均",
Precision = sum(class_metrics$Precision * weights),
Recall = sum(class_metrics$Recall * weights),
F1 = sum(class_metrics$F1 * weights)
)
# 汇总结果
all_metrics <- rbind(class_metrics, macro_avg, weighted_avg)
print(all_metrics)
## 类别 Precision Recall F1
## Precision A 0.8000000 0.7417219 0.7697595
## Precision1 B 0.6593407 0.6593407 0.6593407
## Precision2 C 0.6376812 0.7586207 0.6929134
## 1 宏平均 0.6990073 0.7198944 0.7073378
## 11 加权平均 0.7259517 0.7200000 0.7214088
# 使用caret计算多分类混淆矩阵
multi_confusion <- caret::confusionMatrix(
factor(multi_pred, levels = c("A", "B", "C")),
factor(multi_actual, levels = c("A", "B", "C"))
)
print(multi_confusion$byClass)
## Sensitivity Specificity Pos Pred Value Neg Pred Value Precision
## Class: A 0.7417219 0.8120805 0.8000000 0.7562500 0.8000000
## Class: B 0.6593407 0.8516746 0.6593407 0.8516746 0.6593407
## Class: C 0.7586207 0.8966942 0.6376812 0.9393939 0.6376812
## Recall F1 Prevalence Detection Rate Detection Prevalence
## Class: A 0.7417219 0.7697595 0.5033333 0.3733333 0.4666667
## Class: B 0.6593407 0.6593407 0.3033333 0.2000000 0.3033333
## Class: C 0.7586207 0.6929134 0.1933333 0.1466667 0.2300000
## Balanced Accuracy
## Class: A 0.7769012
## Class: B 0.7555077
## Class: C 0.8276575
对于时间序列数据,不能随机划分,必须按时间顺序。
# 时间序列交叉验证示意
n <- 100
time_series <- 1:n
# 创建时间序列CV的划分
ts_splits <- data.frame()
for (i in 1:5) {
train_end <- 20 * i
test_start <- train_end + 1
test_end <- min(train_end + 20, n)
ts_splits <- rbind(ts_splits, data.frame(
fold = i,
type = "训练集",
start = 1,
end = train_end
))
ts_splits <- rbind(ts_splits, data.frame(
fold = i,
type = "测试集",
start = test_start,
end = test_end
))
}
# 可视化
ts_splits$fold <- factor(ts_splits$fold)
ts_splits$type <- factor(ts_splits$type, levels = c("训练集", "测试集"))
ggplot2::ggplot(ts_splits, ggplot2::aes(x = start, xend = end, y = fold, yend = fold, color = type)) +
ggplot2::geom_segment(linewidth = 3) +
ggplot2::scale_color_manual(values = c("训练集" = "steelblue", "测试集" = "salmon")) +
ggplot2::labs(
title = "时间序列交叉验证",
subtitle = "训练集在前,测试集在后,保持时间顺序",
x = "时间", y = "折数"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
当数据有分组结构时(如同一患者的多次测量),需要确保同一组的数据不会同时出现在训练集和测试集。
# 分组交叉验证示意
set.seed(42)
# 创建分组数据(5个患者,每个患者10次测量)
group_data <- data.frame(
patient_id = rep(1:5, each = 10),
measurement = 1:50
)
# 使用caret创建分组划分
group_folds <- caret::createFolds(group_data$patient_id, k = 5)
# 可视化分组CV
group_viz <- data.frame()
for (i in 1:5) {
test_patients <- unique(group_data$patient_id)[group_folds[[i]]]
test_idx <- group_data$patient_id %in% test_patients
group_viz <- rbind(group_viz, data.frame(
fold = i,
patient = group_data$patient_id,
type = ifelse(test_idx, "测试集", "训练集")
))
}
group_viz$fold <- factor(group_viz$fold)
group_viz$type <- factor(group_viz$type, levels = c("训练集", "测试集"))
ggplot2::ggplot(group_viz, ggplot2::aes(x = patient, y = fold, fill = type)) +
ggplot2::geom_tile(color = "white", linewidth = 0.5) +
ggplot2::scale_fill_manual(values = c("训练集" = "steelblue", "测试集" = "salmon")) +
ggplot2::labs(
title = "分组交叉验证",
subtitle = "同一患者的数据要么全在训练集,要么全在测试集",
x = "患者ID", y = "折数"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
核心意义:交叉验证变体针对特殊数据结构设计,确保评估结果可靠,避免数据泄露。
举例:预测患者再入院风险
为什么需要选择正确的交叉验证方法?
患者数据有特殊结构:同一患者可能有多次住院记录,时间序列数据有先后顺序。如果用普通K折交叉验证,同一患者的数据可能同时出现在训练集和测试集,导致过于乐观的评估结果;或者用未来的数据预测过去,不符合实际应用场景。
研究目的是什么?
可能得到什么样的结果?
使用正确的交叉验证方法后,模型评估结果可能比普通CV低5-10%,但这才是真实的泛化性能。例如,普通CV显示AUC=0.85,分组CV显示AUC=0.78,后者更接近实际部署效果。
需要注意什么?
选择交叉验证方法时,要考虑数据结构和实际应用场景。如果模型将用于预测新患者,必须使用分组CV;如果用于时间序列预测,必须使用时间序列CV。
比较两个模型在多次交叉验证中的性能差异。
# 模拟两个模型的交叉验证结果
set.seed(42)
k <- 10 # 10折交叉验证
# 模型A的准确率
model_a_acc <- rnorm(k, mean = 0.85, sd = 0.03)
# 模型B的准确率(略差)
model_b_acc <- rnorm(k, mean = 0.82, sd = 0.04)
# 配对t检验
t_test_result <- t.test(model_a_acc, model_b_acc, paired = TRUE)
print(t_test_result)
##
## Paired t-test
##
## data: model_a_acc and model_b_acc
## t = 2.1428, df = 9, p-value = 0.06075
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.002950963 0.108865307
## sample estimates:
## mean difference
## 0.05295717
# 可视化比较
compare_df <- data.frame(
fold = 1:k,
模型A = model_a_acc,
模型B = model_b_acc
)
compare_long <- tidyr::pivot_longer(compare_df, cols = c(模型A, 模型B),
names_to = "模型", values_to = "准确率")
ggplot2::ggplot(compare_long, ggplot2::aes(x = fold, y = 准确率, color = 模型)) +
ggplot2::geom_point(size = 3) +
ggplot2::geom_line(ggplot2::aes(group = 模型)) +
ggplot2::labs(
title = "两模型交叉验证准确率比较",
subtitle = paste("配对t检验 p值 =", round(t_test_result$p.value, 4)),
x = "折数", y = "准确率"
) +
ggplot2::theme_minimal() +
ggplot2::theme(legend.position = "bottom")
比较两个分类模型在相同样本上的预测结果。
# McNemar检验示例
set.seed(42)
n <- 200
# 模拟两个模型的预测
actual <- sample(c(0, 1), n, replace = TRUE)
model1_pred <- ifelse(runif(n) < 0.8, actual, 1 - actual) # 80%准确率
model2_pred <- ifelse(runif(n) < 0.75, actual, 1 - actual) # 75%准确率
# 创建列联表
# b01: 模型1正确,模型2错误
# b10: 模型1错误,模型2正确
b01 <- sum(model1_pred == actual & model2_pred != actual)
b10 <- sum(model1_pred != actual & model2_pred == actual)
contingency_table <- matrix(
c(sum(model1_pred == actual & model2_pred == actual), b01,
b10, sum(model1_pred != actual & model2_pred != actual)),
nrow = 2,
dimnames = list(模型2 = c("正确", "错误"), 模型1 = c("正确", "错误"))
)
print(contingency_table)
## 模型1
## 模型2 正确 错误
## 正确 121 31
## 错误 38 10
# McNemar检验
mcnemar_result <- mcnemar.test(contingency_table)
print(mcnemar_result)
##
## McNemar's Chi-squared test with continuity correction
##
## data: contingency_table
## McNemar's chi-squared = 0.52174, df = 1, p-value = 0.4701
遍历所有参数组合,找到最优参数。
# 网格搜索示例:决策树
set.seed(42)
# 定义参数网格
tune_grid <- expand.grid(
maxdepth = c(3, 5, 7, 10),
minsplit = c(5, 10, 20)
)
# 使用caret进行网格搜索
tune_control <- caret::trainControl(
method = "cv",
number = 5
)
# 训练模型
model_tune <- caret::train(
Species ~ .,
data = iris,
method = "rpart",
trControl = tune_control,
tuneGrid = data.frame(cp = seq(0.001, 0.1, length.out = 10))
)
# 查看调优结果
print(model_tune)
## CART
##
## 150 samples
## 4 predictor
## 3 classes: 'setosa', 'versicolor', 'virginica'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 120, 120, 120, 120, 120
## Resampling results across tuning parameters:
##
## cp Accuracy Kappa
## 0.001 0.94 0.91
## 0.012 0.94 0.91
## 0.023 0.94 0.91
## 0.034 0.94 0.91
## 0.045 0.94 0.91
## 0.056 0.94 0.91
## 0.067 0.94 0.91
## 0.078 0.94 0.91
## 0.089 0.94 0.91
## 0.100 0.94 0.91
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.1.
# 可视化调优过程
ggplot2::ggplot(model_tune$results, ggplot2::aes(x = cp, y = Accuracy)) +
ggplot2::geom_point(size = 2) +
ggplot2::geom_line() +
ggplot2::geom_errorbar(ggplot2::aes(ymin = Accuracy - AccuracySD,
ymax = Accuracy + AccuracySD),
width = 0.005) +
ggplot2::labs(
title = "网格搜索:不同cp参数的模型性能",
x = "复杂度参数(cp)", y = "准确率"
) +
ggplot2::theme_minimal()
随机搜索(Random Search)
随机采样参数组合,更高效。
# 随机搜索
set.seed(42)
random_control <- caret::trainControl(
method = "cv",
number = 5,
search = "random"
)
model_random <- caret::train(
Species ~ .,
data = iris,
method = "rpart",
trControl = random_control,
tuneLength = 10 # 随机尝试10组参数
)
print(model_random)
## CART
##
## 150 samples
## 4 predictor
## 3 classes: 'setosa', 'versicolor', 'virginica'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 120, 120, 120, 120, 120
## Resampling results across tuning parameters:
##
## cp Accuracy Kappa
## 0.00 0.9400000 0.91
## 0.44 0.7600000 0.64
## 0.50 0.3333333 0.00
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.
使用贝叶斯方法智能搜索最优参数。
# 贝叶斯优化示意(概念演示)
# 实际应用可使用rBayesianOptimization包
# 创建模拟的目标函数
set.seed(42)
objective_function <- function(x) {
# 假设最优参数在x=0.7附近
-((x - 0.7)^2) + rnorm(1, 0, 0.1)
}
# 模拟贝叶斯优化的迭代过程
iterations <- 10
results <- data.frame()
for (i in 1:iterations) {
# 在早期随机探索,后期利用已知信息
if (i <= 3) {
x <- runif(1)
} else {
# 在已知最优点附近探索
best_x <- results$x[which.max(results$score)]
x <- best_x + rnorm(1, 0, 0.1)
x <- max(0, min(1, x)) # 限制在[0,1]范围
}
score <- objective_function(x)
results <- rbind(results, data.frame(iteration = i, x = x, score = score))
}
# 可视化优化过程
ggplot2::ggplot(results, ggplot2::aes(x = iteration, y = score)) +
ggplot2::geom_point(size = 2, color = "steelblue") +
ggplot2::geom_line(color = "steelblue") +
ggplot2::labs(
title = "贝叶斯优化过程示意",
subtitle = "逐步逼近最优参数",
x = "迭代次数", y = "目标函数值"
) +
ggplot2::theme_minimal()
线性回归是机器学习中最基础的模型,假设目标变量与特征之间存在线性关系:
\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p + \epsilon\]
其中: - \(\beta_0\) 是截距 - \(\beta_1, \beta_2, ..., \beta_p\) 是回归系数 - \(\epsilon\) 是误差项
目标:最小化残差平方和(OLS):
\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]
# 线性回归示例:医疗费用预测
set.seed(42)
n <- 200
# 模拟数据
age <- round(runif(n, 18, 80))
bmi <- rnorm(n, 25, 5)
smoker <- sample(c(0, 1), n, replace = TRUE, prob = c(0.8, 0.2))
# 医疗费用与各因素相关
medical_cost <- 1000 + 50 * age + 200 * bmi + 5000 * smoker + rnorm(n, 0, 1000)
# 创建数据框
medical_data <- data.frame(
age = age,
bmi = bmi,
smoker = factor(smoker, levels = c(0, 1), labels = c("非吸烟", "吸烟")),
cost = medical_cost
)
# 拟合线性回归模型
lm_model <- lm(cost ~ age + bmi + smoker, data = medical_data)
# 查看模型摘要
summary(lm_model)
##
## Call:
## lm(formula = cost ~ age + bmi + smoker, data = medical_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2418.3 -631.9 -11.5 720.7 3277.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 980.584 441.075 2.223 0.0273 *
## age 49.585 4.153 11.939 <2e-16 ***
## bmi 199.889 15.498 12.898 <2e-16 ***
## smoker吸烟 4922.772 191.751 25.673 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1048 on 196 degrees of freedom
## Multiple R-squared: 0.8437, Adjusted R-squared: 0.8413
## F-statistic: 352.7 on 3 and 196 DF, p-value: < 2.2e-16
# 提取回归系数
coef_df <- data.frame(
变量 = names(coef(lm_model)),
系数 = round(coef(lm_model), 2),
标准误 = round(summary(lm_model)$coefficients[, 2], 2),
t值 = round(summary(lm_model)$coefficients[, 3], 2),
p值 = round(summary(lm_model)$coefficients[, 4], 4)
)
print(coef_df)
## 变量 系数 标准误 t值 p值
## (Intercept) (Intercept) 980.58 441.07 2.22 0.0273
## age age 49.59 4.15 11.94 0.0000
## bmi bmi 199.89 15.50 12.90 0.0000
## smoker吸烟 smoker吸烟 4922.77 191.75 25.67 0.0000
# 预测值 vs 实际值
medical_data$predicted <- predict(lm_model)
ggplot2::ggplot(medical_data, ggplot2::aes(x = cost, y = predicted, color = smoker)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
ggplot2::labs(
title = "线性回归:预测值 vs 实际值",
x = "实际医疗费用", y = "预测医疗费用"
) +
ggplot2::theme_minimal()
# 残差诊断
medical_data$residuals <- residuals(lm_model)
ggplot2::ggplot(medical_data, ggplot2::aes(x = predicted, y = residuals)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
ggplot2::geom_smooth(method = "loess", se = FALSE, color = "steelblue") +
ggplot2::labs(
title = "残差图",
subtitle = "检查残差是否随机分布(无模式)",
x = "预测值", y = "残差"
) +
ggplot2::theme_minimal()
逻辑回归用于二分类问题,通过sigmoid函数将线性组合映射到[0,1]区间:
\[P(y=1|x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + ... + \beta_p x_p)}}\]
对数几率(Log-odds):
\[\log\frac{P(y=1)}{P(y=0)} = \beta_0 + \beta_1 x_1 + ... + \beta_p x_p\]
# 逻辑回归示例:疾病诊断
set.seed(42)
n <- 300
# 模拟数据
age <- round(runif(n, 30, 80))
bmi <- rnorm(n, 26, 4)
blood_pressure <- rnorm(n, 120, 15)
# 疾病风险与因素相关
log_odds <- -10 + 0.1 * age + 0.15 * bmi + 0.03 * blood_pressure
prob <- 1 / (1 + exp(-log_odds))
disease <- rbinom(n, 1, prob)
# 创建数据框
disease_data <- data.frame(
age = age,
bmi = bmi,
blood_pressure = blood_pressure,
disease = factor(disease, levels = c(0, 1), labels = c("健康", "患病"))
)
# 拟合逻辑回归模型
logit_model <- glm(disease ~ age + bmi + blood_pressure,
data = disease_data,
family = binomial(link = "logit"))
# 查看模型摘要
summary(logit_model)
##
## Call:
## glm(formula = disease ~ age + bmi + blood_pressure, family = binomial(link = "logit"),
## data = disease_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.70130 2.59968 -2.578 0.00994 **
## age 0.06596 0.01692 3.897 9.73e-05 ***
## bmi 0.16000 0.06133 2.609 0.00909 **
## blood_pressure 0.01355 0.01424 0.952 0.34124
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 195.05 on 299 degrees of freedom
## Residual deviance: 169.26 on 296 degrees of freedom
## AIC: 177.26
##
## Number of Fisher Scoring iterations: 6
# 预测概率
disease_data$prob <- predict(logit_model, type = "response")
# 默认决策边界:0.5
disease_data$pred_class <- ifelse(disease_data$prob >= 0.5, "患病", "健康")
# 可视化概率分布
ggplot2::ggplot(disease_data, ggplot2::aes(x = prob, fill = disease)) +
ggplot2::geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
ggplot2::geom_vline(xintercept = 0.5, linetype = "dashed", color = "red") +
ggplot2::labs(
title = "逻辑回归:预测概率分布",
subtitle = "虚线为默认决策边界(0.5)",
x = "预测概率", y = "频数", fill = "实际状态"
) +
ggplot2::theme_minimal()
# 决策边界可视化(以年龄和BMI为例)
ggplot2::ggplot(disease_data, ggplot2::aes(x = age, y = bmi, color = disease)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::labs(
title = "决策边界示意",
x = "年龄", y = "BMI", color = "实际状态"
) +
ggplot2::theme_minimal()
# 计算混淆矩阵
conf_matrix <- table(实际 = disease_data$disease, 预测 = disease_data$pred_class)
print(conf_matrix)
## 预测
## 实际 患病
## 健康 30
## 患病 270
# 计算ROC和AUC
roc_obj <- pROC::roc(as.numeric(disease_data$disease) - 1,
disease_data$prob, quiet = TRUE)
ggplot2::ggplot(data.frame(
FPR = 1 - roc_obj$specificities,
TPR = roc_obj$sensitivities
), ggplot2::aes(x = FPR, y = TPR)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
ggplot2::annotate("text", x = 0.75, y = 0.25,
label = paste("AUC =", round(pROC::auc(roc_obj), 3)),
size = 5) +
ggplot2::labs(
title = "逻辑回归ROC曲线",
x = "假阳性率", y = "真阳性率"
) +
ggplot2::theme_minimal()
当特征数量多或存在多重共线性时,普通线性回归容易过拟合。正则化通过添加惩罚项来约束系数。
岭回归(L2正则化):
\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} \beta_j^2\]
LASSO回归(L1正则化):
\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} |\beta_j|\]
弹性网(Elastic Net):
\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \sum_{j=1}^{p} \beta_j^2\]
# 正则化示例
set.seed(42)
n <- 200
p <- 50 # 50个特征
# 模拟高维数据
X <- matrix(rnorm(n * p), n, p)
# 只有前5个特征有真实效应
true_beta <- c(3, 2, 1.5, 1, 0.5, rep(0, p - 5))
y <- X %*% true_beta + rnorm(n, 0, 1)
# 创建数据框
reg_data <- data.frame(y = y, X)
colnames(reg_data) <- c("y", paste0("X", 1:p))
# 使用caret进行正则化回归
# 岭回归(alpha = 0)
ridge_model <- caret::train(
y ~ .,
data = reg_data,
method = "glmnet",
trControl = caret::trainControl(method = "cv", number = 5),
tuneGrid = expand.grid(alpha = 0, lambda = 10^seq(-3, 1, length.out = 20))
)
# LASSO回归(alpha = 1)
lasso_model <- caret::train(
y ~ .,
data = reg_data,
method = "glmnet",
trControl = caret::trainControl(method = "cv", number = 5),
tuneGrid = expand.grid(alpha = 1, lambda = 10^seq(-3, 1, length.out = 20))
)
# 弹性网(alpha = 0.5)
elastic_model <- caret::train(
y ~ .,
data = reg_data,
method = "glmnet",
trControl = caret::trainControl(method = "cv", number = 5),
tuneGrid = expand.grid(alpha = 0.5, lambda = 10^seq(-3, 1, length.out = 20))
)
# 比较模型性能
model_compare <- data.frame(
模型 = c("岭回归", "LASSO", "弹性网"),
最优lambda = c(ridge_model$bestTune$lambda,
lasso_model$bestTune$lambda,
elastic_model$bestTune$lambda),
RMSE = c(min(ridge_model$results$RMSE),
min(lasso_model$results$RMSE),
min(elastic_model$results$RMSE))
)
print(model_compare)
## 模型 最优lambda RMSE
## 1 岭回归 0.2069138 1.262869
## 2 LASSO 0.1274275 1.078834
## 3 弹性网 0.2069138 1.098975
# 提取LASSO系数
lasso_coef <- coef(lasso_model$finalModel, lasso_model$bestTune$lambda)
lasso_coef_df <- data.frame(
特征 = rownames(lasso_coef),
系数 = as.numeric(lasso_coef)
)
lasso_coef_df <- lasso_coef_df[lasso_coef_df$特征 != "(Intercept)", ]
# 可视化非零系数
nonzero_coef <- lasso_coef_df[lasso_coef_df$系数 != 0, ]
nonzero_coef <- nonzero_coef[order(abs(nonzero_coef$系数), decreasing = TRUE), ]
ggplot2::ggplot(head(nonzero_coef, 15), ggplot2::aes(x = reorder(特征, abs(系数)), y = 系数)) +
ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
ggplot2::coord_flip() +
ggplot2::labs(
title = "LASSO选择的特征(非零系数)",
x = "特征", y = "系数值"
) +
ggplot2::theme_minimal()
支持向量机(SVM)寻找最优超平面,使两类样本之间的间隔最大化。
线性可分SVM:
\[\min_{w,b} \frac{1}{2}\|w\|^2\]
约束条件:\(y_i(w^T x_i + b) \geq 1\)
# 线性SVM示例
set.seed(42)
# 创建线性可分数据
n <- 100
x1 <- c(rnorm(n/2, mean = 2, sd = 0.5), rnorm(n/2, mean = 4, sd = 0.5))
x2 <- c(rnorm(n/2, mean = 2, sd = 0.5), rnorm(n/2, mean = 4, sd = 0.5))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))
svm_data <- data.frame(x1 = x1, x2 = x2, y = y)
# 训练线性SVM
svm_linear <- e1071::svm(y ~ x1 + x2,
data = svm_data,
type = "C-classification",
kernel = "linear",
cost = 1)
# 可视化决策边界
plot(svm_linear, svm_data, x1 ~ x2,
main = "线性SVM决策边界",
xlab = "特征1", ylab = "特征2")
硬间隔:要求所有样本正确分类,容易过拟合。
软间隔:允许部分样本误分类,引入松弛变量\(\xi_i\):
\[\min_{w,b,\xi} \frac{1}{2}\|w\|^2 + C \sum_{i=1}^{n} \xi_i\]
C参数:控制误分类的惩罚程度。
# 比较不同C值
set.seed(42)
# 创建有噪声的数据
n <- 100
x1 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
x2 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))
noise_data <- data.frame(x1 = x1, x2 = x2, y = y)
# 不同C值的SVM
svm_small_c <- e1071::svm(y ~ x1 + x2, data = noise_data,
kernel = "linear", cost = 0.1)
svm_large_c <- e1071::svm(y ~ x1 + x2, data = noise_data,
kernel = "linear", cost = 100)
# 比较支持向量数量
sv_compare <- data.frame(
C值 = c("0.1(软间隔)", "100(硬间隔)"),
支持向量数 = c(length(svm_small_c$index), length(svm_large_c$index))
)
print(sv_compare)
## C值 支持向量数
## 1 0.1(软间隔) 37
## 2 100(硬间隔) 15
感知机(Perceptron)是最简单的神经网络,用于线性可分数据的二分类:
\[f(x) = \text{sign}(w^T x + b)\]
学习规则:对于误分类样本\((x_i, y_i)\):
\[w \leftarrow w + \eta y_i x_i\]
# 感知机示例
set.seed(42)
# 创建线性可分数据
n <- 100
perceptron_data <- data.frame(
x1 = c(rnorm(n/2, 2, 0.5), rnorm(n/2, 4, 0.5)),
x2 = c(rnorm(n/2, 2, 0.5), rnorm(n/2, 4, 0.5)),
y = factor(c(rep(-1, n/2), rep(1, n/2)))
)
# 使用nnet包的单层感知机
perceptron_model <- nnet::nnet(y ~ x1 + x2,
data = perceptron_data,
size = 0, # 无隐藏层
skip = TRUE, # 直接连接
maxit = 100,
trace = FALSE)
# 预测
perceptron_data$pred <- predict(perceptron_model, type = "class")
# 计算准确率
accuracy <- mean(perceptron_data$pred == perceptron_data$y)
print(paste("感知机准确率:", round(accuracy, 3)))
## [1] "感知机准确率: 1"
感知机只能解决线性可分问题,对于异或(XOR)等问题无法求解。
# XOR问题示例
xor_data <- data.frame(
x1 = c(0, 0, 1, 1),
x2 = c(0, 1, 0, 1),
y = factor(c(0, 1, 1, 0))
)
ggplot2::ggplot(xor_data, ggplot2::aes(x = x1, y = x2, color = y)) +
ggplot2::geom_point(size = 5) +
ggplot2::labs(
title = "XOR问题:线性不可分",
subtitle = "感知机无法找到线性决策边界",
x = "x1", y = "x2"
) +
ggplot2::theme_minimal()
核心意义:感知机是神经网络的基础,理解感知机有助于理解深度学习的原理。
举例:简单医学诊断规则
为什么学习感知机?
虽然感知机在实际应用中很少单独使用,但它是理解神经网络的基础。感知机的学习规则(权重更新)是现代深度学习的核心思想。在医学领域,一些简单的线性可分问题可以用感知机快速解决,如根据体温和白细胞计数初步判断是否感染。
研究目的是什么?
可能得到什么样的结果?
对于线性可分数据,感知机可以找到完美分类边界;对于线性不可分数据,感知机无法收敛,需要使用多层感知机或其他非线性方法。
需要注意什么?
感知机对初始权重敏感,可能收敛到不同的解。对于实际医学问题,建议使用更强大的方法如逻辑回归或SVM,它们能处理线性不可分数据并提供概率输出。
线性判别分析(LDA)是一种生成式分类方法,假设各类数据服从高斯分布,寻找使类间距离最大、类内距离最小的投影方向。
Fisher判别准则:
\[J(w) = \frac{w^T S_B w}{w^T S_W w}\]
其中\(S_B\)是类间散度矩阵,\(S_W\)是类内散度矩阵。
# LDA示例
set.seed(42)
# 使用iris数据集(取两个类别)
iris_binary <- iris[iris$Species %in% c("setosa", "versicolor"), ]
iris_binary$Species <- factor(iris_binary$Species)
# 训练LDA模型
lda_model <- MASS::lda(Species ~ Sepal.Length + Sepal.Width +
Petal.Length + Petal.Width,
data = iris_binary)
# 查看模型结果
print(lda_model)
## Call:
## lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
## data = iris_binary)
##
## Prior probabilities of groups:
## setosa versicolor
## 0.5 0.5
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.006 3.428 1.462 0.246
## versicolor 5.936 2.770 4.260 1.326
##
## Coefficients of linear discriminants:
## LD1
## Sepal.Length -0.300458
## Sepal.Width -1.773845
## Petal.Length 2.142260
## Petal.Width 3.035726
# 预测
lda_pred <- predict(lda_model, iris_binary)
# 混淆矩阵
table(实际 = iris_binary$Species, 预测 = lda_pred$class)
## 预测
## 实际 setosa versicolor
## setosa 50 0
## versicolor 0 50
# LDA投影可视化
lda_values <- predict(lda_model, iris_binary)$x
iris_binary$LD1 <- lda_values[, 1]
ggplot2::ggplot(iris_binary, ggplot2::aes(x = LD1, fill = Species)) +
ggplot2::geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
ggplot2::labs(
title = "LDA投影结果",
subtitle = "将多维数据投影到一维判别方向",
x = "线性判别函数值", y = "频数"
) +
ggplot2::theme_minimal()
决策树通过一系列规则将数据递归划分,形成树状结构:
主要算法:
| 算法 | 分裂准则 | 特点 |
|---|---|---|
| ID3 | 信息增益 | 只能处理分类特征 |
| C4.5 | 增益率 | 可处理连续特征 |
| CART | 基尼系数 | 可用于分类和回归 |
# 决策树示例:疾病诊断
set.seed(42)
# 使用模拟数据
n <- 300
tree_data <- data.frame(
age = round(runif(n, 20, 80)),
bmi = rnorm(n, 25, 5),
blood_pressure = rnorm(n, 120, 20),
cholesterol = factor(sample(c("正常", "偏高", "高"), n, replace = TRUE, prob = c(0.5, 0.3, 0.2)))
)
# 生成目标变量(疾病风险)
risk_score <- 0.02 * tree_data$age + 0.05 * tree_data$bmi +
0.01 * tree_data$blood_pressure - 5
risk_score <- risk_score + ifelse(tree_data$cholesterol == "高", 2,
ifelse(tree_data$cholesterol == "偏高", 1, 0))
prob <- 1 / (1 + exp(-risk_score))
tree_data$disease <- factor(ifelse(runif(n) < prob, "患病", "健康"))
# 训练决策树(CART算法)
tree_model <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
data = tree_data,
method = "class",
control = rpart::rpart.control(maxdepth = 4))
# 查看树结构
print(tree_model)
## n= 300
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 300 93 健康 (0.3100000 0.6900000)
## 2) cholesterol=高,偏高 144 63 健康 (0.4375000 0.5625000)
## 4) age>=46.5 89 41 患病 (0.5393258 0.4606742)
## 8) bmi>=19.53831 77 32 患病 (0.5844156 0.4155844) *
## 9) bmi< 19.53831 12 3 健康 (0.2500000 0.7500000) *
## 5) age< 46.5 55 15 健康 (0.2727273 0.7272727) *
## 3) cholesterol=正常 156 30 健康 (0.1923077 0.8076923)
## 6) bmi>=32.1859 9 3 患病 (0.6666667 0.3333333) *
## 7) bmi< 32.1859 147 24 健康 (0.1632653 0.8367347) *
# 可视化决策树
rpart.plot::rpart.plot(tree_model,
main = "决策树:疾病诊断",
extra = 104, # 显示每个节点的类别比例
box.palette = "GnBu")
信息增益(Information Gain):
\[IG(D, A) = H(D) - \sum_{v \in Values(A)} \frac{|D_v|}{|D|} H(D_v)\]
其中\(H(D)\)是熵:
\[H(D) = -\sum_{k=1}^{K} p_k \log_2 p_k\]
基尼系数(Gini Index):
\[Gini(D) = 1 - \sum_{k=1}^{K} p_k^2\]
# 比较不同分裂准则
# 基尼系数(默认)
tree_gini <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
data = tree_data, method = "class")
# 信息增益
tree_info <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
data = tree_data, method = "class",
parms = list(split = "information"))
# 比较变量重要性
importance_compare <- data.frame(
变量 = names(tree_gini$variable.importance),
基尼系数 = tree_gini$variable.importance,
信息增益 = tree_info$variable.importance
)
# 标准化
importance_compare$基尼系数 <- importance_compare$基尼系数 / max(importance_compare$基尼系数)
importance_compare$信息增益 <- importance_compare$信息增益 / max(importance_compare$信息增益)
importance_long <- tidyr::pivot_longer(importance_compare,
cols = c(基尼系数, 信息增益),
names_to = "准则", values_to = "重要性")
ggplot2::ggplot(importance_long, ggplot2::aes(x = 变量, y = 重要性, fill = 准则)) +
ggplot2::geom_bar(stat = "identity", position = "dodge") +
ggplot2::labs(
title = "不同分裂准则的变量重要性",
x = "变量", y = "标准化重要性"
) +
ggplot2::theme_minimal() +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))
预剪枝:在构建过程中限制树的生长。
# 预剪枝参数
pre_prune <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
data = tree_data,
method = "class",
control = rpart::rpart.control(
maxdepth = 3, # 最大深度
minsplit = 20, # 分裂所需最小样本数
minbucket = 10 # 叶节点最小样本数
))
print(paste("预剪枝树的叶节点数:", length(unique(pre_prune$where))))
## [1] "预剪枝树的叶节点数: 5"
后剪枝:先构建完整树,再根据复杂度参数剪枝。
# 后剪枝:使用cp参数
printcp(tree_model)
##
## Classification tree:
## rpart::rpart(formula = disease ~ age + bmi + blood_pressure +
## cholesterol, data = tree_data, method = "class", control = rpart::rpart.control(maxdepth = 4))
##
## Variables actually used in tree construction:
## [1] age bmi cholesterol
##
## Root node error: 93/300 = 0.31
##
## n= 300
##
## CP nsplit rel error xerror xstd
## 1 0.037634 0 1.00000 1.0000 0.086136
## 2 0.032258 3 0.86022 1.1183 0.088634
## 3 0.010000 4 0.82796 1.0860 0.088012
# 可视化cp与误差的关系
plotcp(tree_model)
# 选择最优cp
best_cp <- tree_model$cptable[which.min(tree_model$cptable[, "xerror"]), "CP"]
# 剪枝
pruned_tree <- rpart::prune(tree_model, cp = best_cp)
print(paste("剪枝前叶节点数:", length(unique(tree_model$where))))
## [1] "剪枝前叶节点数: 5"
print(paste("剪枝后叶节点数:", length(unique(pruned_tree$where))))
## [1] "剪枝后叶节点数: 1"
Bagging(Bootstrap Aggregating):通过自助采样训练多个模型,再聚合结果。
随机森林在Bagging基础上增加特征随机性:
# 随机森林示例
set.seed(42)
# 训练随机森林
rf_model <- randomForest::randomForest(
disease ~ age + bmi + blood_pressure + cholesterol,
data = tree_data,
ntree = 500, # 树的数量
mtry = 2, # 每次分裂考虑的特征数
importance = TRUE
)
# 查看模型结果
print(rf_model)
##
## Call:
## randomForest(formula = disease ~ age + bmi + blood_pressure + cholesterol, data = tree_data, ntree = 500, mtry = 2, importance = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 32%
## Confusion matrix:
## 患病 健康 class.error
## 患病 29 64 0.6881720
## 健康 32 175 0.1545894
OOB(Out-of-Bag)误差:每棵树约有1/3的样本未被用于训练,可用作验证集。
# OOB误差随树数量的变化
oob_error <- data.frame(
ntree = 1:500,
OOB误差 = rf_model$err.rate[, "OOB"]
)
ggplot2::ggplot(oob_error, ggplot2::aes(x = ntree, y = OOB误差)) +
ggplot2::geom_line(color = "steelblue") +
ggplot2::labs(
title = "随机森林OOB误差",
subtitle = "误差随树数量增加趋于稳定",
x = "树的数量", y = "OOB误差率"
) +
ggplot2::theme_minimal()
# 变量重要性
importance_df <- data.frame(
变量 = rownames(rf_model$importance),
平均精度下降 = rf_model$importance[, "MeanDecreaseAccuracy"],
基尼系数下降 = rf_model$importance[, "MeanDecreaseGini"]
)
# 可视化
importance_long <- tidyr::pivot_longer(importance_df,
cols = c(平均精度下降, 基尼系数下降),
names_to = "指标", values_to = "值")
ggplot2::ggplot(importance_long, ggplot2::aes(x = reorder(变量, 值), y = 值, fill = 指标)) +
ggplot2::geom_bar(stat = "identity", position = "dodge") +
ggplot2::coord_flip() +
ggplot2::labs(
title = "随机森林变量重要性",
x = "变量", y = "重要性"
) +
ggplot2::theme_minimal()
梯度提升:通过迭代训练决策树,每棵树拟合前一棵树的残差。
\[F_m(x) = F_{m-1}(x) + \eta \cdot h_m(x)\]
其中\(h_m(x)\)是第\(m\)棵树,\(\eta\)是学习率。
# GBDT示例
set.seed(42)
# 为gbm准备数据(将因子转换为数值)
gbm_data <- tree_data
gbm_data$cholesterol_num <- as.numeric(gbm_data$cholesterol)
# 训练GBDT模型
gbm_model <- gbm::gbm(
as.numeric(disease == "患病") ~ age + bmi + blood_pressure + cholesterol_num,
data = gbm_data,
distribution = "bernoulli",
n.trees = 500,
interaction.depth = 4,
shrinkage = 0.01,
cv.folds = 5,
n.cores = 1
)
## CV: 1
## CV: 2
## CV: 3
## CV: 4
## CV: 5
# 查看最优树数量
best_iter <- gbm::gbm.perf(gbm_model, method = "cv")
print(paste("最优迭代次数:", best_iter))
## [1] "最优迭代次数: 122"
# 变量重要性
gbm_summary <- summary(gbm_model, n.trees = best_iter, plotit = FALSE)
ggplot2::ggplot(gbm_summary, ggplot2::aes(x = reorder(var, rel.inf), y = rel.inf)) +
ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
ggplot2::coord_flip() +
ggplot2::labs(
title = "GBDT变量重要性",
x = "变量", y = "相对影响"
) +
ggplot2::theme_minimal()
XGBoost在GBDT基础上引入:
目标函数:
\[Obj = \sum_{i=1}^{n} L(y_i, \hat{y}_i) + \sum_{k=1}^{K} \Omega(f_k)\]
# XGBoost示例
set.seed(42)
# 准备数据
library(xgboost)
# 转换为数值矩阵
tree_data_num <- tree_data
tree_data_num$cholesterol <- as.numeric(factor(tree_data_num$cholesterol)) - 1
tree_data_num$disease <- as.numeric(tree_data_num$disease == "患病")
# 创建xgb矩阵
xgb_matrix <- xgboost::xgb.DMatrix(
data = as.matrix(tree_data_num[, c("age", "bmi", "blood_pressure", "cholesterol")]),
label = tree_data_num$disease
)
# 训练XGBoost
xgb_model <- xgboost::xgb.train(
data = xgb_matrix,
objective = "binary:logistic",
nrounds = 100,
max_depth = 4,
eta = 0.1,
subsample = 0.8,
colsample_bytree = 0.8,
verbose = 0
)
# 预测
xgb_pred <- predict(xgb_model, xgb_matrix)
xgb_pred_class <- ifelse(xgb_pred > 0.5, "患病", "健康")
# 准确率
accuracy <- mean(xgb_pred_class == tree_data$disease)
print(paste("XGBoost准确率:", round(accuracy, 3)))
## [1] "XGBoost准确率: 0.913"
# 变量重要性
xgb_importance <- xgboost::xgb.importance(
model = xgb_model,
feature_names = c("age", "bmi", "blood_pressure", "cholesterol")
)
ggplot2::ggplot(xgb_importance, ggplot2::aes(x = reorder(Feature, Gain), y = Gain)) +
ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
ggplot2::coord_flip() +
ggplot2::labs(
title = "XGBoost变量重要性",
x = "特征", y = "增益"
) +
ggplot2::theme_minimal()
堆叠(Stacking):用多个基模型的预测结果作为新特征,训练元模型。
# Stacking示例
set.seed(42)
# 划分训练集和测试集
train_idx <- sample(1:nrow(tree_data), 0.7 * nrow(tree_data))
train_data <- tree_data[train_idx, ]
test_data <- tree_data[-train_idx, ]
# 训练多个基模型
# 模型1:决策树
model_tree <- rpart::rpart(disease ~ ., data = train_data, method = "class")
# 模型2:随机森林
model_rf <- randomForest::randomForest(disease ~ ., data = train_data)
# 模型3:逻辑回归
model_logit <- glm(disease ~ age + bmi + blood_pressure +
as.numeric(cholesterol),
data = train_data, family = binomial)
# 获取基模型预测概率
pred_tree <- predict(model_tree, test_data, type = "prob")[, 2]
pred_rf <- predict(model_rf, test_data, type = "prob")[, 2]
pred_logit <- predict(model_logit, test_data, type = "response")
# 创建元特征
meta_features <- data.frame(
tree_prob = pred_tree,
rf_prob = pred_rf,
logit_prob = pred_logit,
disease = test_data$disease
)
# 训练元模型(逻辑回归)
meta_model <- glm(disease ~ tree_prob + rf_prob + logit_prob,
data = meta_features, family = binomial)
# 最终预测
final_pred <- predict(meta_model, meta_features, type = "response")
final_pred_class <- ifelse(final_pred > 0.5, "患病", "健康")
# 评估
stacking_accuracy <- mean(final_pred_class == test_data$disease)
print(paste("Stacking准确率:", round(stacking_accuracy, 3)))
## [1] "Stacking准确率: 0.2"
# 投票分类器
# 硬投票
hard_vote <- ifelse(
(predict(model_tree, test_data, type = "class") == "患病") +
(predict(model_rf, test_data, type = "class") == "患病") +
(ifelse(pred_logit > 0.5, "患病", "健康") == "患病") >= 2,
"患病", "健康"
)
# 软投票
soft_vote_prob <- (pred_tree + pred_rf + pred_logit) / 3
soft_vote <- ifelse(soft_vote_prob > 0.5, "患病", "健康")
# 比较
voting_compare <- data.frame(
方法 = c("决策树", "随机森林", "逻辑回归", "硬投票", "软投票"),
准确率 = c(
mean(predict(model_tree, test_data, type = "class") == test_data$disease),
mean(predict(model_rf, test_data, type = "class") == test_data$disease),
mean(ifelse(pred_logit > 0.5, "患病", "健康") == test_data$disease),
mean(hard_vote == test_data$disease),
mean(soft_vote == test_data$disease)
)
)
print(voting_compare)
## 方法 准确率
## 1 决策树 0.7555556
## 2 随机森林 0.6888889
## 3 逻辑回归 0.2333333
## 4 硬投票 0.7111111
## 5 软投票 0.2666667
最大间隔分类器寻找能够正确划分两类样本且间隔最大的超平面。
间隔:超平面到最近样本点的距离。
支持向量:距离超平面最近的样本点,决定超平面位置。
优化目标:
\[\max_{w,b} \gamma \quad \text{s.t.} \quad y_i(w^T x_i + b) \geq \gamma\]
等价于:
\[\min_{w,b} \frac{1}{2}\|w\|^2 \quad \text{s.t.} \quad y_i(w^T x_i + b) \geq 1\]
# 最大间隔分类器示意
set.seed(42)
# 创建线性可分数据
n <- 50
x1 <- c(rnorm(n/2, mean = 1, sd = 0.5), rnorm(n/2, mean = 3, sd = 0.5))
x2 <- c(rnorm(n/2, mean = 1, sd = 0.5), rnorm(n/2, mean = 3, sd = 0.5))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))
svm_data <- data.frame(x1 = x1, x2 = x2, y = y)
# 训练硬间隔SVM
hard_svm <- e1071::svm(y ~ x1 + x2,
data = svm_data,
type = "C-classification",
kernel = "linear",
cost = 1e5) # 大cost近似硬间隔
# 可视化
plot(hard_svm, svm_data, x1 ~ x2,
main = "最大间隔分类器",
xlab = "特征1", ylab = "特征2",
svSymbol = 16) # 支持向量用方形标记
当数据线性不可分时,可以通过映射\(\phi(x)\)将数据映射到高维空间,使其线性可分。
核技巧:不需要显式计算\(\phi(x)\),直接计算核函数:
\[K(x_i, x_j) = \phi(x_i)^T \phi(x_j)\]
| 核函数 | 公式 | 特点 |
|---|---|---|
| 线性核 | \(K(x, z) = x^T z\) | 线性可分数据 |
| 多项式核 | \(K(x, z) = (\gamma x^T z + r)^d\) | 多项式边界 |
| RBF核 | \(K(x, z) = \exp(-\gamma\|x-z\|^2)\) | 通用、最常用 |
| Sigmoid核 | \(K(x, z) = \tanh(\gamma x^T z + r)\) | 类神经网络 |
# 不同核函数比较
set.seed(42)
# 创建非线性可分数据(环形数据)
n <- 200
theta <- runif(n, 0, 2 * pi)
r <- ifelse(sample(c(0, 1), n, replace = TRUE) == 0,
runif(n, 0, 0.5), runif(n, 1, 1.5))
x1 <- r * cos(theta)
x2 <- r * sin(theta)
y <- factor(ifelse(r < 0.7, 0, 1))
nonlinear_data <- data.frame(x1 = x1, x2 = x2, y = y)
# 可视化原始数据
ggplot2::ggplot(nonlinear_data, ggplot2::aes(x = x1, y = x2, color = y)) +
ggplot2::geom_point(alpha = 0.7) +
ggplot2::labs(
title = "非线性可分数据",
x = "特征1", y = "特征2"
) +
ggplot2::theme_minimal()
# 线性核
svm_linear <- e1071::svm(y ~ x1 + x2,
data = nonlinear_data,
kernel = "linear",
cost = 1)
# 预测准确率
pred_linear <- predict(svm_linear, nonlinear_data)
acc_linear <- mean(pred_linear == nonlinear_data$y)
print(paste("线性核准确率:", round(acc_linear, 3)))
## [1] "线性核准确率: 0.71"
# 多项式核
svm_poly <- e1071::svm(y ~ x1 + x2,
data = nonlinear_data,
kernel = "polynomial",
degree = 3,
coef0 = 1,
cost = 1)
pred_poly <- predict(svm_poly, nonlinear_data)
acc_poly <- mean(pred_poly == nonlinear_data$y)
print(paste("多项式核准确率:", round(acc_poly, 3)))
## [1] "多项式核准确率: 1"
核技巧的核心思想:在高维空间中进行内积运算,而不需要显式地计算高维映射。
Mercer定理:只要\(K(x, z)\)是对称正定的,就存在映射\(\phi\)使得\(K(x, z) = \phi(x)^T \phi(z)\)。
# 核技巧示意:2D到3D的映射
# 原始数据
x <- c(1, 2, 3)
z <- c(2, 3, 4)
# 多项式核 K(x,z) = (x^T z)^2
# 等价于映射 phi(x) = (x1^2, sqrt(2)*x1*x2, x2^2)
# 直接计算核函数
kernel_direct <- (sum(x * z))^2
print(paste("直接计算核函数:", kernel_direct))
## [1] "直接计算核函数: 400"
# 通过映射计算
phi_x <- c(x[1]^2, sqrt(2) * x[1] * x[2], x[2]^2)
phi_z <- c(z[1]^2, sqrt(2) * z[1] * z[2], z[2]^2)
kernel_mapped <- sum(phi_x * phi_z)
print(paste("通过映射计算:", kernel_mapped))
## [1] "通过映射计算: 64"
核心意义:核技巧让线性方法能够处理非线性问题,是SVM等算法的核心技术。
举例:医学图像分类
为什么需要核技巧?
医学图像数据(如病理切片、MRI图像)在原始特征空间中往往线性不可分。例如,良性和恶性肿瘤在特征空间中可能形成复杂的非线性边界。核技巧通过隐式映射到高维空间,使原本线性不可分的数据变得线性可分,而无需显式计算高维映射。
研究目的是什么?
可能得到什么样的结果?
使用RBF核的SVM可能在医学图像分类任务上达到90%以上的准确率,而线性SVM可能只有70%。核技巧的代价是需要调参(如RBF核的gamma参数)和更多计算资源。
需要注意什么?
核函数的选择和参数调优对结果影响很大。RBF核是最常用的选择,但对于特定问题,其他核函数可能更合适。核方法在大数据集上计算较慢,需要考虑近似方法或使用线性核。
支持向量回归将SVM扩展到回归问题,寻找在\(\epsilon\)-带内的超平面:
\[\min_{w,b} \frac{1}{2}\|w\|^2 + C \sum_{i=1}^{n} (\xi_i + \xi_i^*)\]
约束条件:\(|y_i - w^T x_i - b| \leq \epsilon + \xi_i\)
# SVR示例
set.seed(42)
# 创建回归数据
n <- 100
x <- runif(n, 0, 10)
y <- sin(x) + rnorm(n, 0, 0.2)
svr_data <- data.frame(x = x, y = y)
# 训练SVR模型
svr_model <- e1071::svm(y ~ x,
data = svr_data,
type = "eps-regression",
kernel = "radial",
epsilon = 0.1,
cost = 10)
# 预测
svr_pred <- predict(svr_model, svr_data)
# 可视化
ggplot2::ggplot(svr_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_line(ggplot2::aes(y = svr_pred), color = "steelblue", linewidth = 1) +
ggplot2::labs(
title = "支持向量回归(SVR)",
subtitle = "蓝线为SVR预测曲线",
x = "x", y = "y"
) +
ggplot2::theme_minimal()
# 比较不同epsilon值
svr_small_eps <- e1071::svm(y ~ x, data = svr_data,
type = "eps-regression", kernel = "radial",
epsilon = 0.01, cost = 10)
svr_large_eps <- e1071::svm(y ~ x, data = svr_data,
type = "eps-regression", kernel = "radial",
epsilon = 0.5, cost = 10)
# 比较支持向量数量
sv_compare <- data.frame(
epsilon = c("0.01", "0.1", "0.5"),
支持向量数 = c(length(svr_small_eps$index),
length(svr_model$index),
length(svr_large_eps$index))
)
print(sv_compare)
## epsilon 支持向量数
## 1 0.01 95
## 2 0.1 67
## 3 0.5 10
当数据存在噪声或异常值时,允许部分样本误分类:
\[\min_{w,b,\xi} \frac{1}{2}\|w\|^2 + C \sum_{i=1}^{n} \xi_i\]
C参数:控制误分类惩罚程度。
# C参数影响示例
set.seed(42)
# 创建有噪声的数据
n <- 100
x1 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
x2 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))
noisy_data <- data.frame(x1 = x1, x2 = x2, y = y)
# 不同C值
svm_c01 <- e1071::svm(y ~ x1 + x2, data = noisy_data,
kernel = "radial", cost = 0.1)
svm_c1 <- e1071::svm(y ~ x1 + x2, data = noisy_data,
kernel = "radial", cost = 1)
svm_c100 <- e1071::svm(y ~ x1 + x2, data = noisy_data,
kernel = "radial", cost = 100)
# 比较支持向量数量
c_compare <- data.frame(
C值 = c("0.1", "1", "100"),
支持向量数 = c(length(svm_c01$index), length(svm_c1$index), length(svm_c100$index))
)
print(c_compare)
## C值 支持向量数
## 1 0.1 63
## 2 1 29
## 3 100 17
核主成分分析将PCA扩展到非线性情况,通过核函数在高维空间进行主成分分析。
# KPCA示例
set.seed(42)
# 创建同心圆数据
n <- 200
theta <- runif(n, 0, 2 * pi)
r <- ifelse(sample(c(0, 1), n, replace = TRUE) == 0,
runif(n, 0.3, 0.5), runif(n, 0.8, 1))
x1 <- r * cos(theta)
x2 <- r * sin(theta)
y <- factor(ifelse(r < 0.6, 0, 1))
kpca_data <- data.frame(x1 = x1, x2 = x2, y = y)
# 使用kernlab包进行KPCA
if (requireNamespace("kernlab", quietly = TRUE)) {
kpca_result <- kernlab::kpca(~ x1 + x2,
data = kpca_data,
kernel = "rbfdot",
kpar = list(sigma = 5))
# 提取主成分
pc_scores <- kernlab::pcv(kpca_result)
kpca_data$PC1 <- pc_scores[, 1]
kpca_data$PC2 <- pc_scores[, 2]
# 可视化KPCA结果
ggplot2::ggplot(kpca_data, ggplot2::aes(x = PC1, y = PC2, color = y)) +
ggplot2::geom_point(alpha = 0.7) +
ggplot2::labs(
title = "核主成分分析(KPCA)结果",
x = "第一主成分", y = "第二主成分"
) +
ggplot2::theme_minimal()
}
RBF核的两个关键参数:
# 使用caret进行SVM调参
set.seed(42)
# 定义参数网格
tune_grid <- expand.grid(
sigma = c(0.01, 0.1, 1, 10),
C = c(0.1, 1, 10, 100)
)
# 使用caret调参
svm_tune <- caret::train(
y ~ x1 + x2,
data = nonlinear_data,
method = "svmRadial",
trControl = caret::trainControl(method = "cv", number = 5),
tuneGrid = tune_grid
)
# 查看调参结果
print(svm_tune)
## Support Vector Machines with Radial Basis Function Kernel
##
## 200 samples
## 2 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 161, 160, 159, 160, 160
## Resampling results across tuning parameters:
##
## sigma C Accuracy Kappa
## 0.01 0.1 0.5350094 0.0000000
## 0.01 1.0 0.5350094 0.0000000
## 0.01 10.0 0.9900000 0.9795918
## 0.01 100.0 1.0000000 1.0000000
## 0.10 0.1 0.9008474 0.7952015
## 0.10 1.0 1.0000000 1.0000000
## 0.10 10.0 1.0000000 1.0000000
## 0.10 100.0 1.0000000 1.0000000
## 1.00 0.1 1.0000000 1.0000000
## 1.00 1.0 1.0000000 1.0000000
## 1.00 10.0 1.0000000 1.0000000
## 1.00 100.0 1.0000000 1.0000000
## 10.00 0.1 0.8900969 0.7823357
## 10.00 1.0 1.0000000 1.0000000
## 10.00 10.0 1.0000000 1.0000000
## 10.00 100.0 1.0000000 1.0000000
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 1 and C = 0.1.
# 可视化调参结果
ggplot2::ggplot(svm_tune$results, ggplot2::aes(x = sigma, y = Accuracy, color = factor(C))) +
ggplot2::geom_point(size = 3) +
ggplot2::geom_line() +
ggplot2::labs(
title = "SVM参数调优结果",
x = "sigma(gamma)", y = "准确率", color = "C"
) +
ggplot2::theme_minimal()
K近邻(K-Nearest Neighbors, KNN)是一种基于实例的学习方法:
# KNN分类示例
set.seed(42)
# 创建模拟数据
n <- 200
knn_data <- data.frame(
feature1 = c(rnorm(n/2, mean = 2, sd = 1), rnorm(n/2, mean = 5, sd = 1)),
feature2 = c(rnorm(n/2, mean = 2, sd = 1), rnorm(n/2, mean = 5, sd = 1)),
class = factor(c(rep("A", n/2), rep("B", n/2)))
)
# 划分训练集和测试集
train_idx <- sample(1:n, 0.7 * n)
train_data <- knn_data[train_idx, ]
test_data <- knn_data[-train_idx, ]
# 使用class包进行KNN分类
knn_pred <- class::knn(
train = train_data[, c("feature1", "feature2")],
test = test_data[, c("feature1", "feature2")],
cl = train_data$class,
k = 5
)
# 计算准确率
accuracy <- mean(knn_pred == test_data$class)
print(paste("KNN准确率:", round(accuracy, 3)))
## [1] "KNN准确率: 0.983"
# 混淆矩阵
table(实际 = test_data$class, 预测 = knn_pred)
## 预测
## 实际 A B
## A 30 1
## B 0 29
# 创建网格用于可视化决策边界
grid_x1 <- seq(min(knn_data$feature1) - 1, max(knn_data$feature1) + 1, by = 0.1)
grid_x2 <- seq(min(knn_data$feature2) - 1, max(knn_data$feature2) + 1, by = 0.1)
grid <- expand.grid(feature1 = grid_x1, feature2 = grid_x2)
# 预测网格点的类别
grid_pred <- class::knn(
train = train_data[, c("feature1", "feature2")],
test = grid,
cl = train_data$class,
k = 5
)
grid$class <- grid_pred
# 可视化
ggplot2::ggplot() +
ggplot2::geom_tile(data = grid, ggplot2::aes(x = feature1, y = feature2, fill = class), alpha = 0.3) +
ggplot2::geom_point(data = train_data, ggplot2::aes(x = feature1, y = feature2, color = class), size = 2) +
ggplot2::labs(
title = "KNN决策边界(K=5)",
x = "特征1", y = "特征2"
) +
ggplot2::theme_minimal()
欧氏距离(Euclidean Distance):
\[d(x, y) = \sqrt{\sum_{i=1}^{p}(x_i - y_i)^2}\]
曼哈顿距离(Manhattan Distance):
\[d(x, y) = \sum_{i=1}^{p}|x_i - y_i|\]
闵可夫斯基距离(Minkowski Distance):
\[d(x, y) = \left(\sum_{i=1}^{p}|x_i - y_i|^q\right)^{1/q}\]
当\(q=2\)时为欧氏距离,\(q=1\)时为曼哈顿距离。
# 比较不同距离度量
# 欧氏距离(默认)
knn_euclidean <- class::knn(
train = train_data[, c("feature1", "feature2")],
test = test_data[, c("feature1", "feature2")],
cl = train_data$class,
k = 5
)
# 使用caret比较不同距离度量
knn_manhattan <- caret::knn3(
class ~ feature1 + feature2,
data = train_data,
k = 5
)
# 预测
pred_manhattan <- predict(knn_manhattan, test_data[, c("feature1", "feature2")], type = "class")
# 比较
distance_compare <- data.frame(
距离度量 = c("欧氏距离", "曼哈顿距离"),
准确率 = c(
mean(knn_euclidean == test_data$class),
mean(pred_manhattan == test_data$class)
)
)
print(distance_compare)
## 距离度量 准确率
## 1 欧氏距离 0.9833333
## 2 曼哈顿距离 0.9833333
加权KNN:距离越近的邻居权重越大,减少K值选择的影响。
常用权重函数:
\[w_i = \frac{1}{d(x, x_i)^2}\]
# 加权KNN示例
# 使用kknn包
if (requireNamespace("kknn", quietly = TRUE)) {
library(kknn)
# 训练加权KNN
weighted_knn <- kknn::kknn(
class ~ feature1 + feature2,
train = train_data,
test = test_data,
k = 5,
distance = 2, # 欧氏距离
kernel = "gaussian" # 高斯核权重
)
# 预测
weighted_pred <- fitted(weighted_knn)
# 比较加权与非加权
compare_df <- data.frame(
方法 = c("普通KNN", "加权KNN"),
准确率 = c(
mean(knn_euclidean == test_data$class),
mean(weighted_pred == test_data$class)
)
)
print(compare_df)
}
KNN基于距离计算,如果特征尺度差异大,大尺度特征会主导距离计算。
# 特征缩放示例
set.seed(42)
# 创建尺度差异大的数据
scale_data <- data.frame(
age = round(runif(200, 20, 80)), # 范围:20-80
income = round(runif(200, 20000, 200000)), # 范围:20000-200000
class = factor(sample(c("A", "B"), 200, replace = TRUE))
)
# 不缩放的KNN
train_idx <- sample(1:200, 150)
train_scale <- scale_data[train_idx, ]
test_scale <- scale_data[-train_idx, ]
knn_no_scale <- class::knn(
train = train_scale[, c("age", "income")],
test = test_scale[, c("age", "income")],
cl = train_scale$class,
k = 5
)
# 标准化后
train_scaled <- train_scale
test_scaled <- test_scale
# 标准化(使用训练集的均值和标准差)
age_mean <- mean(train_scaled$age)
age_sd <- sd(train_scaled$age)
income_mean <- mean(train_scaled$income)
income_sd <- sd(train_scaled$income)
train_scaled$age <- (train_scaled$age - age_mean) / age_sd
train_scaled$income <- (train_scaled$income - income_mean) / income_sd
test_scaled$age <- (test_scaled$age - age_mean) / age_sd
test_scaled$income <- (test_scaled$income - income_mean) / income_sd
knn_scaled <- class::knn(
train = train_scaled[, c("age", "income")],
test = test_scaled[, c("age", "income")],
cl = train_scaled$class,
k = 5
)
# 比较
scale_compare <- data.frame(
方法 = c("未缩放", "标准化后"),
准确率 = c(
mean(knn_no_scale == test_scale$class),
mean(knn_scaled == test_scale$class)
)
)
print(scale_compare)
## 方法 准确率
## 1 未缩放 0.48
## 2 标准化后 0.50
# 不同K值比较
set.seed(42)
k_values <- c(1, 3, 5, 7, 10, 15, 20, 30)
accuracies <- numeric(length(k_values))
for (i in seq_along(k_values)) {
knn_pred <- class::knn(
train = train_data[, c("feature1", "feature2")],
test = test_data[, c("feature1", "feature2")],
cl = train_data$class,
k = k_values[i]
)
accuracies[i] <- mean(knn_pred == test_data$class)
}
k_compare <- data.frame(
K值 = k_values,
准确率 = round(accuracies, 3)
)
print(k_compare)
## K值 准确率
## 1 1 0.983
## 2 3 0.983
## 3 5 0.983
## 4 7 0.983
## 5 10 0.983
## 6 15 0.983
## 7 20 0.967
## 8 30 0.983
# 可视化
ggplot2::ggplot(k_compare, ggplot2::aes(x = K值, y = 准确率)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_point(size = 3, color = "steelblue") +
ggplot2::labs(
title = "K值对KNN准确率的影响",
x = "K值", y = "准确率"
) +
ggplot2::theme_minimal()
# 使用caret进行交叉验证选择K
set.seed(42)
tune_control <- caret::trainControl(method = "cv", number = 5)
knn_tune <- caret::train(
class ~ feature1 + feature2,
data = train_data,
method = "knn",
trControl = tune_control,
tuneGrid = data.frame(k = 1:20)
)
print(knn_tune)
## k-Nearest Neighbors
##
## 140 samples
## 2 predictor
## 2 classes: 'A', 'B'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 112, 112, 112, 113, 111
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 1 0.9785532 0.9569958
## 2 0.9714103 0.9427758
## 3 0.9925926 0.9851240
## 4 0.9925926 0.9851240
## 5 0.9925926 0.9851240
## 6 0.9925926 0.9851240
## 7 0.9925926 0.9851240
## 8 0.9856960 0.9712815
## 9 0.9925926 0.9851240
## 10 0.9856960 0.9712815
## 11 0.9856960 0.9712815
## 12 0.9856960 0.9712815
## 13 0.9856960 0.9712815
## 14 0.9925926 0.9851240
## 15 0.9856960 0.9712815
## 16 0.9925926 0.9851240
## 17 0.9925926 0.9851240
## 18 0.9925926 0.9851240
## 19 0.9925926 0.9851240
## 20 0.9925926 0.9851240
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 20.
# 可视化
ggplot2::ggplot(knn_tune$results, ggplot2::aes(x = k, y = Accuracy)) +
ggplot2::geom_line(color = "steelblue", linewidth = 1) +
ggplot2::geom_point(size = 2, color = "steelblue") +
ggplot2::geom_errorbar(ggplot2::aes(ymin = Accuracy - AccuracySD,
ymax = Accuracy + AccuracySD),
width = 0.5) +
ggplot2::labs(
title = "交叉验证选择最优K值",
x = "K值", y = "准确率"
) +
ggplot2::theme_minimal()
对于大规模数据,暴力搜索效率低。K-D树和球树是常用的加速结构。
K-D树:递归地将空间划分为超矩形区域。
球树:使用超球体划分空间,对高维数据更有效。
# K-D树加速示意(概念演示)
# R中可以使用RANN包进行快速最近邻搜索
if (requireNamespace("RANN", quietly = TRUE)) {
library(RANN)
# 使用nn2函数快速查找最近邻
nn_result <- RANN::nn2(
data = train_data[, c("feature1", "feature2")],
query = test_data[, c("feature1", "feature2")],
k = 5
)
# nn_result$nn.idx 包含最近邻索引
# nn_result$nn.dists 包含距离
print("最近邻索引(前5个测试样本):")
print(head(nn_result$nn.idx))
}
## [1] "最近邻索引(前5个测试样本):"
## [,1] [,2] [,3] [,4] [,5]
## [1,] 45 18 57 76 44
## [2,] 23 127 16 52 128
## [3,] 24 34 104 96 71
## [4,] 23 127 16 52 20
## [5,] 57 77 18 45 76
## [6,] 95 106 49 35 100
KNN也可用于回归问题,预测值为K个邻居的平均值或加权平均值。
# KNN回归示例
set.seed(42)
# 创建回归数据
n <- 200
reg_data <- data.frame(
x = runif(n, 0, 10),
y = sin(runif(n, 0, 10)) + rnorm(n, 0, 0.3)
)
# 划分训练集和测试集
train_idx <- sample(1:n, 150)
train_reg <- reg_data[train_idx, ]
test_reg <- reg_data[-train_idx, ]
# 使用caret进行KNN回归
knn_reg <- caret::knnreg(
y ~ x,
data = train_reg,
k = 5
)
# 预测
pred_reg <- predict(knn_reg, test_reg)
# 计算RMSE
rmse <- sqrt(mean((pred_reg - test_reg$y)^2))
print(paste("KNN回归RMSE:", round(rmse, 3)))
## [1] "KNN回归RMSE: 0.896"
# 可视化
test_reg$pred <- pred_reg
ggplot2::ggplot() +
ggplot2::geom_point(data = train_reg, ggplot2::aes(x = x, y = y), alpha = 0.3) +
ggplot2::geom_point(data = test_reg, ggplot2::aes(x = x, y = pred), color = "steelblue", size = 2) +
ggplot2::labs(
title = "KNN回归",
subtitle = "蓝点为预测值",
x = "x", y = "y"
) +
ggplot2::theme_minimal()
局部加权回归(Locally Weighted Regression)在每个预测点附近拟合一个局部回归模型。
# 局部加权回归示例
set.seed(42)
# 创建数据
n <- 100
lowess_data <- data.frame(
x = runif(n, 0, 10),
y = sin(runif(n, 0, 10)) + rnorm(n, 0, 0.2)
)
# 使用loess进行局部加权回归
loess_model <- loess(y ~ x, data = lowess_data, span = 0.5)
# 预测
x_sorted <- sort(lowess_data$x)
y_pred <- predict(loess_model, newdata = data.frame(x = x_sorted))
# 可视化
ggplot2::ggplot(lowess_data, ggplot2::aes(x = x, y = y)) +
ggplot2::geom_point(alpha = 0.6) +
ggplot2::geom_line(data = data.frame(x = x_sorted, y = y_pred),
ggplot2::aes(x = x, y = y), color = "steelblue", linewidth = 1) +
ggplot2::labs(
title = "局部加权回归(LOESS)",
subtitle = "蓝线为LOESS拟合曲线",
x = "x", y = "y"
) +
ggplot2::theme_minimal()
# 比较不同span值
span_values <- c(0.2, 0.5, 0.8)
colors <- c("red", "steelblue", "green")
plot_data <- data.frame()
for (i in seq_along(span_values)) {
model <- loess(y ~ x, data = lowess_data, span = span_values[i])
pred <- predict(model, newdata = data.frame(x = x_sorted))
plot_data <- rbind(plot_data, data.frame(
x = x_sorted,
y = pred,
span = paste("span =", span_values[i])
))
}
ggplot2::ggplot() +
ggplot2::geom_point(data = lowess_data, ggplot2::aes(x = x, y = y), alpha = 0.3) +
ggplot2::geom_line(data = plot_data, ggplot2::aes(x = x, y = y, color = span), linewidth = 1) +
ggplot2::labs(
title = "不同span参数的LOESS拟合",
x = "x", y = "y", color = "参数"
) +
ggplot2::theme_minimal()
贝叶斯定理描述了条件概率之间的关系:
\[P(A|B) = \frac{P(B|A) \cdot P(A)}{P(B)}\]
在分类问题中:
\[P(y|x) = \frac{P(x|y) \cdot P(y)}{P(x)}\]
朴素贝叶斯假设特征之间相互独立:
\[P(x|y) = \prod_{j=1}^{p} P(x_j|y)\]
分类决策:
\[\hat{y} = \arg\max_y P(y) \prod_{j=1}^{p} P(x_j|y)\]
| 模型 | 特征类型 | 分布假设 |
|---|---|---|
| 高斯朴素贝叶斯 | 连续特征 | 高斯分布 |
| 多项式朴素贝叶斯 | 计数特征 | 多项式分布 |
| 伯努利朴素贝叶斯 | 二值特征 | 伯努利分布 |
# 高斯朴素贝叶斯示例
set.seed(42)
# 使用iris数据集
data(iris)
# 划分训练集和测试集
train_idx <- sample(1:nrow(iris), 0.7 * nrow(iris))
train_iris <- iris[train_idx, ]
test_iris <- iris[-train_idx, ]
# 训练高斯朴素贝叶斯
nb_gaussian <- e1071::naiveBayes(Species ~ ., data = train_iris)
# 预测
nb_pred <- predict(nb_gaussian, test_iris)
# 混淆矩阵
conf_matrix <- table(实际 = test_iris$Species, 预测 = nb_pred)
print(conf_matrix)
## 预测
## 实际 setosa versicolor virginica
## setosa 12 0 0
## versicolor 0 14 1
## virginica 0 1 17
# 准确率
accuracy <- mean(nb_pred == test_iris$Species)
print(paste("高斯朴素贝叶斯准确率:", round(accuracy, 3)))
## [1] "高斯朴素贝叶斯准确率: 0.956"
# 查看各类别的先验概率
print("先验概率:")
## [1] "先验概率:"
print(nb_gaussian$apriori)
## Y
## setosa versicolor virginica
## 38 35 32
# 查看各特征的条件概率参数(均值和标准差)
print("条件概率参数:")
## [1] "条件概率参数:"
print(nb_gaussian$tables)
## $Sepal.Length
## Sepal.Length
## Y [,1] [,2]
## setosa 5.015789 0.3405293
## versicolor 5.874286 0.5135476
## virginica 6.606250 0.6983564
##
## $Sepal.Width
## Sepal.Width
## Y [,1] [,2]
## setosa 3.426316 0.4084805
## versicolor 2.745714 0.3229850
## virginica 2.996875 0.3560395
##
## $Petal.Length
## Petal.Length
## Y [,1] [,2]
## setosa 1.447368 0.1389864
## versicolor 4.217143 0.5176222
## virginica 5.543750 0.5529612
##
## $Petal.Width
## Petal.Width
## Y [,1] [,2]
## setosa 0.2473684 0.1156342
## versicolor 1.3228571 0.2001260
## virginica 2.0125000 0.2836968
# 多项式朴素贝叶斯示例:文本分类
set.seed(42)
# 创建模拟文本数据
n <- 200
text_data <- data.frame(
word1 = rpois(n, lambda = sample(c(2, 5, 1), n, replace = TRUE)),
word2 = rpois(n, lambda = sample(c(3, 1, 4), n, replace = TRUE)),
word3 = rpois(n, lambda = sample(c(1, 3, 2), n, replace = TRUE)),
category = factor(sample(c("医疗", "体育", "科技"), n, replace = TRUE))
)
# 划分训练集和测试集
train_idx <- sample(1:n, 150)
train_text <- text_data[train_idx, ]
test_text <- text_data[-train_idx, ]
# 训练多项式朴素贝叶斯
nb_multinomial <- e1071::naiveBayes(category ~ ., data = train_text, type = "raw")
# 预测
pred_text <- predict(nb_multinomial, test_text)
# 准确率
accuracy <- mean(pred_text == test_text$category)
print(paste("多项式朴素贝叶斯准确率:", round(accuracy, 3)))
## [1] "多项式朴素贝叶斯准确率: 0.34"
# 伯努利朴素贝叶斯示例
set.seed(42)
# 创建二值特征数据
n <- 200
binary_data <- data.frame(
symptom1 = sample(c(0, 1), n, replace = TRUE),
symptom2 = sample(c(0, 1), n, replace = TRUE),
symptom3 = sample(c(0, 1), n, replace = TRUE),
symptom4 = sample(c(0, 1), n, replace = TRUE),
disease = factor(sample(c("健康", "患病"), n, replace = TRUE, prob = c(0.7, 0.3)))
)
# 划分数据
train_idx <- sample(1:n, 150)
train_binary <- binary_data[train_idx, ]
test_binary <- binary_data[-train_idx, ]
# 训练伯努利朴素贝叶斯
nb_bernoulli <- e1071::naiveBayes(disease ~ ., data = train_binary)
# 预测
pred_binary <- predict(nb_bernoulli, test_binary)
# 准确率
accuracy <- mean(pred_binary == test_binary$disease)
print(paste("伯努利朴素贝叶斯准确率:", round(accuracy, 3)))
## [1] "伯努利朴素贝叶斯准确率: 0.7"
高斯判别分析假设各类别的数据服从多元高斯分布:
\[P(x|y=k) = \frac{1}{(2\pi)^{p/2}|\Sigma_k|^{1/2}} \exp\left(-\frac{1}{2}(x-\mu_k)^T \Sigma_k^{-1} (x-\mu_k)\right)\]
# GDA示例(使用MASS包的lda函数)
set.seed(42)
# 使用iris数据
gda_model <- MASS::lda(Species ~ ., data = train_iris)
# 查看模型
print(gda_model)
## Call:
## lda(Species ~ ., data = train_iris)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3619048 0.3333333 0.3047619
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.015789 3.426316 1.447368 0.2473684
## versicolor 5.874286 2.745714 4.217143 1.3228571
## virginica 6.606250 2.996875 5.543750 2.0125000
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 1.084659 -0.1357075
## Sepal.Width 1.221953 -2.0950103
## Petal.Length -2.431518 0.8233400
## Petal.Width -2.643108 -2.4996356
##
## Proportion of trace:
## LD1 LD2
## 0.9916 0.0084
# 预测
gda_pred <- predict(gda_model, test_iris)
# 混淆矩阵
table(实际 = test_iris$Species, 预测 = gda_pred$class)
## 预测
## 实际 setosa versicolor virginica
## setosa 12 0 0
## versicolor 0 15 0
## virginica 0 1 17
# 准确率
accuracy <- mean(gda_pred$class == test_iris$Species)
print(paste("GDA准确率:", round(accuracy, 3)))
## [1] "GDA准确率: 0.978"
贝叶斯网络是一种概率图模型,用有向无环图(DAG)表示变量之间的依赖关系。
组成: - 节点:随机变量 - 边:依赖关系 - 条件概率表(CPT):每个节点的条件概率分布
# 贝叶斯网络示例(概念演示)
# 使用bnlearn包(如果可用)
if (requireNamespace("bnlearn", quietly = TRUE)) {
library(bnlearn)
# 创建一个简单的贝叶斯网络结构
# 假设:吸烟 -> 肺癌 -> 死亡;吸烟 -> 心脏病 -> 死亡
# 定义网络结构
dag <- model2network("[吸烟][肺癌|吸烟][心脏病|吸烟][死亡|肺癌:心脏病]")
# 可视化网络结构
graphviz.chart(dag, main = "贝叶斯网络结构示例")
} else {
# 手动绘制贝叶斯网络概念图
bn_nodes <- data.frame(
node = c("吸烟", "肺癌", "心脏病", "死亡"),
x = c(1, 0.5, 1.5, 1),
y = c(3, 2, 2, 1)
)
ggplot2::ggplot() +
ggplot2::geom_segment(ggplot2::aes(x = 1, xend = 0.5, y = 2.8, yend = 2.3),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
ggplot2::geom_segment(ggplot2::aes(x = 1, xend = 1.5, y = 2.8, yend = 2.3),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
ggplot2::geom_segment(ggplot2::aes(x = 0.5, xend = 1, y = 1.7, yend = 1.3),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
ggplot2::geom_segment(ggplot2::aes(x = 1.5, xend = 1, y = 1.7, yend = 1.3),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
ggplot2::geom_point(data = bn_nodes, ggplot2::aes(x = x, y = y), size = 12, color = "steelblue") +
ggplot2::geom_text(data = bn_nodes, ggplot2::aes(x = x, y = y, label = node), vjust = -2, size = 4) +
ggplot2::labs(
title = "贝叶斯网络结构示例",
subtitle = "吸烟 -> 肺癌 -> 死亡;吸烟 -> 心脏病 -> 死亡"
) +
ggplot2::theme_void()
}
隐马尔可夫模型用于建模序列数据,包含:
# HMM概念示意
# 创建一个简单的疾病状态转移模型
# 状态:健康(H)、疾病(S)
# 观测:正常检查结果(N)、异常检查结果(A)
# 转移概率矩阵
transition_matrix <- matrix(
c(0.9, 0.1, # H -> H, H -> S
0.3, 0.7), # S -> H, S -> S
nrow = 2, byrow = TRUE
)
rownames(transition_matrix) <- c("健康", "疾病")
colnames(transition_matrix) <- c("健康", "疾病")
print("转移概率矩阵:")
## [1] "转移概率矩阵:"
print(transition_matrix)
## 健康 疾病
## 健康 0.9 0.1
## 疾病 0.3 0.7
# 发射概率矩阵
emission_matrix <- matrix(
c(0.8, 0.2, # H -> N, H -> A
0.3, 0.7), # S -> N, S -> A
nrow = 2, byrow = TRUE
)
rownames(emission_matrix) <- c("健康", "疾病")
colnames(emission_matrix) <- c("正常", "异常")
print("发射概率矩阵:")
## [1] "发射概率矩阵:"
print(emission_matrix)
## 正常 异常
## 健康 0.8 0.2
## 疾病 0.3 0.7
# 初始概率
initial_prob <- c(健康 = 0.9, 疾病 = 0.1)
print("初始概率:")
## [1] "初始概率:"
print(initial_prob)
## 健康 疾病
## 0.9 0.1
# 可视化HMM结构
hmm_states <- data.frame(
state = c("健康", "疾病"),
x = c(1, 3),
y = c(2, 2)
)
hmm_obs <- data.frame(
obs = c("正常", "异常"),
x = c(1, 3),
y = c(1, 1)
)
ggplot2::ggplot() +
ggplot2::geom_segment(ggplot2::aes(x = 1.3, xend = 2.7, y = 2, yend = 2),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
ggplot2::geom_segment(ggplot2::aes(x = 2.7, xend = 1.3, y = 1.85, yend = 1.85),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
ggplot2::geom_segment(ggplot2::aes(x = 1, xend = 1, y = 1.7, yend = 1.3),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "steelblue", linetype = "dashed") +
ggplot2::geom_segment(ggplot2::aes(x = 3, xend = 3, y = 1.7, yend = 1.3),
arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "steelblue", linetype = "dashed") +
ggplot2::geom_point(data = hmm_states, ggplot2::aes(x = x, y = y),
size = 12, color = "steelblue", shape = 15) +
ggplot2::geom_point(data = hmm_obs, ggplot2::aes(x = x, y = y),
size = 10, color = "salmon", shape = 16) +
ggplot2::geom_text(data = hmm_states, ggplot2::aes(x = x, y = y, label = state),
vjust = -2, size = 4) +
ggplot2::geom_text(data = hmm_obs, ggplot2::aes(x = x, y = y, label = obs),
vjust = -2, size = 4) +
ggplot2::annotate("text", x = 2, y = 2.15, label = "转移", size = 3, color = "gray30") +
ggplot2::annotate("text", x = 1.5, y = 1.5, label = "发射", size = 3, color = "steelblue") +
ggplot2::annotate("text", x = 3.5, y = 1.5, label = "发射", size = 3, color = "steelblue") +
ggplot2::labs(
title = "隐马尔可夫模型结构",
subtitle = "方形:隐藏状态;圆形:观测值;实线:转移;虚线:发射"
) +
ggplot2::theme_void()
高斯过程(Gaussian Process, GP)是一种非参数贝叶斯方法,定义了函数的分布。
定义:对于任意输入点集合\(\{x_1, ..., x_n\}\),函数值\(\{f(x_1), ..., f(x_n)\}\)服从多元高斯分布。
核心要素: - 均值函数 \(m(x)\) - 协方差函数(核函数) \(k(x, x')\)
# 高斯过程回归示例
set.seed(42)
# 创建数据
n_train <- 20
x_train <- sort(runif(n_train, 0, 10))
y_train <- sin(x_train) + rnorm(n_train, 0, 0.1)
# 测试点
x_test <- seq(0, 10, length.out = 100)
# 定义RBF核函数
rbf_kernel <- function(x1, x2, sigma = 1, length_scale = 1) {
sigma^2 * exp(-0.5 * (outer(x1, x2, "-") / length_scale)^2)
}
# 计算协方差矩阵
K <- rbf_kernel(x_train, x_train)
K_star <- rbf_kernel(x_test, x_train)
K_star_star <- rbf_kernel(x_test, x_test)
# 添加噪声
sigma_n <- 0.1
K_noise <- K + sigma_n^2 * diag(n_train)
# 高斯过程预测
L <- chol(K_noise)
alpha <- solve(t(L), solve(L, y_train))
y_pred <- K_star %*% alpha
# 计算预测方差
v <- solve(L, t(K_star))
y_var <- diag(K_star_star) - colSums(v^2)
y_sd <- sqrt(pmax(y_var, 0))
# 可视化
gp_data <- data.frame(
x = x_test,
y_mean = as.vector(y_pred),
y_lower = as.vector(y_pred - 1.96 * y_sd),
y_upper = as.vector(y_pred + 1.96 * y_sd)
)
ggplot2::ggplot() +
ggplot2::geom_ribbon(data = gp_data, ggplot2::aes(x = x, ymin = y_lower, ymax = y_upper),
fill = "steelblue", alpha = 0.3) +
ggplot2::geom_line(data = gp_data, ggplot2::aes(x = x, y = y_mean), color = "steelblue", linewidth = 1) +
ggplot2::geom_point(data = data.frame(x = x_train, y = y_train),
ggplot2::aes(x = x, y = y), size = 2) +
ggplot2::labs(
title = "高斯过程回归",
subtitle = "阴影区域为95%置信区间",
x = "x", y = "y"
) +
ggplot2::theme_minimal()
# 比较不同长度尺度
length_scales <- c(0.5, 1, 2)
gp_compare <- data.frame()
for (l in length_scales) {
K <- rbf_kernel(x_train, x_train, length_scale = l)
K_star <- rbf_kernel(x_test, x_train, length_scale = l)
K_star_star <- rbf_kernel(x_test, x_test, length_scale = l)
K_noise <- K + sigma_n^2 * diag(n_train)
L <- chol(K_noise)
alpha <- solve(t(L), solve(L, y_train))
y_pred <- K_star %*% alpha
gp_compare <- rbind(gp_compare, data.frame(
x = x_test,
y = as.vector(y_pred),
length_scale = paste("l =", l)
))
}
ggplot2::ggplot() +
ggplot2::geom_point(data = data.frame(x = x_train, y = y_train),
ggplot2::aes(x = x, y = y), size = 2) +
ggplot2::geom_line(data = gp_compare, ggplot2::aes(x = x, y = y, color = length_scale), linewidth = 1) +
ggplot2::labs(
title = "不同长度尺度的高斯过程",
x = "x", y = "y", color = "长度尺度"
) +
ggplot2::theme_minimal()