eco-hydro / phenofit

R package: A state-of-the-art Vegetation Phenology extraction package, phenofit
http://phenofit.top
GNU General Public License v2.0
73 stars 32 forks source link

Beck doublelogistic parametrization typo #11

Closed SveBuch closed 2 years ago

SveBuch commented 2 years ago

Hello,

I'm currently trying to fit a double logistic curve to Sentinel derived NDVI values. When predicting values, based on the fitted curve, the curve was always biased low. It seems as if there is a typo in the formula returned by curvefit. I tracked the error down to the function doubleLog.Beck.

In there, the double-log curve is parametrized as mn + (mx - mn) * (1/(1 + exp(-rsp * (t - sos))) + 1/(1 + exp(rau * (t - eos)))) - 1 while in Beck et al. 2006 (https://doi.org/10.1016/j.rse.2005.10.021) it is parametrized as mn + (mx - mn) * (1/(1 + exp(-rsp * (t - sos))) + 1/(1 + exp(rau * (t - eos))) - 1)

What I don't understand, is, that if I predict with the corrected formula but using the parameters fitted with the one with the typo, the resulting prediction looks quite okay. Maybe you can help me there?

All the best, Sven

Here is a reproducible example:

## doublelog parametrization as in output from curvefit
dlogfun <- function(par, t) {
  mn <- par$mn; mx <- par$mx; sos <- par$sos; rsp <- par$rsp; eos <- par$eos; rau <- par$rau
  mn + (mx - mn) * (1/(1 + exp(-rsp * (t - sos))) + 1/(1 + exp(rau * (t - eos)))) - 1
}

## doublelog parametrization as in Beck et al.
dlogfun2 <- function(par, t) {
  mn <- par$mn; mx <- par$mx; sos <- par$sos; rsp <- par$rsp; eos <- par$eos; rau <- par$rau
  mn + (mx - mn) * (1/(1 + exp(-rsp * (t - sos))) + 1/(1 + exp(rau * (t - eos))) - 1)
}

## modified example from phenofit help
par  = c( mn  = 0.1, mx  = 0.7, sos = 50, rsp = 0.1, eos = 250, rau = 0.1)
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
y <- doubleLog.Beck(par, t)  + rnorm(length(t), 0, 0.1)
fFITs <- curvefit(y, t, tout, "Beck")

## double log-parametrisation from curvefit
fFITs$model

(l_param   <- get_param(list(fFITs))$Beck)

ggplot() +
  geom_point(aes(x = t, y = y)) +
  geom_line(aes(x = tout, y = dlogfun(l_param, tout)), col="red") +
  geom_line(aes(x = tout, y = dlogfun2(l_param, tout))) +
  theme_bw()
kongdd commented 2 years ago

The equation is right in https://github.com/eco-hydro/phenofit/blob/master/R/doubleLogistics_R.R#L148. Could you tell which line do you find the error equation?

kongdd commented 2 years ago

The output is normal:

library(phenofit)
library(ggplot2)

## modified example from phenofit help
par  = c( mn  = 0.1, mx  = 0.7, sos = 50, rsp = 0.1, eos = 250, rau = 0.1)
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
y <- doubleLog.Beck(par, t)  + rnorm(length(t), 0, 0.1)

fit <- curvefit(y, t, tout, "Beck")

## double log-parametrisation from curvefit
par_opt = fit$model$Beck$par
# l_par <- get_param(fit)$Beck # v0.3.5 install from github

ypred = doubleLog.Beck(l_par %>% as.numeric(), tout)

ggplot() +
    geom_point(aes(x = t, y = y)) +
    geom_line(aes(x = tout, y = ypred), col="red") +
    # geom_line(aes(x = tout, y = dlogfun2(l_param, tout))) +
    theme_bw()

image

SveBuch commented 2 years ago

Ah, I see. The typo is in line 152 https://github.com/eco-hydro/phenofit/blob/master/R/doubleLogistics_R.R#L152, meaning only the attribute of the function is affected. So the fitting is done correctly, but since I used the attribute in a formula, the curve was shifted. Thank you

kongdd commented 2 years ago

Thanks for your report. This bug will influence phenology extraction. I will fix it.