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

[Tim Heaton]

Hi,

 Thanks for the reply but I think you are confusing monotonic and
strictly increasing/decreasing. I also just used the y-value of the last
knot as a simple example as that is not the bit where it goes wrong. It
will still produce a non-monotonic spline if you use for example

x <- 1:8
y <- c(-12, -10, 3.5, 4.45, 4.5, 140, 142, 143)

I am pretty sure that it's a bug in the way that the alpha's and beta's
are modified in the code itself which does not guarantee (if there are
two overlapping sections which need their alphas and betas modifying)
that after modification they satisfy the constraints as explained in the
original Fritsch and Carlson paper. The original paper is quite vague
about how to deal with multiple sections which need modifying --- should
one do it in order (in which case one might get a different result if
one entered the data in the opposite direction and moving one knot would
no longer guarantee that the curve changes in only a finite region) or
possibly shrink the coefficients twice (which would create a flatter
spline than necessary but would give a finite effect and the same curve
if the data were entered in the opposite direction).

Tim


David Winsemius wrote:
> 
> 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
>