Bugzilla – Bug 14215

Error in coding for splinefun(method = "monoH.FC")

Last modified: 2014-02-16 11:42:57 UTC

From: t.heaton@shef.ac.uk Full_Name: Tim Heaton Version: 2.8.1 OS: linux-gnu Submission from: (NULL) (143.167.4.162) Hi, In my version of R, the stats package splinefun code for fitting a Fritsch and Carlson monotonic spline does not 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 posed the question as to whether this was a known error on the R help but got no reply, have also 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)

Audit (from Jitterbug): Thu Feb 25 11:37:46 2010 ripley moved from incoming to Analyses

I think this was known, but the fix is not (and it is not obvious that it is worth the work it might entail).

Well, the problem was known to Fritsch & Carlson, and I think our R code must be fixed here.

Fixed, by using *sequential* code --- and hence for speed reasons C --- instead of the "sometimes not quite correct" vectorized R code. For R-devel (2.12.0), to be backported to "R 2.11.0 patched" (-> R 2.11.1)

(spam comment removed)