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