Nothing
#' Fine fitting functions
#'
#' double logistics, piecewise logistics and many other functions to
#' curve fit VI time-series.
#'
#' * `Logistic` The traditional simplest logistic function. It can
#' be only used in half growing season, i.e. vegetation green-up or senescence
#' period.
#' * `doubleLog.Zhang` Piecewise logistics, (Zhang Xiaoyang, RSE, 2003).
#' * `doubleAG` Asymmetric Gaussian.
#' * `doubleLog.Beck` Beck logistics.
#' * `doubleLog.Gu` Gu logistics.
#' * `doubleLog.Elmore` Elmore logistics.
#' * `doubleLog.Klos` Klos logistics.
#'
#' All of those function have `par` and `formula` attributes for the
#' convenience for analytical D1 and D2
#'
#' @param par A vector of parameters
#' @param t A `Date` or numeric vector
#'
#' @references
#' 1. Beck, P.S.A., Atzberger, C., Hogda, K.A., Johansen, B., Skidmore, A.K.,
#' 2006. Improved monitoring of vegetation dynamics at very high latitudes:
#' A new method using MODIS NDVI. Remote Sens. Environ.
#' https://doi.org/10.1016/j.rse.2005.10.021.
#' 2. Elmore, A.J., Guinn, S.M., Minsley, B.J., Richardson, A.D., 2012.
#' Landscape controls on the timing of spring, autumn, and growing season
#' length in mid-Atlantic forests. Glob. Chang. Biol. 18, 656-674.
#' https://doi.org/10.1111/j.1365-2486.2011.02521.x. \cr
#'
#' 3. Gu, L., Post, W.M., Baldocchi, D.D., Black, TRUE.A., Suyker, A.E., Verma,
#' S.B., Vesala, TRUE., Wofsy, S.C., 2009. Characterizing the Seasonal Dynamics
#' of Plant Community Photosynthesis Across a Range of Vegetation Types,
#' in: Noormets, A. (Ed.), Phenology of Ecosystem Processes: Applications
#' in Global Change Research. Springer New York, New York, NY, pp. 35-58.
#' https://doi.org/10.1007/978-1-4419-0026-5_2. \cr
#' 4. Peter M. Atkinson, et al., 2012, RSE, 123:400-417
#'
#' 5. https://github.com/cran/phenopix/blob/master/R/FitDoubleLogGu.R
#' @example R/examples/ex-FitDL.R
#' @rdname logistics
#' @export
Logistic = function(par, t){
mn = par[1]
mx = par[2]
sos = par[3]
rsp = par[4]
pred = (mx - mn)/(1 + exp(-rsp*(t - sos))) + mn
# pred = c/(1 + exp(a + b * t)) + d
return(pred)
}
attr(Logistic, 'par') = c("mn", "mx", "sos", "rsp")
attr(Logistic, 'formula') = expression((mx - mn)/(1 + exp(-rsp*(t - sos))) + mn)
# piecewise function
#' @rdname logistics
#' @export
doubleLog.Zhang = function(par, t){
t0 = par[1]
mn = par[2]
mx = par[3]
sos = par[4]
rsp = par[5]
eos = par[6]
rau = par[7]
if (t0 - sos <= 1 || t0 - eos >= -1) return(rep(99.0, length(t)))
# In order to make sure the shape of S curve, should be satisfy:
# t0 < eos, t0 > sos
# xpred1 = (mx - mn)/(1 + exp(-rsp*(t - sos))) + mn
# xpred2 = (mx - mn)/(1 + exp(rau*(t - eos))) + mn
# pred = xpred1*(t <= t0) + xpred2*(t > t0)
# above expressions cost 1.5 times as the below
pred = (mx - mn)/c(1 + exp(-rsp*(t[t <= t0] - sos)),
1 + exp( rau*(t[t > t0] - eos))) + mn
return(pred)
}
attr(doubleLog.Zhang, 'par') = c("t0", "mn", "mx", "sos", "rsp", "eos", "rau")
attr(doubleLog.Zhang, 'formula') = expression( mn + (mx - mn)/(1 + exp(-rsp*(t - sos))),
mn + (mx - mn)/(1 + exp( rau*(t - eos))) )
#' @rdname logistics
#' @export
doubleLog.AG = function(par, t){
t0 = par[1]
mn = par[2]
mx = par[3]
rsp = par[4]
a3 = par[5]
rau = par[6]
a5 = par[7]
pred = mn + (mx - mn)*exp(- c(
((t0 - t[t <= t0])*rsp) ^a3,
((t[t > t0] - t0)*rau) ^a5) )
return(pred)
}
# a3, a5 should be greater than 1
attr(doubleLog.AG, 'par') = c("t0", "mn", "mx", "rsp", "a3", "rau", "a5")
attr(doubleLog.AG, 'formula') = expression( mn + (mx - mn)*exp(- ((t0 - t)*rsp) ^a3 ),
mn + (mx - mn)*exp(- ((t - t0)*rau) ^a5 ))
#' @rdname logistics
#' @export
doubleLog.AG2 = function(par, t){
t0 = par[1]
mn_l = par[2]
mn_r = par[3]
mx = par[4]
rsp = par[5]
a3 = par[6]
rau = par[7]
a5 = par[8]
t_l = t[t <= t0]
t_r = t[t > t0]
pred_l = mn_l + (mx - mn_l)*exp(- ((t0 - t_l)*rsp) ^a3)
pred_r = mn_r + (mx - mn_r)*exp(- ((t_r - t0)*rau) ^a5)
# pred = mn + (mx - mn)*exp(- c( ((t0 - t[t <= t0])*rsp) ^a3,
# ((t[t > t0] - t0)*rau) ^a5) )
# return(pred)
c(pred_l, pred_r)
}
# a3, a5 should be greater than 1
attr(doubleLog.AG2, 'par') = c("t0", "mn_l", "mn_r", "mx", "rsp", "a3", "rau", "a5")
attr(doubleLog.AG2, 'formula') = expression(
mn_l + (mx - mn_l)*exp(- ((t0 - t)*rsp) ^a3),
mn_r + (mx - mn_r)*exp(- ((t - t0)*rau) ^a5)
)
#' @rdname logistics
#' @export
doubleLog.Beck = function(par, t) {
mn = par[1]
mx = par[2]
sos = par[3]
rsp = par[4]
eos = par[5]
rau = par[6]
# if (sos >= eos) return(rep(9999, length(t)))
# if (eos < sos) return(rep(99.0, length(t)))
if (!all(is.finite(par)) || eos < sos) return(rep(99.0, length(t)))
pred = mn + (mx - mn)*(1/(1 + exp(-rsp*(t - sos))) + 1/(1 + exp(rau*(t - eos))) - 1)
return(pred)
}
attr(doubleLog.Beck, 'par') = c("mn", "mx", "sos", "rsp", "eos", "rau")
attr(doubleLog.Beck, "formula") = expression(mn + (mx - mn) * (1 / (1 + exp(-rsp * (t - sos))) + 1 / (1 + exp(rau * (t - eos))) - 1))
#' @rdname logistics
#' @export
doubleLog.Elmore = function(par, t) {
mn = par[1]
mx = par[2]
# m3 = par[3]
# m4 = par[4]
# m5 = par[5]
# m6 = par[6]
# m3l = m3/m4
# m4l = 1/m4
# m5l = m5/m6
# m6l = 1/m6
sos = par[3] # SOS
rsp = par[4] # 1/rsp
eos = par[5] # EOS
rau = par[6] # 1/rau
m7 = par[7]
pred = mn + (mx - m7*t)*( 1/(1 + exp(-rsp*(t-sos))) - 1/(1 + exp(-rau*(t-eos))) )
return(pred)
}
attr(doubleLog.Elmore, 'par') = c("mn", "mx", "sos", "rsp", "eos", "rau", "m7")
attr(doubleLog.Elmore, 'formula') = expression( mn + (mx - m7*t)*( 1/(1 + exp(-rsp*(t-sos))) - 1/(1 + exp(-rau*(t-eos))) ) )
#' @rdname logistics
#' @export
doubleLog.Gu = function(par, t) {
y0 = par[1]
a1 = par[2]
a2 = par[3]
sos = par[4]
rsp = par[5]
eos = par[6]
rau = par[7]
# t1 = par[4]
# t2 = par[5]
# b1 = par[6]
# b2 = par[7]
c1 = par[8]
c2 = par[9]
# pred = y0 + (a1/(1 + exp(-(t - t1)/b1))^c1) - (a2/(1 + exp(-(t - t2)/b2))^c2)
pred = y0 + (a1/(1 + exp(-rsp*(t - sos)))^c1) - (a2/(1 + exp(-rau*(t - eos)))^c2)
return(pred)
}
attr(doubleLog.Gu, 'par') = c('y0', 'a1', 'a2', 'sos', 'rsp', 'eos', 'rau', 'c1', 'c2')
attr(doubleLog.Gu, 'formula') = expression(y0 + (a1/(1 + exp(-rsp*(t - sos)))^c1) - (a2/(1 + exp(-rau*(t - eos)))^c2))
#' @rdname logistics
#' @export
doubleLog.Klos = function(par, t) {
a1 = par[1]
a2 = par[2]
b1 = par[3]
b2 = par[4]
c = par[5]
B1 = par[6]
B2 = par[7]
m1 = par[8]
m2 = par[9]
q1 = par[10]
q2 = par[11]
v1 = par[12]
v2 = par[13]
pred = (a1*t + b1) + (a2*t^2 + b2*t + c) *
(1/(1 + q1 * exp(-B1 * (t - m1)))^v1 -
1/(1 + q2 * exp(-B2 * (t - m2)))^v2 )
return(pred)
}
attr(doubleLog.Klos, 'par') = c('a1', 'a2', 'b1', 'b2', 'c', 'B1', 'B2',
'm1', 'm2', 'q1', 'q2', 'v1', 'v2')
attr(doubleLog.Klos, 'formula') = expression((a1*t + b1) + (a2*t^2 + b2*t + c) * (1/(1 + q1 * exp(-B1 * (t - m1)))^v1
- 1/(1 + q2 * exp(-B2 * (t - m2)))^v2))
.qr.solve = function(a, b, tol = 1e-07, LAPACK = TRUE) {
if (!is.qr(a)) a = qr(a, tol = tol, LAPACK = LAPACK)
nc = ncol(a$qr)
nr = nrow(a$qr)
if (a$rank != min(nc, nr)) stop("singular matrix 'a' in solve")
if (missing(b)) {
if (nc != nr) stop("only square matrices can be inverted")
b = diag(1, nc)
}
res = qr.coef(a, b)
res[is.na(res)] = 0
res
}
# vc = .qr.solve(opt$hessian)
# npar = nrow(vc)
# s2 = opt.df$cost[best]^2 / (n - npar)
# std.errors = sqrt(diag(vc) * s2) # standard errors
# attach gradient and hessian analytical function to curve fitting functions
.dls = lapply(
c("doubleLog.Beck", "doubleLog.Elmore", "doubleLog.Gu",
"doubleLog.Klos", "doubleLog.Zhang", "doubleLog.AG"),
function (FUN){
# FUN = deparse(substitute(fun))
fun = get(FUN)
attr(fun, 'gradient') = gradf_t(fun) # gradient
attr(fun, 'hessian') = hessf_t(fun) # hessian
# print(environment(fun))
assign(FUN, fun, envir = environment(fun)) #environment("namespace:phenofit"))#
# fun
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.