[R] re cursive root finding

baptiste auguie ba208 at exeter.ac.uk
Sat Aug 9 01:46:26 CEST 2008 

Hi,

my problem is how to get rid off the "visual decision" step and have  
the procedure automatically search the interval in some way until a  
desired number of roots is found.

Thanks,

baptiste

On 8 Aug 2008, at 20:22, Ravi Varadhan wrote:

> 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: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org 
> ] On
> Behalf Of baptiste auguie
> Sent: Friday, August 08, 2008 1:25 PM
> To: Hans W. Borchers
> Cc: r-help at 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 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
>
> ______________________________________________
> 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

More information about the R-help mailing list