要求样条曲线是凸的

Question

我需要对一组数据进行样条插值，并要求得到的函数单调递减和凸性。我将这组数据传递给splinefun，可以保证这些属性，但无法保证得到的函数是凸的。有没有办法对一组数据进行样条插值并要求得到的函数是凸的呢？

- ruser

---

嗨，我想在SO的子站点之一，如Gamedev、MathOverflow或Cross Validated上问这个问题。抱歉，我是Matlab的人，现在是睡觉时间。只是一个提示。 - Stígandr

你能分享一些示例代码、数据或期望的结果来简化回答这个问题吗？ - thelatemail

顺便问一下，你看过这个吗：http://www.stat.colostate.edu/~meyer/penspl.htm？ - thelatemail

你尝试过使用splinefun()函数，方法为"monoH.FC"吗？ - fang

Answer 1

首先提供一些示例数据：

x = c(0,1,2,3,4,5,6)
y = c(2,1, 0.59, 0.27, 0.25, -0.23, -0.45)
dat <- data.frame(x=x,y=y)

单调样条

我们可以使用splinefun(x,y,"monoH.FC")来做一条单调的样条，正如 @fang 建议的那样。

# Setting up Monotonic Spline
MonoSpline = splinefun(x,y,"monoH.FC")
#Getting Ready for plotting Monotonic Spline
xArray = seq(0,6,0.01)
MonoResult = MonoSpline(xArray)

单调凸样条（使用scam包）

要创建一个单调凸样条，您需要使用scam包。接下来我们可以：

# Setting up Monotonic Convex Spline
# install.packages("scam")
require(scam)
MonoConvexSpline <- scam(y~s(x,k=4,bs="mdcx",m=1),data=dat)
MonoConvexSplinePredict =function(Test){
  predict.scam(MonoConvexSpline,data.frame(x = Test))
} 
#Getting Ready for plotting Monotonic Convex Spline
MonoConvexSplineResult = MonoConvexSplinePredict(xArray)

请注意以下内容：

1. 选项bs="mdcx"表示我们需要一个凸降样条。如果您需要增加凸性，降低凹性等，则可以在此处查找适当的bs值。
2. 如果您将非单调数据放入splinefun(x,y,"monoH.FC")函数中，则会收到错误提示。
3. 如果您将非凸数据放入scam函数中，则仍将得到样条。数据会被更改，使小凸度变为凸度。没有任何警告提示，因此请小心，因为您的数据可能看起来完全不同。例如，下面的图是使用与上面相同的代码制作的，只是我们使用了bs="mdcv"用于减少凹函数：

单调凸样条（使用cobs包）

# Convex Cobs Spline
library(cobs)
spCobs = cobs(x , y, constraint = c("decrease", "convex"), nknots = 8)
spCobsResults = predict(spCobs, xArray)[,2] 

绘制一切

然后使用以下方式进行绘制

Plot = qplot(xlab = "x", ylab = "y")
Plot = Plot + geom_line(aes(xArray,MonoResult            , colour = "Monotonic Spline"            ))
Plot = Plot + geom_line(aes(xArray,MonoConvexSplineResult, colour = "Monotonic Convex scam Spline"))
Plot = Plot + geom_line(aes(xArray,spCobsResults         , colour = "Monotonic Convex cobs Spline"))
Plot

速度

• 使用scam函数预测点比使用单调样条或cobs样条需要更长时间。您可以从下面的微型基准测试中看到这一点。
• 但是，相对于普通单调样条，cobs样条和scam样条最初计算需要更长时间。

.

# Prediction
library(microbenchmark)
    microbenchmark(
      MonoSpline(xArray),
      predict.scam(MonoConvexSpline,data.frame(x = xArray)),
      predict(spCobs, xArray)[,2] 
    )

Unit: microseconds
                                                   expr      min        lq      mean    median        uq      max neval
                                     MonoSpline(xArray)  141.540  147.8175  223.3695  156.9490  167.9830 1593.456   100
 predict.scam(MonoConvexSpline, data.frame(x = xArray)) 2778.655 2838.0095 3161.2282 2914.8665 3153.4285 6168.741   100
                           predict(spCobs, xArray)[, 2]  125.179  133.1690  155.1226  145.1535  162.2755  366.784   100

# Calculating Spline
library(microbenchmark)
microbenchmark(
  splinefun(x,y,"monoH.FC"),
  scam(y~s(x,k=4,bs="mdcx",m=1),data=dat),
  cobs(x , y, constraint = c("decrease", "convex"), nknots = 8) 
)

Unit: microseconds
                                                         expr        min         lq        mean      median         uq       max neval
                                  splinefun(x, y, "monoH.FC")     90.175    127.462    411.6407    153.7155    198.993  24877.47   100
        scam(y ~ s(x, k = 4, bs = "mdcx", m = 1), data = dat) 166769.270 196719.139 231631.5321 224372.7940 265074.525 355734.37   100
 cobs(x, y, constraint = c("decrease", "convex"), nknots = 8) 145511.335 172887.618 203786.0940 202997.4795 228688.607 347661.29   100

- Stuart

---

请注意，还有conSpline软件包可以拟合凹或凸样条，并且这比上述选项更快（优化是在Rcpp中完成的）。还有cgam软件包，但速度也相当慢... - Tom Wenseleers

请注意，如果您使用spCobs <- cobs(x, y, constraint = c("decrease", "convex"), knots=seq(min(x),max(x),length.out=6), nknots = 6)的方式指定所需的节点，cobs的速度会快大约10倍，因为这将禁用耗时的自适应节点选择... - Tom Wenseleers

Answer 2

我的另一个回答展示了单调样条和cobs以及scam形状约束样条。这些形状约束样条的问题在于它们相当慢，并且不一定插值所有数据点。

Schumaker样条

我发布了一个实现Schumaker样条的软件包，如果数据是单调和凸/凹的，则该样条是单调和凸/凹的。它速度快并且插值所有数据点。

例如：

#install.packages("schumaker")
library(schumaker)

x = seq(1,10)
y = -log(x)
xarray = seq(1,10,0.01)

SchumSpline = schumaker::Schumaker(x,y)
Schum0 = SchumSpline$Spline(xarray)
Schum1 = SchumSpline$DerivativeSpline(xarray)
Schum2 = SchumSpline$SecondDerivativeSpline(xarray)

plot(xarray, Schum0, type = "l", col = 4, ylim = c(-3,1), main = "Schumaker Spline and first two derivatives",
     ylab = "Spline and derivatives", xlab = "x")
points(x,y)
lines(xarray, Schum1, col = 2)
lines(xarray, Schum2, col = 3)
abline(h = 0, col = 1)
text(x=rep(8,8,8), y=c(-2, -0.5,+0.2), pos=4, labels=c('Spline', 'First Derivative', 'Second Derivative'))

你可以看到二阶导数始终为正（普通单调样条不是这样。请参阅软件包的vignette）。

请注意，只有当数据全局凸/凹时，样条才会全局凸/凹。这是无法避免的，因为这是一个插值样条。

该样条比cobs和scam样条更快。它创建比单调样条慢，但评估速度更快。完整的速度测试可在vignette中找到。

- Stuart