[Dplr-commits] r717 - in pkg/dplR: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 27 11:28:43 CET 2013
Author: mvkorpel
Date: 2013-11-27 11:28:42 +0100 (Wed, 27 Nov 2013)
New Revision: 717
Modified:
pkg/dplR/ChangeLog
pkg/dplR/DESCRIPTION
pkg/dplR/R/detrend.R
pkg/dplR/R/detrend.series.R
pkg/dplR/man/detrend.Rd
pkg/dplR/man/detrend.series.Rd
pkg/dplR/man/dplR-package.Rd
Log:
Fixed a bug in detrend.series() where RWI could go negative.
Thanks to Jacob Cecile (added to DESCRIPTION as a contributor).
Also, detrend() and detrend.series() have a new argument 'constrain.modnegexp'.
Modified: pkg/dplR/ChangeLog
===================================================================
--- pkg/dplR/ChangeLog 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/ChangeLog 2013-11-27 10:28:42 UTC (rev 717)
@@ -1,5 +1,24 @@
* CHANGES IN dplR VERSION 1.5.8
+File: DESCRIPTION
+
+- Added Jacob Cecile as a contributor
+
+File: detrend.R
+---------------
+
+- Adjusted for the changes in detrend.series(),
+ i.e. added the argument constrain.modnegexp.
+
+File: detrend.series.R
+----------------------
+
+- Fixed a bug where RWI could go negative. Thanks to Jacob Cecile for
+ reporting the bug and contributing a proposed solution.
+- A new argument: constrain.modnegexp. It is now possible to constrain
+ the modified negative exponential function to non-negative values at
+ infinity.
+
File: redfit.R
--------------
Modified: pkg/dplR/DESCRIPTION
===================================================================
--- pkg/dplR/DESCRIPTION 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/DESCRIPTION 2013-11-27 10:28:42 UTC (rev 717)
@@ -3,7 +3,7 @@
Type: Package
Title: Dendrochronology Program Library in R
Version: 1.5.8
-Date: 2013-11-09
+Date: 2013-11-27
Authors at R: c(person("Andy", "Bunn", role = c("aut", "cph",
"cre", "trl"), email = "andy.bunn at wwu.edu"), person("Mikko",
"Korpela", role = c("aut", "trl")), person("Franco", "Biondi",
@@ -12,8 +12,9 @@
"cph")), person("Manfred", "Mudelsee", role = "aut"),
person("Fares", "Qeadan", role = c("aut", "cph")),
person("Michael", "Schulz", role = "aut"), person("Christian",
- "Zang", role = c("aut", "cph")))
-Author: Andy Bunn [aut, cph, cre, trl], Mikko Korpela [aut, trl], Franco Biondi [aut, cph], Filipe Campelo [aut, cph], Pierre Mérian [aut, cph], Manfred Mudelsee [aut], Fares Qeadan [aut, cph], Michael Schulz [aut], Christian Zang [aut, cph]
+ "Zang", role = c("aut", "cph")), person("Jacob", "Cecile",
+ role = "ctb"))
+Author: Andy Bunn [aut, cph, cre, trl], Mikko Korpela [aut, trl], Franco Biondi [aut, cph], Filipe Campelo [aut, cph], Pierre Mérian [aut, cph], Manfred Mudelsee [aut], Fares Qeadan [aut, cph], Michael Schulz [aut], Christian Zang [aut, cph], Jacob Cecile [ctb]
Copyright: Authors and Aalto University (for work of M. Korpela)
Maintainer: Andy Bunn <andy.bunn at wwu.edu>
Depends: R (>= 2.15.0)
Modified: pkg/dplR/R/detrend.R
===================================================================
--- pkg/dplR/R/detrend.R 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/R/detrend.R 2013-11-27 10:28:42 UTC (rev 717)
@@ -1,9 +1,13 @@
`detrend` <-
function(rwl, y.name = names(rwl), make.plot = FALSE,
method=c("Spline", "ModNegExp", "Mean"),
- nyrs = NULL, f = 0.5, pos.slope = FALSE)
+ nyrs = NULL, f = 0.5, pos.slope = FALSE,
+ constrain.modnegexp = c("never", "when.fail", "always"))
{
+ stopifnot(identical(make.plot, TRUE) || identical(make.plot, FALSE),
+ identical(pos.slope, FALSE) || identical(pos.slope, TRUE))
known.methods <- c("Spline", "ModNegExp", "Mean")
+ constrain2 <- match.arg(constrain.modnegexp)
method2 <- match.arg(arg = method,
choices = known.methods,
several.ok = TRUE)
@@ -36,7 +40,9 @@
fits <- detrend.series(rwl.i, make.plot=FALSE,
method=method2,
nyrs=nyrs, f=f,
- pos.slope=pos.slope)
+ pos.slope=pos.slope,
+ constrain.modnegexp=
+ constrain2)
if(is.data.frame(fits))
row.names(fits) <- rn
fits
@@ -47,7 +53,8 @@
fits <- detrend.series(rwl[[i]], y.name=y.name[i],
make.plot=make.plot,
method=method2, nyrs=nyrs, f=f,
- pos.slope=pos.slope)
+ pos.slope=pos.slope,
+ constrain.modnegexp=constrain2)
if(is.data.frame(fits))
row.names(fits) <- rn
out[[i]] <- fits
Modified: pkg/dplR/R/detrend.series.R
===================================================================
--- pkg/dplR/R/detrend.series.R 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/R/detrend.series.R 2013-11-27 10:28:42 UTC (rev 717)
@@ -1,11 +1,13 @@
`detrend.series` <-
function(y, y.name = "", make.plot = TRUE,
method = c("Spline", "ModNegExp", "Mean"),
- nyrs = NULL, f = 0.5, pos.slope = FALSE)
+ nyrs = NULL, f = 0.5, pos.slope = FALSE,
+ constrain.modnegexp = c("never", "when.fail", "always"))
{
stopifnot(identical(make.plot, TRUE) || identical(make.plot, FALSE),
identical(pos.slope, FALSE) || identical(pos.slope, TRUE))
known.methods <- c("Spline", "ModNegExp", "Mean")
+ constrain2 <- match.arg(constrain.modnegexp)
method2 <- match.arg(arg = method,
choices = known.methods,
several.ok = TRUE)
@@ -24,27 +26,54 @@
if("ModNegExp" %in% method2){
## Nec or lm
- nec.func <- function(Y) {
+ nec.func <- function(Y, constrain) {
a <- mean(Y[seq_len(floor(length(Y) * 0.1))])
b <- -0.01
k <- mean(Y[floor(length(Y) * 0.9):length(Y)])
- nec <- nls(formula = Y ~ a * exp(b * seq_along(Y)) + k,
- start = list(a=a, b=b, k=k))
- if(coef(nec)[2] >= 0) stop()
- fits <- predict(nec)
- if(fits[1] < fits[length(fits)]) stop()
- if(fits[length(fits)] < 0) stop()
+ nlsForm <- Y ~ a * exp(b * seq_along(Y)) + k
+ nlsStart <- list(a=a, b=b, k=k)
+ checked <- FALSE
+ if (constrain == "never") {
+ nec <- nls(formula = nlsForm, start = nlsStart)
+ } else if (constrain == "always") {
+ nec <- nls(formula = nlsForm, start = nlsStart,
+ lower = c(a=0, b=-Inf, k=0),
+ upper = c(a=Inf, b=0, k=Inf),
+ algorithm = "port")
+ } else {
+ nec <- nls(formula = nlsForm, start = nlsStart)
+ if(coef(nec)[2] >= 0) stop()
+ fits <- predict(nec)
+ if(fits[1] < fits[length(fits)]) stop()
+ if(fits[length(fits)] > 0) {
+ checked <- TRUE
+ } else {
+ nec <- nls(formula = nlsForm, start = nlsStart,
+ lower = c(a=0, b=-Inf, k=0),
+ upper = c(a=Inf, b=0, k=Inf),
+ algorithm = "port")
+ }
+ }
+ if (!checked) {
+ if(coef(nec)[2] >= 0) stop()
+ fits <- predict(nec)
+ if(fits[1] < fits[length(fits)]) stop()
+ if(fits[length(fits)] <= 0) stop()
+ }
fits
}
- ModNegExp <- try(nec.func(y2), silent=TRUE)
+ ModNegExp <- try(nec.func(y2, constrain2), silent=TRUE)
if(class(ModNegExp)=="try-error") {
## Straight line via linear regression
tm <- cbind(1, seq_along(y2))
lm1 <- lm.fit(tm, y2)
coefs <- lm1[["coefficients"]]
+ ModNegExp <- NULL
if (all(is.finite(coefs)) && (coefs[2] <= 0 || pos.slope)) {
ModNegExp <- drop(tm %*% coefs)
- } else {
+ }
+ if (is.null(ModNegExp) ||
+ ModNegExp[1] <= 0 || ModNegExp[length(y2)] <= 0) {
ModNegExp <- rep(mean(y2), length(y2))
}
}
@@ -64,6 +93,9 @@
else
nyrs2 <- nyrs
Spline <- ffcsaps(y=y2, x=seq_along(y2), nyrs=nyrs2, f=f)
+ if (any(Spline <= 0)) {
+ Spline <- rep(mean(y2), length(y2))
+ }
resids$Spline <- y2 / Spline
do.spline <- TRUE
} else {
Modified: pkg/dplR/man/detrend.Rd
===================================================================
--- pkg/dplR/man/detrend.Rd 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/man/detrend.Rd 2013-11-27 10:28:42 UTC (rev 717)
@@ -8,7 +8,8 @@
\usage{
detrend(rwl, y.name = names(rwl), make.plot = FALSE,
method = c("Spline", "ModNegExp", "Mean"), nyrs = NULL,
- f = 0.5, pos.slope = FALSE)
+ f = 0.5, pos.slope = FALSE,
+ constrain.modnegexp = c("never", "when.fail", "always"))
}
\arguments{
@@ -38,6 +39,9 @@
slope to be used in method \code{"ModNegExp"}. If \code{FALSE} the
line will be horizontal. }
+ \item{constrain.modnegexp}{ a \code{character} string which controls
+ the constraints of the \code{"ModNegExp"} model. }
+
}
\details{
See \code{\link{detrend.series}} for details on detrending
Modified: pkg/dplR/man/detrend.series.Rd
===================================================================
--- pkg/dplR/man/detrend.series.Rd 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/man/detrend.series.Rd 2013-11-27 10:28:42 UTC (rev 717)
@@ -8,7 +8,8 @@
\usage{
detrend.series(y, y.name = "", make.plot = TRUE,
method = c("Spline", "ModNegExp", "Mean"),
- nyrs = NULL, f = 0.5, pos.slope = FALSE)
+ nyrs = NULL, f = 0.5, pos.slope = FALSE,
+ constrain.modnegexp = c("never", "when.fail", "always"))
}
\arguments{
@@ -35,6 +36,16 @@
\item{pos.slope}{ a \code{logical} flag. Will allow for a positive
slope to be used in method \code{"ModNegExp"}. If \code{FALSE} the
line will be horizontal. }
+
+ \item{constrain.modnegexp}{ a \code{character} string which controls
+ the constraints of the \code{"ModNegExp"} model. The value is an
+ answer to the question: When should the parameters of the modified
+ negative exponential function be constrained? The options are
+ \code{"never"}: do not constrain (the default), \code{"when.fail"}:
+ only compute the constrained solution if the unconstrained fit
+ contains other than positive values, and \code{"always"}: return the
+ constrained solution, even if the unconstrained one would have been
+ valid. See \sQuote{Details}. }
}
\details{
@@ -57,16 +68,25 @@
stand effects.
The \code{"ModNegExp"} approach attempts to fit a classic nonlinear
- model of biological growth of the form \code{\var{Y} ~ \var{a} *
- exp(\var{b} * 1:length(\var{Y})) + \var{k}} using
- \code{\link{nls}}. See Fritts (2001) for details about the
- parameters. If a nonlinear model cannot be fit then a linear model is
- fit. That linear model can have a positive slope unless
- \code{\var{pos.slope}} is \code{FALSE} in which case method
+ model of biological growth of the form \eqn{f(t) = a e^{b t} + k}{f(t)
+ = a exp(b t) + k}, where the argument of the function is time, using
+ \code{\link{nls}}. See Fritts (2001) for details about the
+ parameters. Option \code{\var{constrain.modnegexp}} gives a
+ possibility to constrain the parameters of the modified negative
+ exponential function. If the constraints are enabled, the nonlinear
+ optimization algorithm is instructed to keep the parameters in the
+ following ranges: \eqn{a \ge 0}{a >= 0}, \eqn{b \le 0}{b <= 0} and
+ \eqn{k \ge 0}{k >= 0}. If a suitable nonlinear model cannot be fit
+ (function is non-decreasing or some values are not positive) then a
+ linear model is fit. That linear model can have a positive slope
+ unless \code{\var{pos.slope}} is \code{FALSE} in which case method
\code{"Mean"} is used.
The \code{"Mean"} approach fits a horizontal line using the mean of
- the series.
+ the series. This method is the fallback solution in cases where the
+ \code{"Spline"} or the linear fit (also a fallback solution itself)
+ contains zeros or negative values, which would lead to invalid
+ ring-width indices.
These methods are chosen because they are commonly used in
dendrochronology. It is, of course, up to the user to determine the
@@ -86,7 +106,8 @@
Fritts, H. C. (2001) \emph{Tree Rings and Climate}.
Blackburn. \acronym{ISBN-13}: 978-1-930665-39-2.
}
-\author{ Andy Bunn. Patched and improved by Mikko Korpela. }
+\author{ Andy Bunn. Patched and improved by Mikko Korpela. A bug fix
+ related to negative output values is based on work by Jacob Cecile. }
\seealso{ \code{\link{detrend}} }
\examples{library(stats)
## Using a plausible representation of a tree-ring series
Modified: pkg/dplR/man/dplR-package.Rd
===================================================================
--- pkg/dplR/man/dplR-package.Rd 2013-11-09 21:59:24 UTC (rev 716)
+++ pkg/dplR/man/dplR-package.Rd 2013-11-27 10:28:42 UTC (rev 717)
@@ -29,7 +29,8 @@
Campelo, Pierre \enc{Mérian}{Merian}, Fares Qeadan and Christian
Zang. Function \code{\link{redfit}} is an improved translation of
program REDFIT which is original work of Manfred Mudelsee and Michael
- Schulz.
+ Schulz. Jacob Cecile contributed a bug fix to
+ \code{\link{detrend.series}}.
}
\references{
Cook, E. R. and Kairiukstis, L. A. (1990) \emph{Methods of
More information about the Dplr-commits
mailing list