要求样条曲线凸出

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)

请注意以下事项:

1. 该选项bs="mdcx"意味着我们需要减少凸花键.如果你想增加凸面,减少凹面等,那么在这里查找适当的bs值.
2. 如果将非单调数据放入splinefun(x,y,"monoH.FC")函数中,则会出现错误.
3. 如果您将非凸数据放入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样条曲线和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