[R] Monotonic interpolation

[Simon Wood]

> Has anyone got a function for smooth monotonic interpolation of a
> univariate function?  I'm after something like the NAG function PCHIM
> which does monotonic Hermite interpolation. Alternatively, montononic
> cubic spline interpolation.

cm.splinefun() pasted below fits an interpolating spline and filters the
coefficients for (co)monotonicity following Hyman's 1983 algorithm. It
works pretty much as base function cm.splinefun, from which it is
modified.

Simon

  ______________________________________________________________________
> Simon Wood  snw at st-and.ac.uk  http://www.ruwpa.st-and.ac.uk/simon.html
> The Mathematical Institute, North Haugh, St. Andrews, Fife KY16 9SS UK
> Direct telephone: (0)1334 463799          Indirect fax: (0)1334 463748 



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
}

cm.splinefun<-function(x, y = NULL, method = "fmm",gulim=0) 
# modification of base 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 base 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 (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 = "base")
    
    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 = "base")$y
    }
}


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._