**R**'s

**spectrum()**function) I've been trying to visualize in an intuitive, aesthetically appealing way. At first, I just used lattice's

**bwplot**, but the spacing of the X-axis here really matters. The spectra's frequencies aren't regularly-spaced categories, which is the default of

**bwplot**. If all the series are of the same length, then the coefficients of the are all estimated at the same frequencies, and one can use

**plyr**'s split-apply-combine magic to compute the distribution of coefficients at each frequency.

I spent some time trying to faithfully reproduce the plotting layout of the figures produced by, e.g.

**spectrum(rnorm(1e3))**with respect to the axes. This was way more annoying than i expected...

To include a generic example, I started looking around for easily generated, interesting time series. The logistic map is so damned easy to compute and generally pretty interesting, capable of producing a range of frequencies. Still, I was rather surprised by the final output. Playing with these parameters produces a range of serious weirdness -- apparent ghosting, interesting asymmetries, etc..

Note that I've just used the logistic map for an interesting, self-contained example. You can use the exact same code to generate your own

**spectrum()**box-and-whisker plot by substituting your matrix of time series for

**mytimeseries**below. Do note here that series are by row, which is not exactly standard. For series by column, just change the margin from 1 to 2 in the call to

**adply**below.

## box-and-whisker plot of FFT power spectrum by frequency/period ## vary these parameters -- r as per logistic map ## see http://en.wikipedia.org/wiki/Logistic_map for details logmap = function(n, r, x0=0.5) { ret = rep(x0,n) for (ii in 2:n) { ret[ii] = r*ret[ii-1]*(1-ret[ii-1]) } return(ret) } ## modify these for interesting behavior differences rfrom = 3.4 rto = 3.7 rsteps=200 seqlen=1e4 require(plyr) require(ggplot2) mytimeseries = aaply(seq(from=rfrom, to=rto, length.out=rsteps), 1, function(x) { logmap(seqlen, x) }) ## you can plug in any array here for mytimeseries ## each row is a timeseries ## for series by column, change the margin from 1 to 2, below logspec = adply( mytimeseries, 1, function(x) { ## change spec.pgram parameters as per goals ret = spectrum(x, taper=0, demean=T, pad=0, fast=T, plot=F) return( data.frame(freq=ret$freq, spec=ret$spec)) }) ## boxplot stats at each unique frequency logspec.bw = ddply(logspec, 'freq', function(x) { ret = boxplot.stats(log10(x$spec)); ## sometimes $out is empty, bind NA to keep dimensions correct return(data.frame(t(10^ret$stats), out=c(NA,10^ret$out))) }) ## plot -- outliers are dots ## be careful with spacing of hard returns -- ggplot has odd rules ## as close to default "spectrum" plot as possible ggplot(logspec.bw, aes(x=1/(freq))) + geom_ribbon(aes(ymin=X1, ymax=X5), alpha=0.35, fill='green') + geom_ribbon(aes(ymin=X2, ymax=X4), alpha=0.5, fill='blue') + geom_line(aes(y=X3)) + geom_point(aes(y=out), color='darkred') + scale_x_continuous( trans=c('inverse'), name='Period') + scale_y_continuous(name='Power', trans='log10')

Following my nose further on the logistic map example, I also used the

**animate**package to make a movie that walks through the power spectrum of the map for a range of r values. Maybe one day I'll set this to music and post it to Youtube! Though I think some strategic speeding up and slowing down in important parameter regions is warranted for a truly epic tour of the logistic map's power spectrum.

require(animation) saveVideo({ ani.options(interval = 1/5, ani.type='png', ani.width=1600, ani.height=1200) for (ii in c(seq(from=3.2, to=3.5, length.out=200), seq(from=3.5, to=4, length.out=1200))) { spectrum(logmap(1e4, ii), taper=0, demean=T, pad=0, fast=T, sub=round(ii, 3)) } }, video.name = "logmap2.mp4", other.opts = "-b 600k")

## No comments:

## Post a Comment