Hi, Here is one way you can locate the peaks and troughs of a smoothed function estimate (using the example data from smooth.spline() demo):
##-- example from smooth.spline() y18 <- c(1:3,5,4,7:3,2*(2:5),rep(10,4)) xx <- seq(1,length(y18), len=201) s2 <- smooth.spline(y18) # GCV d1 <- predict(s2, xx, der=1) # We plot the smooth and its first derivative par(mfrow=c(2,1)) plot(y18, main=deparse(s2$call), col.main=2) lines(predict(s2, xx), col = 2) plot(d1$x, d1$y, type="l") abline(h=0) # We can visually pick intervals where first derivative is zero # Using uniroot() to locate the zeros of derivative deriv.x <- function(x, sobj, deriv) predict(sobj, x=x, deriv=deriv)$y uniroot(deriv.x, interval=c(5,8), sobj=s2, deriv=1)$root uniroot(deriv.x, interval=c(8,11), sobj=s2, deriv=1)$root uniroot(deriv.x, interval=c(15,18), sobj=s2, deriv=1)$root Hope this helps, Ravi. -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of baptiste auguie Sent: Friday, August 08, 2008 1:25 PM To: Hans W. Borchers Cc: r-help@r-project.org Subject: Re: [R] re cursive root finding On 8 Aug 2008, at 16:44, Hans W. Borchers wrote: > > As your curve is defined by its points, I don't see any reason to > artificially apply functions such as 'uniroot' or 'optim' (being a > real overkill in this situation). > I probably agree with this, although the process of using a spline leaves me with a "real" function. > First smooth the curve with splines, Savitsky-Golay, or Whittacker > smoothing, etc., then loop through the sequence of points and compute > the gradient by hand, single-sided, two-sided, or both. > > At the same time, mark those indices where the gradient is zero or > changes its sign; these will be the solutions you looked for. > I guess my question is really how to find those indices. predict.smooth.spline can give me the derivative for any x values I want, so I don't need to compute the gradient numerically. However, I still don't know how to locate the zeros automatically. How did you find these values? smooth <- smooth.spline(values$x, values$y) predict(smooth, der=1) -> test which(test$y == 0 ) #gives (obviously) no answer, while which(test$y < 1e-5) # gives many ... > With your example, I immediate got as maxima or minima: > > x1 = 1.626126 > x2 = 4.743243 > x3 = 7.851351 > > // Hans Werner Borchers Thanks, baptiste > > > Any comments? Maybe the problem was not clear or looked too specific. > I'll add a more "bare bones" example, if only to simulate discussion: > >> x <- seq(1,10,length=1000) >> set.seed(11) # so that you get the same randomness y <- >> jitter(sin(x),a = 0.2) values <- data.frame(x= x, y = y) >> >> findZero <- function (guess, neighbors, values) { >> >> smooth <- smooth.spline(values$x, values$y) >> >> obj <- function(x) { >> (predict(smooth, x)$y) ^2 >> } >> minimum <- which.min(abs(values$x - guess)) >> optimize(obj, interval = c(values$x[minimum - neighbors], >> values$x[minimum + neighbors])) # uniroot could be used >> instead i suppose >> >> } >> >> test <- findZero(guess = 6 , neigh = 50, values = values) >> plot(x,y) >> abline(h=0) >> abline(v=test$minimum, col="red") > > Now, I would like to find all (N=)3 roots, without a priori knowledge > of their location in the interval. I considered several approaches: > > 1) find all the numerical values of the initial data that are close to > zero with a given threshold. Group these values in N sets using cut() > and hist() maybe? I could never get this to work, the factors given by > cut confuse me (i couldn't extract their value). Then, apply the > function given above with the guess given by the center of mass of the > different groups of zeros. > > 2) apply findZero once, store the result, then add something big > (1e10) to the neighboring points and look for a zero again and repeat > the procedure until N are found. This did not work, I assume because > it does not perturb the minimization problem in the way I want. > > 3) once a zero is found, look for zeros on both sides, etc... this > quickly makes a complicated decision tree when the number of zeros > grows and I could not find a clean way to implement it. > > Any thoughts welcome! I feel like I've overlooked an obvious trick. > > Many thanks, > > baptiste > > > On 7 Aug 2008, at 11:49, baptiste auguie wrote: > >> Dear list, >> >> >> I've had this problem for a while and I'm looking for a more general >> and robust technique than I've been able to imagine myself. I need to >> find N (typically N= 3 to 5) zeros in a function that is not a >> polynomial in a specified interval. >> >> The code below illustrates this, by creating a noisy curve with three >> peaks of different position, magnitude, width and asymmetry: >> >>> x <- seq(1, 10, length=500) >>> exampleFunction <- function(x){ # create some data with peaks of >>> different scaling and widths + noise >>> fano <- function (p, x) >>> { >>> y0 <- p[1] >>> A1 <- abs(p[2]) >>> w1 <- abs(p[3]) >>> x01 <- p[4] >>> q <- p[5] >>> B1 <- (2 * A1/pi) * ((q * w1 + x - x01)^2/(4 * (x - x01)^2 + >>> w1^2)) >>> y0 + B1 >>> } >>> p1 <- c(0.1, 1, 1, 5, 1) >>> p2 <- c(0.5, 0.7, 0.2, 4, 1) >>> p3 <- c(0, 0.5, 3, 1.2, 1) >>> y <- fano(p1, x) + fano(p2, x) + fano(p3, x) >>> jitter(y, a=0.005*max(y)) >>> } >>> >>> y <- exampleFunction(x) >>> >>> sample.df <- data.frame(x = x, y = y) >>> >>> with(sample.df, plot(x, y, t="l")) # there are three peaks, located >>> around x=2, 4 ,5 y.spl <- smooth.spline(x, y) # smooth the noise >>> lines(y.spl, col="red") >>> >> >> I wish to obtain the zeros of the first and second derivatives of the >> smoothed spline y.spl. I can use uniroot() or optim() to find one >> root, but I cannot find a reliable way to iterate and find the >> desired number of solutions (3 peaks and 2 inflexion points on each >> side of them). I've used locator() or a guesstimate of the disjoints >> intervals to look for solutions, but I would like to get rid off this >> unnecessary user input and have a robust way of finding a predefined >> number of solutions in the total interval. Something along the lines >> of: >> >> findZeros <- function( f , numberOfZeros = 3, exclusionInterval = >> c(0.1 , 0.2, 0.1) >> { >> # >> # while (number of solutions found is less than numberOfZeros) # >> search for a root of f in the union of intervals excluding a >> neighborhood of the solutions already found (exclusionInterval) # } >> >> I could then apply this procedure to the first and second derivatives >> of y.spl, but it could also be useful for any other function. >> >> Many thanks for any remark of suggestion! >> >> baptiste >> > > -- > View this message in context: > http://www.nabble.com/recursive-root-finding-tp18868013p18894331.html > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > 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. _____________________________ Baptiste AuguiƩ School of Physics University of Exeter Stocker Road, Exeter, Devon, EX4 4QL, UK Phone: +44 1392 264187 http://newton.ex.ac.uk/research/emag ______________________________________________ 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. ______________________________________________ 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.