[R] Non-monotonic spline using splinefun(method = "monoH.FC")

[David Winsemius]

On Feb 15, 2010, at 10:59 AM, Tim Heaton wrote:

> Hi,
>
> In my version of R, the stats package splinefun code for fitting a
> Fritsch and Carlson monotonic spline does not appear to guarantee a
> monotonic result. If two adjoining sections both have over/undershoot
> the way the resulting adjustment of alpha and beta is performed can  
> give
> modified values which still do not satisfy the required constraints. I
> do not think this is due to finite precision arithmetic. Is this a  
> known
> bug? Have had a look through the bug database but couldn't find  
> anything.

IThe help page says that the resulting function will be "monotone  
<bold-on> iff <bold-off> the data are."

y[7] < y[8]  # False
FailMonSpline[7] < FailMonSpline[8] # False, ... , as promised.

>
> Below is an example created to demonstrate this,
>
> ###############################################
> # Create the following data
> # This is created so that their are two adjoining sections which  
> have to
> be adjusted
> x <- 1:8
> y <- c(-12, -10, 3.5, 4.45, 4.5, 140, 142, 142)
>
> # Now run the splinefun() function
>
> FailMonSpline <- splinefun(x, y, method = "mono")
>
> # In theory this should be monotonic increasing but the required
> conditions are not satisfied
>
> # Check values of alpha and beta for this curve
> m <- FailMonSpline(x, deriv = 1)
> nx <- length(x)
> n1 <- nx - 1L
> dy <- y[-1] - y[-nx]
> dx <- x[-1] - x[-nx]
> Sx <- dy/dx
>
> alpha <- m[-nx]/Sx
> beta <- m[-1]/Sx
> a2b3 <- 2 * alpha + beta - 3
> ab23 <- alpha + 2 * beta - 3
> ok <- (a2b3 > 0 & ab23 > 0)
> ok <- ok & (alpha * (a2b3 + ab23) < a2b3^2)
> # If the curve is monotonic then all ok should be FALSE however this  
> is
> not the case
> ok
>
>
> # Alternatively can easily seen to be non-monotonic by plotting the
> region between 4 and 5
>
> t <- seq(4,5, length = 200)
> plot(t, FailMonSpline(t), type = "l")
>
> ########################################################
> The version of R I am running is
>
> platform       x86_64-suse-linux-gnu
> arch           x86_64
> os             linux-gnu
> system         x86_64, linux-gnu
> status
> major          2
> minor          8.1
> year           2008
> month          12
> day            22
> svn rev        47281
> language       R
> version.string R version 2.8.1 (2008-12-22)
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
Heritage Laboratories
West Hartford, CT