[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