redfit {dplR} | R Documentation |
Red-Noise Spectra of Time-Series
Description
Estimate red-noise spectra from a possibly unevenly spaced time-series.
Usage
redfit(x, t, tType = c("time", "age"), nsim = 1000, mctest = TRUE,
ofac = 4, hifac = 1, n50 = 3, rhopre = NULL,
p = c(0.10, 0.05, 0.02), iwin = 2,
txOrdered = FALSE, verbose = FALSE, seed = NULL,
maxTime = 10, nLimit = 10000)
runcrit(n, p = c(0.10, 0.05, 0.02), maxTime = 10, nLimit = 10000)
Arguments
x |
a |
t |
a |
tType |
a |
nsim |
a |
mctest |
a |
ofac |
oversampling factor for Lomb-Scargle Fourier transform.
A |
hifac |
maximum frequency to analyze relative to the Nyquist
frequency. A |
n50 |
number of segments. The segments overlap by about 50 percent. |
rhopre |
a |
p |
a |
iwin |
the type of window used for scaling the values of each
segment of data. A |
txOrdered |
a |
verbose |
a |
seed |
a value to be given to |
maxTime |
a |
nLimit |
a |
n |
an integral value giving the length of the sequence in the number of runs test. |
Details
Function redfit
computes the spectrum of a possibly unevenly
sampled time-series by using the Lomb-Scargle Fourier transform. The
spectrum is bias-corrected using spectra computed from simulated AR1
series and the theoretical AR1 spectrum.
The function duplicates the functionality of program REDFIT by Schulz and Mudelsee. See the manual of that program for more information. The results of this function should be very close to REDFIT. However, some changes have been made:
More precision is used in some constants and computations.
All the data are used: the last segment always contains the last pair of (t, x). There may be small differences between
redfit
and REDFIT with respect to the number of points per segment and the overlap of consecutive segments.The critical values of the runs test (see the description of
runcrit
below) differ betweenredfit
and REDFIT. The approximate equations in REDFIT produce values quite far off from the exact values when the number of frequencies is large.The user can select the significance levels of the runs test.
Most of the window functions have been adjusted.
6 dB bandwidths have been computed for discrete-time windows.
Function runcrit
computes the limits of the acceptance region
of a number of runs test: assuming a sequence of n
i.i.d.
discrete random variables with two possible values a and b
of equal probability (0.5), we are examining the distribution of the
number of runs. A run is an uninterrupted sequence of only a or
only b. The minimum number of runs is 1 (a sequence with only
a or only b) while the maximum number is n
(alternating a and b). See Bradley, p. 253–254,
259–263. The function is also called from redfit
;
see rcnt
in ‘Value’ for the interpretation. In
this case the arguments p
, maxTime
and
nLimit
are passed from redfit
to runcrit
,
and n
is the number of output frequencies.
The results of runcrit
have been essentially precomputed for
some values of p
and n
. If a precomputed
result is not found and n
is not too large
(nLimit
, maxTime
), the exact results are
computed on-demand. Otherwise, or if package "gmp"
is not
installed, the normal distribution is used for approximation.
Value
Function runcrit
returns a list
containing
rcritlo
, rcrithi
and rcritexact
(see below). Function redfit
returns a list
with the
following elements:
varx |
variance of |
rho |
average autocorrelation coefficient, either estimated
from the data or prescribed ( |
tau |
the time scale of an AR1 process corresponding to
|
rcnt |
a |
rcritlo |
a |
rcrithi |
a |
rcritexact |
a |
freq |
the frequencies used. A |
gxx |
estimated spectrum of the data (t, x). A
|
gxxc |
red noise corrected spectrum of the data. A
|
grravg |
average AR1 spectrum over |
gredth |
theoretical AR1 spectrum. A |
corr |
a |
ci80 |
a |
ci90 |
a |
ci95 |
95th percentile red noise spectrum. |
ci99 |
99th percentile red noise spectrum. |
call |
the call to the function. See |
params |
A
|
vers |
version of dplR. See |
seed |
value of the |
t |
if duplicated values of |
x |
if duplicated values of |
Author(s)
Mikko Korpela. Examples by Andy Bunn.
References
Function redfit
is based on the Fortran program
REDFIT
(version 3.8e), which is in the public domain.
Bradley, J. V. (1968) Distribution-Free Statistical Tests. Prentice-Hall.
Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. Computers & Geosciences, 28(3), 421–426.
See Also
Examples
# Create a simulated tree-ring width series that has a red-noise
# background ar1=phi and sd=sigma and an embedded signal with
# a period of 10 and an amplitude of have the rednoise sd.
library(graphics)
library(stats)
runif(1)
rs <- .Random.seed
set.seed(123)
nyrs <- 500
yrs <- 1:nyrs
# Here is an ar1 time series with a mean of 2mm,
# an ar1 of phi, and sd of sigma
phi <- 0.7
sigma <- 0.3
sigma0 <- sqrt((1 - phi^2) * sigma^2)
x <- arima.sim(list(ar = phi), n = nyrs, sd = sigma0) + 2
# Here is a sine wave at f=0.1 to add in with an amplitude
# equal to half the sd of the red noise background
per <- 10
amp <- sigma0 / 2
wav <- amp * sin(2 * pi / per * yrs)
# Add them together so we have signal and noise
x <- x + wav
# Here is the redfit spec
redf.x <- redfit(x, nsim = 500)
# Acceptance region of number of runs test
# (not useful with default arguments of redfit())
runcrit(length(redf.x[["freq"]]))
op <- par(no.readonly = TRUE) # Save to reset on exit
par(tcl = 0.5, mar = rep(2.2, 4), mgp = c(1.1, 0.1, 0))
plot(redf.x[["freq"]], redf.x[["gxxc"]],
ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),
type = "n", ylab = "Spectrum", xlab = "Frequency (1/yr)",
axes = FALSE)
grid()
lines(redf.x[["freq"]], redf.x[["gxxc"]], col = "#1B9E77")
lines(redf.x[["freq"]], redf.x[["ci99"]], col = "#D95F02")
lines(redf.x[["freq"]], redf.x[["ci95"]], col = "#7570B3")
lines(redf.x[["freq"]], redf.x[["ci90"]], col = "#E7298A")
freqs <- pretty(redf.x[["freq"]])
pers <- round(1 / freqs, 2)
axis(1, at = freqs, labels = TRUE)
axis(3, at = freqs, labels = pers)
mtext(text = "Period (yr)", side = 3, line = 1.1)
axis(2); axis(4)
legend("topright", c("x", "CI99", "CI95", "CI90"), lwd = 2,
col = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"),
bg = "white")
box()
## Not run:
# Second example with tree-ring data
# Note the long-term low-freq signal in the data. E.g.,
# crn.plot(cana157)
library(utils)
data(cana157)
yrs <- time(cana157)
x <- cana157[, 1]
redf.x <- redfit(x, nsim = 1000)
plot(redf.x[["freq"]], redf.x[["gxxc"]],
ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),
type = "n", ylab = "Spectrum", xlab = "Frequency (1/yr)",
axes = FALSE)
grid()
lines(redf.x[["freq"]], redf.x[["gxxc"]], col = "#1B9E77")
lines(redf.x[["freq"]], redf.x[["ci99"]], col = "#D95F02")
lines(redf.x[["freq"]], redf.x[["ci95"]], col = "#7570B3")
lines(redf.x[["freq"]], redf.x[["ci90"]], col = "#E7298A")
freqs <- pretty(redf.x[["freq"]])
pers <- round(1 / freqs, 2)
axis(1, at = freqs, labels = TRUE)
axis(3, at = freqs, labels = pers)
mtext(text = "Period (yr)", side = 3, line = 1.1)
axis(2); axis(4)
legend("topright", c("x", "CI99", "CI95", "CI90"), lwd = 2,
col = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"),
bg = "white")
box()
par(op)
## End(Not run)
.Random.seed <- rs