Gabor Grothendieck <ggrothendieck <at> myway.com> writes: : : Dr Carbon <drcarbon <at> gmail.com> writes: : : : : : I have a seasonal time series. I want to calculate the annual mean : : value of the time series at its peak : : : : (say the mean of the three values before the peak, the peak, and the : : three values after the peak). : : : : The peak of the time series might change cycle slightly from year to year. : : : : # E.g., : : nPts <- 254 : : foo <- sin((2 * pi * 1/24) * 1:nPts) : : foo <- foo + rnorm(nPts, 0, 0.05) : : bar <- ts(foo, start = c(1980,3), frequency = 24) : : plot(bar) : : start(bar) : : end(bar) : : : : # I want to find the peak value from each year, and then get the mean : : of the values on either side. : : # So, if the peak value in the year 1981 is : : max.in.1981 <- max(window(bar, start = c(1981,1), end = c(1981,24))) : : # e.g, cycle 7 or 8 : : window(bar, start = c(1981,1), end = c(1981,24)) == max.in.1981 : : # E.g. if the highest value in 1981 is in cycle 8 I want : : mean.in.1981 <- mean(window(bar, start = c(1981,5), end = c(1981,11))) : : plot(bar) : : points(ts(mean.in.1981, start = c(1981,8), frequency = 24), col = : : "red", pch = "+") : : : : Is there a way to "automate" this for each year. : : Calculate the moving average of bar, which we call barma, and define a : function f which takes a two column structure, locates : the largest entry in column 1 and returns the corresponding : entry in column 2. Use 'by' to apply f to the two columns, bar and barma : for each year. Finally convert result back to a ts object. : : barma <- filter(bar, rep(1,7)/7) : f <- function(bar) bar[which.max(bar[,1]),2] : barpeakavg <- by(cbind(bar, barma), floor(time(bar)+.0001), f) : barpeakavg.ts <- ts(barpeakavg, start = start(time(aggregate(bar)))) : : : : : How can I return the cycle of the max value by year? : : : : aggregate(bar, FUN = which.max)
Note that this will be off for the first year if that year does not begin at the beginning of the cycle. If thats a problem then use the previous solution but replace barma with barcycle where barcycle <- cycle(bar) ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
