re cursive root finding

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

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.

With your example, I immediate got as maxima or minima:

   x1 = 1.626126
   x2 = 4.743243
   x3 = 7.851351

//  Hans Werner Borchers

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

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.