Bug 14215

Summary: Error in coding for splinefun(method = "monoH.FC")
Product: R Reporter: Jitterbug compatibility account <jitterbug-import>
Component: AnalysesAssignee: Martin Maechler <maechler>
Status: RESOLVED FIXED    
Severity: normal CC: jackie.rosen, maechler
Priority: P5    
Version: old   
Hardware: All   
OS: Linux   

Description Jitterbug compatibility account 2010-02-18 21:08:11 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)

Comment 1 Jitterbug compatibility account 2010-02-25 17:37:46 UTC
Audit (from Jitterbug):
Thu Feb 25 11:37:46 2010	ripley	moved from incoming to Analyses
Comment 2 Brian Ripley 2010-04-11 17:34:06 UTC
I think this was known, but the fix is not (and it is not obvious
that it is worth the work it might entail).
Comment 3 Martin Maechler 2010-04-19 10:58:23 UTC
Well, the problem was known to Fritsch & Carlson, and I think our R code must be fixed here.
Comment 4 Martin Maechler 2010-04-23 07:33:53 UTC
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)
Comment 5 Jackie Rosen 2014-02-16 11:42:57 UTC
(spam comment removed)