[R] recursive root finding

[baptiste auguie]

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



_____________________________

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