如何在 R 中绘制散点包络曲线，并使用 0.99 分位数回归？

2026-11-04数据小陈

可以用 0.99 分位数回归估计散点的上部边界，也可以按照给出的函数形式做非线性分位数拟合。分位数曲线允许部分观测点落在上方；如果要求包住所有观测点，则需要采用带包络约束的方法。

当前材料没有附上手绘红线，无法确定它的具体形状。下面分别说明“估计上边界”和“包住全部观测点”这两个目标。

先确定包络线需要满足什么条件

上边界常见的目标有两种：

条件高分位数曲线：例如估计给定 x 时 y 的第 99 百分位数，允许少量观测值高于曲线。
观测数据的严格上包络线：要求每个观测点都满足 y_i <= f(x_i)。

0.99 分位数回归对应第一种目标。它不保证训练数据恰好有 99% 的点位于曲线下方，也不保证曲线高于所有点。

你给出的函数为：

[ f(x)=a\exp(bx)\left[1-\left(\frac{x-z}{w}\right)^2\right]. ]

这个函数规定了可拟合的形状。数据的上边界是否符合这种形状，需要结合拟合结果判断。普通最小二乘拟合主要描述均值趋势，不会自动得到上包络线。

用指定函数做 0.99 分位数拟合

分位数回归最小化非对称绝对误差。令残差为 r = y - f(x)，则损失为：

[ \rho_\tau(r)=\begin{cases} \tau r,&r\geq0,\ (\tau-1)r,&r<0. \end{cases} ]

当 tau = 0.99 时，对于大小相同的误差，观测点高于曲线时的惩罚是低于曲线时的 99 倍。这会推动拟合曲线靠近数据上部。

下面用基础 R 的 optim() 进行数值优化。先运行问题中的 points 定义，再运行以下代码。

为减少指数项的数值问题，先把 x 缩放到 [-1, 1]。代码还增加了一个建模条件：曲线在观测区间内保持正值。这个条件适合当前给出的正值数据，但并非原函数必然具有的性质。

library(ggplot2)

tau <- 0.99
x_mid <- mean(range(points$x))
x_scale <- diff(range(points$x)) / 2

# Parameters: log(A), B, Z, log(delta)
# Model in scaled coordinates:
# A * exp(B * u) * (1 - ((u - Z) / W)^2)
predict_env <- function(par, x) {
  u <- (x - x_mid) / x_scale
  A <- exp(par[1])
  B <- par[2]
  Z <- par[3]
  W <- 1 + abs(Z) + exp(par[4])

  A * exp(B * u) * (1 - ((u - Z) / W)^2)
}

quantile_loss <- function(par) {
  fitted <- predict_env(par, points$x)
  if (any(!is.finite(fitted))) return(1e100)

  r <- points$y - fitted
  sum(ifelse(r >= 0, tau * r, (tau - 1) * r))
}

# Try several starting values for the exponential slope.
fits <- lapply(c(-1, 0, 1), function(B_start) {
  optim(
    par = c(log(max(points$y)), B_start, 0, log(0.5)),
    fn = quantile_loss,
    method = "Nelder-Mead",
    control = list(maxit = 10000, reltol = 1e-9)
  )
})

losses <- vapply(fits, function(fit) fit$value, numeric(1))
fit <- fits[[which.min(losses)]]

if (fit$convergence != 0) {
  warning("Optimization did not report convergence; inspect the fit.")
}

pred <- data.frame(
  x = seq(min(points$x), max(points$x), length.out = 400)
)
pred$y <- predict_env(fit$par, pred$x)

ggplot(points, aes(x, y)) +
  geom_point() +
  geom_line(data = pred, colour = "red", linewidth = 1) +
  theme_bw()

这里的 W > 1 + abs(Z) 保证缩放区间内 abs(u - Z) < W，因此括号中的二次项保持为正。

这仍属于原先指定的函数族，只是换了一套便于优化的参数。如果需要换回原始的 a、b、z、w，可以运行：

A <- exp(fit$par[1])
B <- fit$par[2]
Z <- fit$par[3]
W <- 1 + abs(Z) + exp(fit$par[4])

coef_original <- c(
  a = A * exp(-B * x_mid / x_scale),
  b = B / x_scale,
  z = x_mid + x_scale * Z,
  w = x_scale * W
)

coef_original

数值优化可能对初值敏感。尝试多个初值有助于检查稳定性，但不能证明找到了全局最优解；convergence == 0 也不代表模型形状一定合适。上述代码没有预设或声称任何拟合参数结果。

如果必须包住全部观测点

这组数据存在重复的 x。一个直接的方法是取每个 x 对应的最大 y，再按 x 排序连接：

upper <- aggregate(y ~ x, data = points, FUN = max)
upper <- upper[order(upper$x), ]

ggplot(points, aes(x, y)) +
  geom_point() +
  geom_line(data = upper, colour = "red", linewidth = 1) +
  theme_bw()

在每个已观测的 x 位置，这条折线都不低于对应的所有观测值。这是一种简单的样本上边界构造方法，不一定平滑，也不能保证包住未来数据。

如果要求“平滑、采用指定函数、且包住全部点”，就需要在拟合时显式加入约束：

[ f(x_i;a,b,z,w)\geq y_i,\quad\forall i. ]

还需要选择一个目标，避免曲线被任意抬高，例如最小化曲线与观测值之间的总间距。指定函数也可能无法同时满足期望的形状和全部约束，因此需要单独求解这个约束优化问题。

选择 0.99 分位数时要注意什么

这组数据在每个 x 位置的观测数量有限。极端分位数会依赖少数较大的观测值，也会明显受到函数形状的影响。可以分别拟合 tau = 0.95 和 tau = 0.99，观察曲线是否发生很大变化；比较时不要只凭曲线是否贴近最高点判断模型好坏。

原函数在 abs(x - z) > abs(w) 时会产生负值。上面的正值限制只针对观测区间，绘图也限制在该区间内，不宜直接据此外推。

建议将自定义函数命名为 envelope_fun 等名称，避免覆盖 R 自带的 curve()。