[R] Missing "spline_coef" DLL and Rob Hyndmans monotonic interpolator
Paul.Rustomji at csiro.au
Paul.Rustomji at csiro.au
Mon Jun 2 06:04:46 CEST 2008
Hello R help
I have been trying to use Rob Hyndman's monotonically increasing spline
function. But like another user or two seem have a problem with a
missing DLL (namely "spline_coef"). None of the previous help postings
seemed to have any solutions to this problem. As per a Ripley
suggestion I have deleted all previous versions of R and reinstalled R
2.7.0 and the problem persists.
Thanks
Paul.
x <- seq(0,4,l=20)
y <- sort(rnorm(20))
plot(x,y)
lines(spline(x, y, n = 201), col = 2) # Not necessarily monotonic
lines(cm.spline(x, y, n = 201), col = 3)
>Error in .C("spline_coef", method = as.integer(method), n = nx, x = x,
:
C symbol name "spline_coef" not in DLL for package "stats"
Cm.spline code from
http://www-personal.buseco.monash.edu.au/~hyndman/Rlibrary/interpcode.R
#######################################
cm.spline <- function (x, y = NULL, n = 3 * length(x), method = "fmm",
xmin = min(x), xmax = max(x), gulim=0)
# modification of stats package spline to produce co-monotonic
# interpolant by Hyman Filtering. if gulim!=0 then it is taken as the
upper
# limit on the gradient, and interpolant is gradient limited rather than
# monotonic. Modifications from R stats package spline()
# (c) Simon N. Wood & Rob J. Hyndman. 2002.
{
x <- xy.coords(x, y)
y <- x$y
x <- x$x
nx <- length(x)
method <- match(method, c("periodic", "natural", "fmm"))
if (is.na(method))
stop("spline: invalid interpolation method")
dx <- diff(x)
if (any(dx < 0))
{
o <- order(x)
x <- x[o]
y <- y[o]
}
if(any(diff(y)<0))
stop("Data are not monotonic")
if (method == 1 && y[1] != y[nx])
{
warning("spline: first and last y values differ - using y[1] for
both")
y[nx] <- y[1]
}
z <- .C("spline_coef", method = as.integer(method), n = nx,
x = x, y = y, b = double(nx), c = double(nx), d = double(nx),
e = double(if (method == 1) nx else 0), PACKAGE = "stats")
z$y <- z$y-z$x*gulim # trick to impose upper
z$b <- z$b-gulim # limit on interpolator gradient
z <- hyman.filter(z) # filter gradients for co-monotonicity
z$y <- z$y+z$x*gulim # undo trick
z$b <- z$b+gulim # transformation
z <- spl.coef.conv(z) # force other coefficients to consistency
u <- seq(xmin, xmax, length.out = n)
.C("spline_eval", z$method, nu = length(u), x = u, y = double(n),
z$n, z$x, z$y, z$b, z$c, z$d, PACKAGE = "stats")[c("x","y")]
}
cm.splinefun <- function(x, y = NULL, method = "fmm",gulim=0)
# modification of stats package splinefun to produce co-monotonic
#interpolant by Hyman Filtering. if gulim!=0 then it is taken as the
upper
#limit on the gradient, and intrpolant is gradient limited rather than
# monotonic. Modifications from R stats package splinefun()
# (c) Simon N. Wood 2002
{
x <- xy.coords(x, y)
y <- x$y
x <- x$x
n <- length(x)
method <- match(method, c("periodic", "natural", "fmm"))
if (is.na(method))
stop("splinefun: invalid interpolation method")
if (any(diff(x) < 0))
{
z <- order(x)
x <- x[z]
y <- y[z]
}
if(any(diff(y)<0))
stop("Data are not monotonic")
if (method == 1 && y[1] != y[n])
{
warning("first and last y values differ in spline - using y[1]
for both")
y[n] <- y[1]
}
z <- .C("spline_coef", method = as.integer(method), n = n,
x = as.double(x), y = as.double(y), b = double(n), c =
double(n),
d = double(n), e = double(if (method == 1) n else 0),
PACKAGE = "stats")
z$y <- z$y-z$x*gulim # trick to impose upper
z$b <- z$b-gulim # limit on interpolator gradient
z <- hyman.filter(z) # filter gradients for co-monotonicity
z$y <- z$y+z$x*gulim # undo trick
z$b <- z$b+gulim # transformation
z <- spl.coef.conv(z) # force other coefficients to consistency
rm(x, y, n, method)
function(x)
{
.C("spline_eval", z$method, length(x), x = as.double(x), y =
double(length(x)),
z$n, z$x, z$y, z$b, z$c, z$d, PACKAGE = "stats")$y
}
}
spl.coef.conv <- function(z)
# takes an object z containing equally lengthed arrays
# z$x, z$y, z$b, z$c, z$d defining a cubic spline interpolating
# z$x, z$y and forces z$c and z$d to be consistent with z$y and
# z$b (gradient of spline). This is intended for use in conjunction
# with Hyman's monotonicity filter.
# Note that R's spline routine has s''(x)/2 as c and s'''(x)/6 as d.
# (c) Simon N. Wood
{
n <- length(z$x)
h <- z$x[2:n]-z$x[1:(n-1)]
y0 <- z$y[1:(n-1)];y1 <- z$y[2:n]
b0 <- z$b[1:(n-1)];b1 <- z$b[2:n]
cc <- -(3*(y0-y1)+(2*b0+b1)*h)/h^2
c1 <- (3*(y0[n-1]-y1[n-1])+(b0[n-1]+2*b1[n-1])*h[n-1])/h[n-1]^2
dd <- (2*(y0-y1)/h+b0+b1)/h^2
d1 <- dd[n-1]
z$c <- c(cc,c1);z$d <- c(dd,d1)
z
}
hyman.filter <- function(z)
# Filters cubic spline function to yield co-monotonicity in accordance
# with Hyman (1983) SIAM J. Sci. Stat. Comput. 4(4):645-654, z$x is knot
# position z$y is value at knot z$b is gradient at knot. See also
# Dougherty, Edelman and Hyman 1989 Mathematics of Computation
# 52: 471-494. (c) Simon N. Wood
{
n <- length(z$x)
ss <- (z$y[2:n]-z$y[1:(n-1)])/(z$x[2:n]-z$x[1:(n-1)])
S0 <- c(ss[1],ss)
S1 <- c(ss,ss[n-1])
sig <- z$b
ind <- (S0*S1>0)
sig[ind] <- S1[ind]
ind <- (sig>=0)
if(sum(ind))
z$b[ind] <-
pmin(pmax(0,z$b[ind]),3*pmin(abs(S0[ind]),abs(S1[ind])))
ind <- !ind
if(sum(ind))
z$b[ind] <-
pmax(pmin(0,z$b[ind]),-3*pmin(abs(S0[ind]),abs(S1[ind])))
z
}
#######################################
Paul Rustomji
Rivers and Estuaries
CSIRO Land and Water
GPO Box 1666
Canberra ACT 2601
ph +61 2 6246 5810
mobile 0406 375 739
More information about the R-help
mailing list