如何在R中绘制包络曲线？含自定义函数与分位数回归咨询

一、用自定义非线性函数拟合包络曲线 #

你给出的自定义函数是含4个参数的非线性模型，需要先通过非线性最小二乘法拟合参数，再生成曲线数据，最后添加手绘风格的红色线条。

步骤1：拟合非线性模型 #

加载基础包，用nls函数拟合模型，注意要给参数设置合理初始值（可根据数据分布大致估计）：

library(ggplot2)

# 定义自定义曲线函数
custom_curve <- function(x, a, b, z, w) {
  a * exp(b * x) * (1 - ((x - z)/w)^2)
}

# 拟合模型，初始值可根据数据调整
fit <- nls(y ~ custom_curve(x, a, b, z, w), 
           data = points,
           start = list(a = 0.8, b = -0.1, z = 18, w = 10))

# 查看拟合参数结果
summary(fit)

步骤2：生成曲线并绘制手绘风格线条 #

生成覆盖x全范围的序列，预测对应y值，再用两种方式实现手绘效果：

# 生成曲线数据
x_seq <- seq(min(points$x), max(points$x), length.out = 100)
curve_data <- data.frame(
  x = x_seq,
  y = predict(fit, newdata = list(x = x_seq))
)

# 方法1：用ggrough包实现原生手绘风格
install.packages("ggrough")
library(ggrough)

ggplot(points, aes(x, y)) +
  geom_point() +
  geom_rough_line(data = curve_data, aes(x, y), color = "red", linewidth = 1) +
  theme_bw()

# 方法2：不额外装包，用抖动模拟手绘感
ggplot(points, aes(x, y)) +
  geom_point() +
  geom_line(data = curve_data, aes(x, y + rnorm(nrow(curve_data), 0, 0.02)), 
            color = "red", linewidth = 1, alpha = 0.8) +
  theme_bw()

二、用99分位数回归绘制包络曲线 #

99分位数回归可以精准捕捉数据的上边界，非常适合做包络曲线，具体步骤如下：

步骤1：拟合99分位数回归模型 #

使用quantreg包的rq函数拟合分位数为0.99的模型，可选择线性或非线性样条模型：

install.packages("quantreg")
library(quantreg)
library(splines) # 用于非线性样条拟合

# 线性99分位数回归
fit_qr <- rq(y ~ x, data = points, tau = 0.99)

# 非线性样条99分位数回归（更贴合数据非线性趋势）
fit_qr_spline <- rq(y ~ bs(x, df = 5), data = points, tau = 0.99)

步骤2：绘制手绘风格的分位数曲线 #

生成预测数据后，添加手绘风格线条：

# 生成样条模型的预测数据
qr_curve_data <- data.frame(
  x = x_seq,
  y = predict(fit_qr_spline, newdata = list(x = x_seq))
)

# 手绘风格绘制
ggplot(points, aes(x, y)) +
  geom_point() +
  geom_rough_line(data = qr_curve_data, aes(x, y), color = "red", linewidth = 1) +
  theme_bw()

# 抖动模拟手绘感（无额外包）
ggplot(points, aes(x, y)) +
  geom_point() +
  geom_line(data = qr_curve_data, aes(x, y + rnorm(nrow(qr_curve_data), 0, 0.02)), 
            color = "red", linewidth = 1, alpha = 0.8) +
  theme_bw()

关键注意点 #

• 自定义非线性模型的初始值直接影响拟合成功率，若报错可根据数据峰值位置调整z、峰值大小调整a；
• 99分位数回归用线性模型可能无法贴合非线性趋势，优先尝试样条模型；
• 手绘抖动幅度需匹配y值范围，本次数据y在0-1之间，0.02的抖动幅度较为合适。

内容的提问来源于stack exchange，提问作者Jiale