redfit {dplR}  R Documentation 
RedNoise Spectra of TimeSeries
Description
Estimate rednoise spectra from a possibly unevenly spaced timeseries.
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 LombScargle 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 timeseries by using the LombScargle Fourier transform. The
spectrum is biascorrected 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 discretetime 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 ondemand. 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) DistributionFree Statistical Tests. PrenticeHall.
Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating rednoise spectra directly from unevenly spaced paleoclimatic time series. Computers & Geosciences, 28(3), 421–426.
See Also
Examples
# Create a simulated treering width series that has a rednoise
# 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 treering data
# Note the longterm lowfreq 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