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@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
______________________________________________
R-help@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.
