Tim Heaton
2010-Feb-15 15:59 UTC
[R] Non-monotonic spline using splinefun(method = "monoH.FC")
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. 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)
David Winsemius
2010-Feb-15 17:44 UTC
[R] Non-monotonic spline using splinefun(method = "monoH.FC")
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
Possibly Parallel Threads
- Error in coding for splinefun(method = "monoH.FC") (PR#14215)
- Anova in 'car': "SSPE apparently deficient rank"
- Type II and III sum of squares (R and SPSS)
- CRAN mirror for R in India: new one at WBUT, how do we get listed in the CRAN website?
- Monotonic interpolation