[R] re cursive root finding

Ravi Varadhan RVaradhan at jhmi.edu
Fri Aug 8 21:22:52 CEST 2008 

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.

More information about the R-help mailing list