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