[R] recursive root finding

[baptiste auguie]

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
>
>
>
> _____________________________
>
> 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 at 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.