[R] which.minimums not which.min
Philippe Grosjean
phgrosjean at sciviews.org
Mon Mar 20 09:18:56 CET 2006
Fred J. wrote:
>
> Philippe Grosjean <phgrosjean at sciviews.org> wrote: What Fred is looking for is local minima/maxima, also known as turning
> points, or pits/peaks in a series. You can look at ?turnpoints in
> pastecs library.
>
> > x <- c(1:4,0:5, 4, 11)
> > x
> [1] 1 2 3 4 0 1 2 3 4 5 4 11
> > tp <- turnpoints(x)
> > summary(tp)
> Turning points for: x
>
> nbr observations : 12
> nbr ex-aequos : 0
> nbr turning points: 4 (first point is a peak)
> E(p) = 6.666667 Var(p) = 1.811111 (theoretical)
>
> point type proba info
> 1 4 peak 0.100000000 3.3219281
> 2 5 pit 0.002380952 8.7142455
> 3 10 peak 0.005952381 7.3923174
> 4 11 pit 0.666666667 0.5849625
> > plot(tp) # Only useful for a longer and more complex series!
> > # Get the position of peaks
> > (1:length(x))[extract(tp, no.tp = FALSE, peak = TRUE, pit = FALSE)]
> [1] 4 10
> Warning message:
> arguments after the first two are ignored in: UseMethod("extract", e, n,
> ...)
> > (1:length(x))[extract(tp, no.tp = FALSE, peak = FALSE, pit = TRUE)]
> [1] 5 11
> Warning message:
> arguments after the first two are ignored in: UseMethod("extract", e, n,
> ...)
> > # By the way, there are warnings although it works well (I ask on R-Help)
>
> Now, you can easily code your which.minima() function using turnpoints:
>
> x <- c(1:4,0:5, 4, 11)
> x
> tp <- turnpoints(x)
> summary(tp)
> plot(tp) # Only useful for a longer and more complex series!
> # Get the position of peaks
> (1:length(x))[extract(tp, no.tp = FALSE, peak = TRUE, pit = FALSE)]
> (1:length(x))[extract(tp, no.tp = FALSE, peak = FALSE, pit = TRUE)]
> # By the way, there are warnings although it works well (I ask on R-Help)
>
> which.minima <- function(x) {
> if (!require(pastecs)) stop("pastecs library is required!")
> x <- as.vector(x)
> (1:length(x))[extract(turnpoints(x), no.tp = FALSE, peak = FALSE, pit =
> TRUE)]
> }
>
> which.minima(x)
>
> Of course, you could optimize this code. This is just a rough solution!
> Best,
>
> Philippe Grosjean
>
> ..............................................<ï¿½}))><........
> ) ) ) ) )
> ( ( ( ( ( Prof. Philippe Grosjean
> ) ) ) ) )
> ( ( ( ( ( Numerical Ecology of Aquatic Systems
> ) ) ) ) ) Mons-Hainaut University, Pentagone (3D08)
> ( ( ( ( ( Academie Universitaire Wallonie-Bruxelles
> ) ) ) ) ) 8, av du Champ de Mars, 7000 Mons, Belgium
> ( ( ( ( (
> ) ) ) ) ) phone: + 32.65.37.34.97, fax: + 32.65.37.30.54
> ( ( ( ( ( email: Philippe.Grosjean at umh.ac.be
> ) ) ) ) )
> ( ( ( ( ( web: http://www.umh.ac.be/~econum
> ) ) ) ) ) http://www.sciviews.org
> ( ( ( ( (
> ..............................................................
>
> Is there a limit as far as the data input to the function turnpoints? It seams to be unable to process a vector with 15000 records with number in the form of 734.2983, the data is known to have picks and pits
>
> The steps:
> Where pp is the vector with data of the above-mentioned style.
>
> >tp <- turnpoints(rev(pp))
> > mx_id_rev <- (1:length(pp))[extract(tp, no.tp = FALSE, peak = TRUE, pit = FALSE)]
> > tp
> Turning points for: rev(pp)
>
> nbr observations : 11745
> nbr ex-aequos : 10273
> nbr turning points: 1231 (first point is a pit)
> E(p) = 7828.667 Var(p) = 2087.678 (theoretical)
> > mx_id_rev
> numeric(0)
>
> Regards,
There is certainly a limit, but not that low. My test gives (R under Win
XP with 1Gb RAM):
> x <- rnorm(15000)
> library(pastecs)
Loading required package: boot
> tp <-turnpoints(rev(x))
> tp
Turning points for: rev(x)
nbr observations : 15000
nbr ex-aequos : 0
nbr turning points: 10079 (first point is a peak)
E(p) = 9998.667 Var(p) = 2666.344 (theoretical)
> res <- (length(x):1)[extract(tp, no.tp = FALSE, peak = TRUE, pit = FALSE)]
> length(res)
[1] 6250
> res <- (length(x):1)[extract(tp, n= 15000, no.tp = FALSE, peak = TRUE, pit = FALSE)]
> length(res)
[1] 6250
[1] 5040
> plot(x[1:500], type = "l")
> points(res, x[res], col = 2)
That seems to work. So, there is a bug when you don't specify 'n'. I
work on a patch for that and upload the new version of pastecs package
soon. Thanks for pointing me this bug.
Best,
Philippe Grosjean
More information about the R-help
mailing list