要求花键凸出

Question

我需要将样条曲线拟合到一组数据，并且需要使结果函数单调递减和凸出。我传递给splinefun的数据可以保证具有这些属性，但这不能保证结果函数是凸的。有什么方法可以将样条曲线拟合到一组数据，并要求结果函数为凸函数？

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包。然后，我们可以:

# 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)

请注意以下几点:

选项bs="mdcx"表示我们希望减少凸的样条曲线。如果要增加凸面，减少凹面等，请查找适当的bs值here。

如果将非Montonic数据放入splinefun(x,y,"monoH.FC")函数中，则会出现错误。

如果将非凸数据放入scam函数中，则仍然会得到样条曲线。数据被改变，使得小的凸面被改变为凸面。对此没有警告，但是请小心，因为您的数据看起来可能完全不同。例如，下面的绘图是使用与上面相同的代码制作的，不同的是我们使用bs="mdcv"来减少凹面函数:

单调凸花键(带穗子包装)

# 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样条和骗局样条的初始计算时间要长得多。

。

# 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