Hello,
I don't see what's wrong with turnpoints() from pastecs. It is easy to
use, and provides additional information for each turnpoints, i.e.,
probability of occurrence against the null hypothesis that the series is
purely random, and the number of bits of information associated with the
points according to Kendall's information theory. See ?turnpoints.
Rewriting the function is a nice exercise, but a small explanation on
how to use turnpoints() is much easier. So, nobody is able to tell that
simply using:
> library(pastecs)
> turnpoints(dat$count)
does the job? Am I the only one interested by the extra information
provided by turnpoints()?
Here is a more extensive example that shows also how to get date or
counts associated to pits/peaks:
txt <- "
y m d count
93 02 07 3974.6
93 02 08 3976.7
93 02 09 3955.2
93 02 10 3955.0
93 02 11 3971.8
93 02 12 3972.8
93 02 13 3961.0
93 02 14 3972.8
93 02 15 4008.0
93 02 16 4004.2
93 02 17 3981.2
93 02 18 3996.8
93 02 19 4028.2
93 02 20 4029.5
93 02 21 3953.4
93 02 22 3857.3
93 02 23 3848.3
93 02 24 3869.8
93 02 25 3898.1
93 02 26 3920.5
93 02 27 3936.7
93 02 28 3931.9
"
con <- textConnection(txt)
dat <- read.table(con, header = TRUE)
close(con)
dat$date <- as.Date(paste(dat$y, dat$m, dat$d), format = "%y %m %d")
library(pastecs)
tp <- turnpoints(dat$count)
tp
summary(tp)
# Indicate which turnpoints are significant (see ?turnpoints)
plot(tp, level = 0.05)
# Another plot
plot(dat$count, type = "l")
lines(tp)
# Get counts for all turnpoints
allcounts <- <- dat$count[extract(tp, no.tp = FALSE, peak = TRUE, pit =
TRUE)]
# Get dates for all turnpoints
alldates <- dat$date[extract(tp, no.tp = FALSE, peak = TRUE, pit = TRUE)]
alldates
# Get dates for "informative" turnpoints (5%) only (see ?turnpoints)
alldates[tp$proba < 0.05]
# Get dates for peaks only
dat$date[extract(tp, no.tp = FALSE, peak = TRUE, pit = FALSE)]
# Etc...
Best,
Philippe Grosjean
..............................................<°}))><........
) ) ) ) )
( ( ( ( ( Prof. Philippe Grosjean
) ) ) ) )
( ( ( ( ( Numerical Ecology of Aquatic Systems
) ) ) ) ) Mons-Hainaut University, Belgium
( ( ( ( (
..............................................................
(Ted Harding) wrote:
On 17-Sep-09 08:10:47, ogbos okike wrote:
Good morning once more. My problem of yesterday has been addressed.
Having learned a few tricks from that, I wish to ask another question
in connection with that. My data is a cosmic ray data consisting of
dates and counts.
When I plot a graph of counts versus dates, the resultant signal
shows a number of maximum and minimum points. These minimum points
(turning points) are of interest to me. Reading these dates and counts
off from the plot is difficult as I am dealing with a large data.
I have been looking at turnpoints function in pastecs library but have
not been able to figure out the appropriate commands that one can use
to find the minima/maxima (turning points) or pits/peaks in a series.
My data is of the form shown below where y stands for year, m month,
d day and finally count. Is there a way I could find these minima
together with the dates they occurred?
I would be indebted to those of you who will show me the way out of
these problem.
Thank you.
Best regards
Ogbos
y m d count
93 02 07 3974.6
93 02 08 3976.7
93 02 09 3955.2
93 02 10 3955.0
93 02 11 3971.8
93 02 12 3972.8
93 02 13 3961.0
93 02 14 3972.8
93 02 15 4008.0
93 02 16 4004.2
93 02 17 3981.2
93 02 18 3996.8
93 02 19 4028.2
93 02 20 4029.5
93 02 21 3953.4
93 02 22 3857.3
93 02 23 3848.3
93 02 24 3869.8
93 02 25 3898.1
93 02 26 3920.5
93 02 27 3936.7
93 02 28 3931.9
The following simple function TP() (for "Turning Point") locates
the positions i where x[i] is greater than both of its immediate
neighbours (local maximum) or less than both of its neighbours
(local minimum).
TP <- function(x){
L <- length(x)
which( ((x[1:(L-2)]<x[2:(N-1)])&(x[2:(L-1)]>x[3:L]))
|((x[1:(L-2)]>x[2:(N-1)])&(x[2:(L-1)]<x[3:L])) ) + 1
}
Applied to your series "count" above:
TP(count)
# [1] 2 4 6 7 9 11 14 17 21
If you assign these values to an index:
ix <- TP(count)
rbind(d[ix],count[ix])
# [1,] 8.0 10 12.0 13 15 17.0 20.0 23.0 27.0
# [2,] 3976.7 3955 3972.8 3961 4008 3981.2 4029.5 3848.3 3936.7
Of course, this is only a very simplistic view of "turning point",
and will pick out everything which is a local minimum or maximum.
The above function can be extended (in a fairly obvious way) to
identify each position i where x[i] is greater than its neighbours
out to 2 on either side, or less than these neighbours; or more
generally out to k on either side.
A lot depends on how you want to interpret "turning point". With
your "count" series, it might be that you were only interested
in identifying the relatively extreme turning points, such as
i=4 (maybe), i=9 (maybe), i=14, i=17, i=21(maybe).
Hoping this helps,
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.hard...@manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 17-Sep-09 Time: 09:47:53
------------------------------ XFMail ------------------------------
______________________________________________
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.
______________________________________________
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.
