Dear list,
I've recently become interested in comparing the spectral estimates
using the different methods ("pgram" and "ar") in the spectrum()
function in the stats package.
With many thanks to the authors of these complicated functions, I
would like to point out what looks to me like a bit of an
inconsistency -- but I would not be surprised if there is good
reasoning that justifies it that I am just not seeing right now. If we
use the lh data, the two methods return similar results:
> spectrum(lh, col = "blue")
> spec.ar(lh, add = TRUE)
But using the ldeaths data:
> spectrum(ldeaths, col = "blue")
> spec.ar(ldeaths, add = TRUE)
the resulting plots do not compare over the same frequency range. This
results because spec.ar defines frequency as
> freq <- seq.int(0, 0.5, length.out = n.freq)
whereas spec.pgram uses
> xfreq <- frequency(x)
> N <- nrow(x)
> Nspec <- floor(N/2)
> freq <- seq.int(from = xfreq/N, by = xfreq/N, length.out = Nspec)
And so the reason the spectral estimates of lh are similar is that
frequency(lh) = 1, whereas frequency(ldeaths) = 12.
The documentation seems more extensive for spec.pgram (and the
pertinent section in MASS focuses on spec.pgram), and I realize that
there is a warning in ?spec.ar that AR spectra can be misleading. But
is there a reason that I am not aware of that the frequencies of the
AR spectra are defined in this way? It seems to me that it would be
desirable for frequency to be defined over the same range as in
spec.pgram. All that would need to be added would be a line to scale
the freq vector using the sampling frequency before it is returned.
Eric Thompson
Graduate Student
Dept. of Civil & Env. Eng.
Tufts University
> sessionInfo()
R version 2.6.0 (2007-10-03)
i686-pc-linux-gnu
locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] datasets utils stats graphics grDevices methods base
other attached packages:
[1] MASS_7.2-37
loaded via a namespace (and not attached):
[1] rcompgen_0.1-17
______________________________________________
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.
