[R] re cursive root finding

[baptiste auguie]

Hi all,

with a lot more fiddling around, I've come up with this solution:

> x <- seq(1, 15, length=100)
> y <- jitter(cos(x), a=0.2)
>
> plot(x, y)
>
> findZero <- function(x, y, guess = mean(range(x)), exclusion=  
> 5*diff(range(x))/length(x))
> 	{
> 	interval <- c(guess - exclusion, guess+exclusion)
> 	print(interval)
>    spl <- smooth.spline(x, y)
> 	obj <- function(xx) {
> 			z <- predict(spl, xx,  der=0)$y
> 			sum((z)^2)
> 			}
> 		optimize(obj, interval)
> 	}
>
> # findZero(x, y) # test
>
> findZeros <- function(x, y, ...){
> 	spl <- smooth.spline(x, y)
> 	y.spl <- predict(spl, x, der=0)$y
> 	diff(diff(abs(y.spl), 1) > 0, 1)->test
> 	guesses <- x[which(test==1)+1]
> 	sapply(guesses, function(guess) findZero(x, y, guess, ...)$minimum)
> }
>
> zeros <- findZeros(x, y)
> abline(v=zeros, col="red")

albeit not elegant or fast, nor very robust I suspect, it seems to  
work for me.

Best wishes

baptiste

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

_____________________________

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