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

[David Winsemius]

On Feb 15, 2010, at 1:06 PM, Tim Heaton wrote:

> 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 see your point.

FailMonSpline(seq(1, 8, by=0.1))-FailMonSpline(0.1+seq(1, 8, by=0.1))
  produces 3 positive numbers at indices 35-37.

-- 
David

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

David Winsemius, MD
Heritage Laboratories
West Hartford, CT