[Dplr-commits] r679 - in branches/redfit: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 4 10:00:54 CEST 2013
Author: mvkorpel
Date: 2013-09-04 10:00:54 +0200 (Wed, 04 Sep 2013)
New Revision: 679
Modified:
branches/redfit/R/redfit.R
branches/redfit/man/redfit.Rd
Log:
In redfit.R:
* rcritlo and rcrithi are returned by redfit(), not computed in print.redfit().
* rcnt is not computed if the conditions of the statistical test are not met.
* small optimizations in print.redfit()
In redfit.Rd:
* Added reference, Bendat & Piersol: Random Data. Still need to find
the book and check facts.
* Documented rcritlo, rcrithi and rhopre.
Modified: branches/redfit/R/redfit.R
===================================================================
--- branches/redfit/R/redfit.R 2013-09-03 15:13:48 UTC (rev 678)
+++ branches/redfit/R/redfit.R 2013-09-04 08:00:54 UTC (rev 679)
@@ -311,16 +311,34 @@
ci99 <- NULL
}
-
## Test equality of theoretical AR1 and estimated spectrum using a
- ## runs test (Bendat and Piersol, 1986, p. 95).
- rcnt <- 1 + sum(diff(sign(gxxc - gredth)) != 0)
+ ## runs test (Bendat and Piersol, 1986, p. 95). The empirical
+ ## equations for calculating critical values for 5-% significance
+ ## were derived from the tabulated critical values in B&P.
+ if (iwin2 == 0 && ofac == 1 && dn50 == 1) {
+ rcnt <- 1 + sum(diff(sign(gxxc - gredth)) != 0)
+ ## dplR: NOTE! Integer division is used in REDFIT. This should be
+ ## checked (by finding a copy of Bendat and Piersol). For now, we
+ ## can assume that real(nout/2) was supposed to be real(nout)/2.
+ ## sqrtHalfNfreq <- sqrt(nfreq %/% 2)
+ sqrtHalfNfreq <- sqrt(nfreq / 2)
+ ## dplR: NOTE! Is round() the right function to use? Maybe floor()
+ ## for the lower limit and ceiling for the higher limit?
+ rcritlo <- round((-0.79557086 + 1.0088719 * sqrtHalfNfreq)^2)
+ rcrithi <- round(( 0.75751462 + 0.9955133 * sqrtHalfNfreq)^2)
+ } else {
+ rcnt <- NULL
+ rcritlo <- NULL
+ rcrithi <- NULL
+ }
## dplR: Elements of the list returned from this function:
## varx data variance estimated from spectrum
## rho average autocorrelation coefficient (estimated or prescribed)
## tau average tau, tau == -avgdt / log(rho)
## rcnt runs count, test of equality of theoretical and data spectrum
+ ## rcritlo critical low value for rcnt
+ ## rcrithi critical high value for rcnt
## freq frequency vector
## gxx autospectrum of input data
## gxxc corrected autospectrum of input data
@@ -331,10 +349,10 @@
## ci90 90% false-alarm level from MC
## ci95 95% false-alarm level from MC
## ci99 99% false-alarm level from MC
- ## call dplR: how the function was called
- ## params dplR: parameters dependent on the command line arguments
- ## vers dplR: version of dplR containing the function
- ## seed dplR: if not NULL, value used for set.seed(seed)
+ ## call how the function was called
+ ## params parameters dependent on the command line arguments
+ ## vers version of dplR containing the function
+ ## seed if not NULL, value used for set.seed(seed)
dplrNS <- tryCatch(getNamespace("dplR"), error = function(...) NULL)
if (!is.null(dplrNS) && exists("redfit", dplrNS) &&
identical(match.fun(as.list(cl)[[1]]), get("redfit", dplrNS))) {
@@ -343,6 +361,7 @@
vers <- NULL
}
res <- list(varx = varx, rho = rho, tau = tau, rcnt = rcnt,
+ rcritlo = rcritlo, rcrithi = rcrithi,
freq = freq, gxx = gxx, gxxc = gxxc, grravg = grravg,
gredth = gredth, corr = corr,
ci80 = ci80, ci90 = ci90, ci95 = ci95, ci99 = ci99,
@@ -396,11 +415,7 @@
params <- x[["params"]]
iwin <- params[["iwin"]]
n50 <- params[["n50"]]
- nseg <- params[["nseg"]]
- ofac <- params[["ofac"]]
- rhopre <- params[["rhopre"]]
mctest <- params[["mctest"]]
- nfreq <- params[["nfreq"]]
gredth <- x[["gredth"]]
## scaling factors for red noise from chi^2 distribution
@@ -413,27 +428,6 @@
fac95 <- qchisq(0.95, dof) / dof
fac99 <- qchisq(0.99, dof) / dof
- ## critical false alarm level after Thomson (1990)
- ## dplR: modified from original REDFIT code to accommodate for
- ## lower / upper tail difference
- alphacrit <- (nseg - 1) / nseg
- faccrit <- qchisq(alphacrit, dof) / dof
-
- ## Test equality of theoretical AR1 and estimated spectrum using a
- ## runs test (Bendat and Piersol, 1986, p. 95). The empirical
- ## equations for calculating critical values for 5-% significance
- ## were derived from the tabulated critical values in B&P.
- ##
- ## dplR: NOTE! Integer division is used in REDFIT. This should be
- ## checked (by finding a copy of Bendat and Piersol). For now, we
- ## can assume that real(nout/2) was supposed to be real(nout)/2.
- ## sqrtHalfNfreq <- sqrt(nfreq %/% 2)
- sqrtHalfNfreq <- sqrt(nfreq / 2)
- ## dplR: NOTE! Is round() the right function to use? Maybe floor()
- ## for the lower limit and ceiling for the higher limit?
- rcritlo <- round((-0.79557086 + 1.0088719 * sqrtHalfNfreq)^2)
- rcrithi <- round(( 0.75751462 + 0.9955133 * sqrtHalfNfreq)^2)
-
if (csv.out || do.table) {
dframe <- c(x[c("freq", "gxx", "gxxc", "gredth", "grravg", "corr")],
list(gredth * fac80, gredth * fac90,
@@ -449,6 +443,16 @@
}
if (!csv.out) {
## dplR: print miscellaneous information AND if (do.table) print(dframe)
+ nseg <- params[["nseg"]]
+ ofac <- params[["ofac"]]
+ rhopre <- params[["rhopre"]]
+
+ ## critical false alarm level after Thomson (1990)
+ ## dplR: modified from original REDFIT code to accommodate for
+ ## lower / upper tail difference
+ alphacrit <- (nseg - 1) / nseg
+ faccrit <- qchisq(alphacrit, dof) / dof
+
precat("redfit()", newline = FALSE)
vers <- x[["vers"]]
if (!is.null(vers)) {
@@ -513,13 +517,16 @@
domain = "R-dplR")
precat(gtxt)
precat(rep.int("-", nchar(gtxt)))
- if (iwin == 0 && ofac == 1 && n50 == 1) {
+ rcnt <- x[["rcnt"]]
+ if (!is.null(rcnt)) {
gtxt <- gettext("5-% acceptance region:", domain = "R-dplR")
precat(gtxt, newline = FALSE)
- cat(" rcritlo = ", format(rcritlo, digits = digits), "\n", sep = "")
+ cat(" rcritlo = ", format(x[["rcritlo"]], digits = digits), "\n",
+ sep = "")
precat(rep.int(" ", nchar(gtxt)), newline = FALSE)
- cat(" rcrithi = ", format(rcrithi, digits = digits), "\n", sep = "")
- precat("r_test = ", format(x[["rcnt"]], digits = digits))
+ cat(" rcrithi = ", format(x[["rcrithi"]], digits = digits), "\n",
+ sep = "")
+ precat("r_test = ", format(rcnt, digits = digits))
} else {
if (iwin != 0) {
precat(gettext("Test requires iwin = 0", domain = "R-dplR"))
Modified: branches/redfit/man/redfit.Rd
===================================================================
--- branches/redfit/man/redfit.Rd 2013-09-03 15:13:48 UTC (rev 678)
+++ branches/redfit/man/redfit.Rd 2013-09-04 08:00:54 UTC (rev 679)
@@ -110,8 +110,22 @@
\item{rcnt }{ a \code{numeric} value giving the number of runs in a
statistical test studying the difference between a theoretical AR1
- spectrum and the bias-corrected spectrum estimated from the data. }
+ spectrum and the bias-corrected spectrum estimated from the data.
+ Requires that \code{\var{iwin} == 0} (\code{"rectangular"}),
+ \code{\var{ofac} == 1} and \code{\var{n50} == 1}. Otherwise the
+ test is not performed and \var{rcnt} is \code{NULL}. See Bendat and
+ Piersol, p. 95. }
+ \item{rcritlo }{ a \code{numeric} critical low value for
+ \code{\var{rcnt}}. Approximately 2.5 percent of the null
+ distribution lies below this value (TO BE CHECKED). Is \code{NULL}
+ when \code{\var{rcnt}} is \code{NULL}. }
+
+ \item{rcritlo }{ a \code{numeric} critical high value for
+ \code{\var{rcnt}}. Approximately 2.5 percent of the null
+ distribution lies above this value (TO BE CHECKED). Is \code{NULL}
+ when \code{\var{rcnt}} is \code{NULL}.}
+
\item{freq }{ the frequencies used. A \code{numeric} vector. The
other numeric vectors have the same length, i.e. one value for each
frequency studied. }
@@ -179,6 +193,8 @@
\item{nsim }{ value of the \code{\var{nsim}} argument. }
\item{mctest }{ value of the \code{\var{mctest}} argument. }
+
+ \item{rhopre }{ value of the \code{\var{rhopre}} argument. }
}
}
@@ -194,10 +210,13 @@
\href{http://www.ncdc.noaa.gov/paleo/softlib/redfit/redfit.html}{REDFIT},
which is in the public domain.
+ Bendat, J. S. and Piersol, A. G. (1986) \emph{Random Data: Analysis
+ and Measurement Procedures.} Wiley. \acronym{ISBN}: 0-471-04000-2.
+
Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating red-noise
spectra directly from unevenly spaced paleoclimatic time series.
\emph{Computers & Geosciences}, 28(3):421\enc{–}{--}426.
-
+
}
\author{
Mikko Korpela
More information about the Dplr-commits
mailing list