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