如何平滑数据并强制单调性

机器算法验证 回归 平滑
2022-02-15 04:56:52

我有一些数据要平滑,以便平滑点单调递减。我的数据急剧下降,然后开始趋于平稳。这是一个使用 R 的示例

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))
ggplot(df, aes(x=x, y=y))+geom_line()

要平滑的数据图

我可以使用什么好的平滑技术?此外,如果我可以强制第一个平滑点接近我观察到的点,那就太好了。

2个回答

您可以通过mgcv包中的mono.con()pcls()函数使用具有单调性约束的惩罚样条来执行此操作。由于这些功能不如 用户友好,因此需要进行一些操作,但步骤如下所示,主要基于来自 的示例,已修改以适合您提供的示例数据:gam()?pcls

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))

## Set up the size of the basis functions/number of knots
k <- 5
## This fits the unconstrained model but gets us smoothness parameters that
## that we will need later
unc <- gam(y ~ s(x, k = k, bs = "cr"), data = df)

## This creates the cubic spline basis functions of `x`
## It returns an object containing the penalty matrix for the spline
## among other things; see ?smooth.construct for description of each
## element in the returned object
sm <- smoothCon(s(x, k = k, bs = "cr"), df, knots = NULL)[[1]]

## This gets the constraint matrix and constraint vector that imposes
## linear constraints to enforce montonicity on a cubic regression spline
## the key thing you need to change is `up`.
## `up = TRUE` == increasing function
## `up = FALSE` == decreasing function (as per your example)
## `xp` is a vector of knot locations that we get back from smoothCon
F <- mono.con(sm$xp, up = FALSE)   # get constraints: up = FALSE == Decreasing constraint!

现在我们需要填写传递给pcls()包含我们想要拟合的惩罚约束模型细节的对象

## Fill in G, the object pcsl needs to fit; this is just what `pcls` says it needs:
## X is the model matrix (of the basis functions)
## C is the identifiability constraints - no constraints needed here
##   for the single smooth
## sp are the smoothness parameters from the unconstrained GAM
## p/xp are the knot locations again, but negated for a decreasing function
## y is the response data
## w are weights and this is fancy code for a vector of 1s of length(y)
G <- list(X = sm$X, C = matrix(0,0,0), sp = unc$sp,
          p = -sm$xp, # note the - here! This is for decreasing fits!
      y = df$y,
          w = df$y*0+1)
G$Ain <- F$A    # the monotonicity constraint matrix
G$bin <- F$b    # the monotonicity constraint vector, both from mono.con
G$S <- sm$S     # the penalty matrix for the cubic spline
G$off <- 0      # location of offsets in the penalty matrix

现在我们终于可以进行拟合了

## Do the constrained fit 
p <- pcls(G)  # fit spline (using s.p. from unconstrained fit)

p包含对应于样条的基函数的系数向量。为了可视化拟合的样条曲线,我们可以从模型中预测 x 范围内的 100 个位置。我们做了 100 个值,以便在绘图上获得一条漂亮的平滑线。

## predict at 100 locations over range of x - get a smooth line on the plot
newx <- with(df, data.frame(x = seq(min(x), max(x), length = 100)))

为了生成预测值,我们使用Predict.matrix(),它生成一个矩阵,使得当乘以系数p从拟合模型中产生预测值时:

fv <- Predict.matrix(sm, newx) %*% p
newx <- transform(newx, yhat = fv[,1])

plot(y ~ x, data = df, pch = 16)
lines(yhat ~ x, data = newx, col = "red")

这会产生:

在此处输入图像描述

我会留给你把数据整理成一个整洁的表格,以便用ggplot绘图......

您可以通过增加 的基函数的维度来强制进行更紧密的拟合(部分回答您关于让第一个数据点更平滑拟合的问题)x例如,设置k等于8( k <- 8) 并重新运行上面的代码,我们得到

在此处输入图像描述

你不能把k这些数据推得更高,你必须小心过度拟合;所做pcls()的只是在给定约束和提供的基函数的情况下解决惩罚最小二乘问题,它没有为您执行平滑度选择 - 我不知道......)

如果您想要插值,请查看?splinefun具有单调约束的 Hermite 样条和三次样条的基本 R 函数。在这种情况下,您不能使用它,因为数据不是严格单调的。

Natalya Pya 最近的骗局包基于Pya & Wood (2015) 的论文“Shape constrainedadditive models”可以使 Gavin 的优秀答案中提到的部分过程变得更加容易。

library(scam)
con <- scam(y ~ s(x, k = k, bs = "mpd"), data = df)
plot(con)

您可以使用许多 bs 函数 - 在上面我使用 mpd 表示“单调递减 P 样条曲线”,但它也具有单独或与单调约束一起强制凸度或凹度的函数。