Nothing
#
# This file contains code related to patch size distribution fitting. These
# functions can fit Power-law (pl), Truncated Power-law (tpl), Lognormal
# (lnorm) and Exponential (exp) distributions using maximum likelihood, as
# per Clauset et al. 's (2007) recommandations.
#
# In addition, it provides the estimation of xmin using the ks-distance for
# power-laws (same, following Clauset et al. 2007)
# Optimisation global options
ITERLIM <- 10000
GRADTOL <- 1e-10
STEPTOL <- 1e-10
STEPMAX <- 5
ifnotfinite <- function(x, otherwise = .Machine$double.xmax) {
ifelse(is.finite(x), x, sign(x) * otherwise)
}
warnbound <- function() {
getOption("spatialwarnings.debug.fit_warn_on_bound", default = FALSE)
}
warnNA <- function() {
getOption("spatialwarnings.debug.fit_warn_on_NA", default = FALSE)
}
# This is a safe version of nlm that returns a sensible result (NaNs) when
# the algorithm fails to converge. This can happen quite often when looking
# for pathological cases (e.g. fitting distribution based on few points in the
# tails, etc.).
optim_safe <- function(f, pars0,
lower = rep(-Inf, length(pars0)),
upper = rep(Inf, length(pars0)),
fit_on_logscale = FALSE, ...) {
if ( fit_on_logscale ) {
optimf <- function(pars) {
f(exp(pars))
}
pars0 <- log(pars0)
lower <- suppressWarnings(ifnotfinite(log(lower)))
upper <- suppressWarnings(ifnotfinite(log(upper)))
} else {
optimf <- f
}
optiresult <- try({
optim(pars0, optimf,
control = list(maxit = ITERLIM),
lower = lower, upper = upper,
method = "L-BFGS-B", ...)
}, silent = TRUE)
# Code results above 3 means a true problem, below 3 the
# solution is either exact or approximate.
# Sometimes L-BFGS-B gets stuck in a very flat area just because our initial
# guess was very good, and L-BFGS-B does not like that (it reports
# abnormal termination of line search). Or L-BFGS-B will try the bounds of
# the parameter space, which will return error. So here we try to use BFGS
# which should report success if our initial guess was not too bad.
if ( inherits(optiresult, "try-error") || optiresult[["convergence"]] > 3 ) {
optiresult_bfgs <- try({
optim(pars0, optimf,
control = list(maxit = ITERLIM),
method = "BFGS", ...)
}, silent = TRUE)
# If success, go with BFGS results
if ( inherits(optiresult_bfgs, "try-error") ) {
optiresult <- optiresult_bfgs
}
}
# If we could not reach a proper solution, report an error
if ( inherits(optiresult, "try-error") ) {
optiresult <- list(value = NaN,
par = rep(NaN, length(pars0)),
convergence = 128)
return(optiresult)
}
if ( optiresult[["convergence"]] > 3 ) {
warning(paste0('optim returned an error (error code:',
optiresult[["convergence"]], ").\n",
"Make sure the results are reasonable using plot_distr"))
}
# Convert the estimated pars back from log scale
if ( fit_on_logscale ) {
optiresult[["par"]] <- exp(optiresult[["par"]])
}
return(optiresult)
}
# Bounds on parameters, these should be large and no observed distribution should
# have values beyond them
# Power-laws lambdas
PLMIN <- 1 + sqrt(.Machine$double.eps)
PLMAX <- 10
# Exponential rates
EXPMIN <- sqrt(.Machine$double.eps) # A very close value to, but not, zero
EXPMAX <- 10
# Bounds for truncated power-laws
TPL_EXPOMIN <- -1 # Taken from Clauset's code
TPL_EXPOMAX <- 10
TPL_RATEMIN <- sqrt(.Machine$double.eps)
TPL_RATEMAX <- 10
# Riemann zeta function with xmin taken into account :
# sum( 1/k^-expo ) for i=xmin to i = inf
# This is vectorized over xmins so that we do not sum things several times.
zeta_w_xmin <- function(expo, xmins) {
perm <- order(xmins)
xmins <- xmins[order(xmins)]
# Compute zeta value
zetaval <- gsl::zeta(expo)
# Initialize
output <- rep(NaN, length(xmins))
current_k <- xmins[1]
output[perm[1]] <- zetaval - sum_all_one_over_k(from = 1, to = xmins[1], expo)
# If there is only one value, we bail now
if ( length(xmins) <= 1) {
return(output)
}
for ( i in 2:length(xmins)) {
next_k <- xmins[i]
if (next_k > current_k) {
output[perm[i]] <- output[perm[i-1]] -
sum_all_one_over_k(from = current_k,
to = next_k, expo)
current_k <- next_k
} else {
output[perm[i]] <- output[perm[i-1]]
}
}
return(output)
}
# PL fitting
# ---------------------------------------
# Normalizing constant for pl with xmin
displnorm <- function(expo, xmin) {
# Adjust constant for threshold (note that this has no effect if xmin == 1,
# as expected)
const <- gsl::zeta(expo)
const - sum_all_one_over_k(from = 1, to = xmin, expo)
}
# PL: P(x=k)
dpl <- function(x, expo, xmin = 1, log = FALSE) {
const <- displnorm(expo, xmin)
# Compute values
if ( ! log ) {
ans <- (1/const) * x^(-expo)
ans[x < xmin] <- NaN
} else {
if ( const < 0 ) {
# Const can be negative as nlm finds its way: the check makes sure
# no warning is produced by the log.
ans <- NaN
} else {
ans <- -expo * log(x) - log(const)
}
}
return(ans)
}
# PL: P(x>=k)
ippl <- function(x, expo, xmin = 1) {
const <- displnorm(expo, xmin)
is_below_xmin <- x < xmin
ps <- zeta_w_xmin(expo, x[!is_below_xmin]) / const
# Values below threshold are NA'ed
ans <- NaN*x
ans[!is_below_xmin] <- ps
return(ans)
}
# PL: Log likelihood
pl_ll <- function(dat, expo, xmin) {
sum( dpl(dat, expo, xmin, log = TRUE) )
}
#' @title Distribution-fitting functions
#'
#' @description These functions fit parametric distributions to a set of
#' discrete values.
#'
#' @param dat The set of values to which the distribution are fit
#'
#' @param xmin The minimum possible value to consider when fitting the
#' distribution
#'
#' @return A list containing at list the following components:
#'
#' \describe{
#' \item{\code{type}}{The type of distribution fitted (as a character string)}
#' \item{\code{method}}{The method used for the fit - here, maximum likelihood, 'll'}
#' \item{\code{ll}}{The log likelihood at the estimated parameter values}
#' \item{\code{xmin}}{The value of xmin used for the fit}
#' \item{\code{npars}}{The number of parameters of the distribution}
#' }
#'
#' Additionally, this list may have one or more of the following elements depending on
#' the type of distribution that has been fitted:
#' \describe{
#' \item{\code{plexpo}}{The exponent of the power-law}
#' \item{\code{cutoff}}{The rate of truncation, for truncated power law and
#' exponential fits}
#' \item{\code{meanlog}}{The mean of the lognormal distribution}
#' \item{\code{sdlog}}{The s.d. of the lognormal distribution}
#' }
#'
#' @details These functions will fit distributions to a set of values using
#' maximum-likelihood estimation. In the context of the 'spatialwarnings'
#' package, they are most-often used to fit parametric distributions on patch
#' size distributions. As a result, these functions assume that the data
#' contains only integer, strictly positive values. The type of distribution
#' depends on the prefix of the function: 'pl' for power-law, 'tpl' for
#' truncated power-law, 'lnorm' for lognormal and 'exp' for an exponential
#' distribution.
#'
#' In the context of distribution-fitting, 'xmin' represents the minimum value
#' that a distribution can take. It is often used to represent the minimum
#' scale at which a power-law model is appropriate (Clauset et al. 2009), and
#' can be estimated on an empirical distribution using
#' \code{\link{xmin_estim}}. Again, please note that the fitting procedure
#' assumes here that xmin is equal or grater than one.
#'
#' Please note that a best effort is made to have the fit converge, but
#' it may sometimes fail when the parameters are far from their usual
#' range, and numerical issues may occur. It is good practice to make
#' sure the fits are sensible when convergence warnings are reported.
#'
#' For reference, the shape of the distributions is as follow:
#'
#' \describe{
#' \item{power-law}{\eqn{x^{-a}}{x^(-a)} where a is the power-law exponent}
#' \item{exponential}{\eqn{exp(-bx)}{exp(-bx)} where b is the truncation rate
#' of the exponential}
#' \item{truncated power-law}{\eqn{x^{-a}exp(-bx)}{x^(-a)exp(-bx)} where a
#' and b are the exponent of the power law and the rate of truncation}
#' }
#'
#' The lognormal form follows the \link[=dlnorm]{standard definition}.
#'
#' The following global options can be used to change the behavior of fitting functions
#' and/or produce more verbose output:
#' \describe{
#' \item{spatialwarnings.constants.reltol}{the relative tolerance to use to compute
#' the power-law normalizing constant
#' \deqn{sum_{k=1}^{\infty} x^{ak}e^{-bk}}{sum( x^(ak)exp(-bk)) for k in 1:Inf}.
#' Increase to increase the precision of this constant, which can be useful in some
#' cases, typically with large sample sizes. Default is 1e-8.}
#' \item{spatialwarnings.constants.maxit}{the maximum number of iterations to compute
#' the normalizing constant of a truncated power-law. Increase if you get a warning
#' that the relative tolerance level (defined above) was not reached. Default is 1e8}
#' \item{spatialwarnings.debug.fit_warn_on_bound}{logical value. Warn if the fit is
#' at the boundary of the valid range for distribution parameter}
#' \item{spatialwarnings.debug.fit_warn_on_NA}{logical value. Warn if the returned fit
#' has \code{NA}/\code{NaN} parameters}
#' }
#'
#' @seealso \code{\link{patchdistr_sews}}, \code{\link{xmin_estim}}
#'
#' @references
#'
#' Clauset, Aaron, Cosma Rohilla Shalizi, and M. E. J. Newman. 2009. “Power-Law
#' Distributions in Empirical Data.” SIAM Review 51 (4): 661–703.
#' https://doi.org/10.1137/070710111.
#'
#' @examples
#'
#' # Fit an exponential model to patch size distribution
#' exp_fit(patchsizes(forestgap[[8]]))
#'
#' # Use the estimated parameters as an indicator function
#' \donttest{
#'
#' get_truncation <- function(mat) {
#' c(exp_cutoff = exp_fit(patchsizes(mat))$cutoff)
#' }
#' trunc_indic <- compute_indicator(forestgap, get_truncation)
#' plot(trunc_indic)
#' plot(indictest(trunc_indic, nulln = 19))
#'
#' }
#'
#'@export
pl_fit <- function(dat, xmin = 1) {
# Cut data to specified range
dat <- dat[dat >= xmin]
# Start with the approximation given in Clauset's
npts <- length(dat)
expo_estim <- 1 + npts / (sum(log(dat)) - npts*log(xmin-.5))
negll <- function(expo) {
result <- - pl_ll(dat, expo, xmin)
if ( is.infinite(result) ) {
return(NaN)
} else {
return(result)
}
}
est <- optim_safe(negll, expo_estim,
lower = PLMIN, upper = PLMAX)
result <- list(type = 'pl',
method = 'll',
plexpo = est[["par"]],
ll = - est[['value']],
xmin = xmin,
npars = 1)
if ( warnNA() && is.na(result[["plexpo"]]) ) {
warning("Fitting of PL returned NA")
return(result)
}
if ( warnbound() && any(abs(result[["plexpo"]] - c(PLMIN, PLMAX)) < 1e-8) ) {
warning("Estimated PL exponent is on the boundary of the valid range")
}
return(result)
}
#' @title Estimate the minimum patch size of a power-law distribution
#'
#' @description When fitting a power-law to a discrete distribution, it might
#' be worth discarding points below a certain threshold (xmin) to improve
#' the fit. This function estimates the optimal xmin based on the
#' Kolmogorov-Smirnoff distance between the fit and the empirical
#' distribution, as suggested by Clauset et al. (2009).
#'
#' @param dat A vector of integer values
#'
#' @param bounds A vector of two values representing the bounds in which
#' the best xmin is searched
#'
#' @return The estimated xmin as an integer value
#'
#' @details The function returns NA if \code{dat} has only three unique values
#' or if the power-law fit failed.
#'
#' @seealso \code{\link{patchdistr_sews}}
#'
#' @references
#'
#' Clauset, A., Shalizi, C. R., & Newman, M. E. (2009).
#' Power-law distributions in empirical data. SIAM review, 51(4), 661-703.
#'
#' @examples
#'
#' \donttest{
#' psd <- patchsizes(forestgap[[5]])
#' xmin_estim(psd)
#' }
#'@export
xmin_estim <- function(dat, bounds = range(dat)) {
# Create a vector of possible values for xmin
xmins <- sort(unique(dat))
# We need at least 3 values for a pl fit, so the last value of xmin
# needs to have three points after it
if ( length(xmins) <= 3 ) {
return(NaN)
}
# We build a vector of possible xmins. The last three values are stripped
# away as they won't allow enough data for a fit
xmins <- head(xmins, length(xmins)-3)
xmins <- xmins[xmins >= min(bounds) & xmins <= max(bounds)]
# Compute all ks-distances
kss <- adply(xmins, 1, get_ks_dist, dat = dat)[ ,2]
if ( all(is.nan(kss)) ) {
return(NaN)
}
# Note that sometimes the fit fails, especially when xmin is around the
# distribution tail -> we need to remove some NAs here
xmin <- xmins[!is.na(kss) & kss == min(kss, na.rm = TRUE)]
# Note that xmin can be NaN
return(xmin)
}
get_ks_dist <- function(xmin, dat) {
# Crop dat to values above xmin and compute cdf
dat <- dat[dat >= xmin]
# Compute empirical (inverse) cdf values
udat <- unique(dat)
cdf_empirical <- rep(NA, length(dat))
for ( val in udat ) {
cdf_empirical[dat == val] <- mean(dat >= val)
}
# Fit and retrieve cdf. Note: here we suppress the warnings because finding
# xmin requires removing patches below a threshold, which often leads to fit
# being done on pathological cases like few unique patch sizes
fit <- suppressWarnings({
pl_fit(dat, xmin = xmin)
})
if ( is.na(fit[['plexpo']]) ) {
# Note: a warning was already produced in this case as it means that the
# fit failed to converge: we do not produce one here again.
return(NaN)
}
cdf_fitted <- ippl(dat, fit[["plexpo"]], fit[["xmin"]])
# # debug
# plot(data.frame(dat, rbinom(length(dat), 1, .5)), type = 'n')
# plot(log10(data.frame(dat, cdf_empirical)))
# points(log10(data.frame(dat, cdf_fitted)), col = 'red')
# browser()
# zeta.fit(dat, xmin)$exponent
# fit$expo
# We return the ks distance
maxks <- max(abs(cdf_empirical - cdf_fitted))
# cat(xmin, ",", max(dat), "->", maxks, "\n" )
return( maxks )
}
# EXP fitting
# ---------------------------------------
# pexp/dexp is already implemented in R
# EXP: P(x = k)
ddisexp <- function(dat, rate, xmin = 1, log = FALSE) {
# sum(P = k) for k = 1 to inf
if ( log ) {
const <- log(1 - exp(-rate)) + rate * xmin
return( ifelse(dat < xmin, NaN, const - rate * dat) )
} else {
const <- (1 - exp(-rate)) * exp(rate*xmin)
return( ifelse(dat < xmin, NaN, const * exp(-rate * dat)) )
}
}
# EXP: P(x>=k)
# Imported and cleaned up from powerRlaw (def_disexp.R)
ipdisexp <- function(x, rate, xmin) {
# p >= k
p <- pexp(x + .5, rate, lower.tail = FALSE)
# p >= 1
const <- 1 - pexp(xmin + .5, rate)
return(p/const)
}
exp_ll <- function(dat, rate, xmin) {
sum( ddisexp(dat, rate, xmin, log = TRUE))
}
#'@rdname pl_fit
#'@export
exp_fit <- function(dat, xmin = 1) {
dat <- dat[dat>=xmin]
rate0 <- 1 / mean(dat)
negll <- function(rate) {
- exp_ll(dat, rate, xmin)
}
est <- optim_safe(negll, rate0,
fit_on_logscale = TRUE,
lower = EXPMIN,
upper = EXPMAX)
result <- list(type = 'exp',
method = 'll',
cutoff = est[['par']],
ll = - est[["value"]],
npars = 1)
if ( warnbound() && any(abs(result[["cutoff"]] - c(EXPMIN, EXPMAX)) < 1e-8) ) {
warning("Estimated EXP exponent is on the boundary of the valid range")
}
return(result)
}
# LNORM fitting
# ---------------------------------------
# LNORM: P(X=k)
ddislnorm <- function(x, meanlog, sdlog, xmin, log = FALSE) {
p_over_thresh <- plnorm(xmin - .5, meanlog, sdlog, lower.tail = FALSE)
p_equals_k <- plnorm(x-.5, meanlog, sdlog, lower.tail = FALSE) -
plnorm(x+.5, meanlog, sdlog, lower.tail = FALSE)
if ( !log ) {
return( ifelse(x<xmin, NaN, p_equals_k / p_over_thresh) )
} else {
return( ifelse(x<xmin, NaN, log(p_equals_k) - log(p_over_thresh)) )
}
}
# LNORM: P(X>=k)
ipdislnorm <- function(x, meanlog, sdlog, xmin) {
px_supto_k <- plnorm(x - .5, meanlog, sdlog, lower.tail = FALSE)
px_supto_xmin <- plnorm(xmin - .5, meanlog, sdlog, lower.tail = FALSE)
ifelse(x<xmin, NaN, px_supto_k / px_supto_xmin)
}
# LNORM: LL
lnorm_ll <- function(x, meanlog, sdlog, xmin) {
x <- x[x>=xmin]
sum( ddislnorm(x, meanlog, sdlog, xmin, log = TRUE) )
}
# LNORM: fit
#'@rdname pl_fit
#'@export
lnorm_fit <- function(dat, xmin = 1) {
# Pars[1] holds mean of log-transformed data
# Pars[2] holds sd
pars0 <- c( mean(log(dat)), sd(log(dat)) )
negll <- function(pars) {
ll <- - lnorm_ll(dat, pars[1], pars[2], xmin)
if ( is.finite(ll) ) ll else 1e10
}
est <- optim_safe(negll, pars0)
result <- list(type = 'lnorm',
method = 'll',
meanlog = est[['par']][1],
sdlog = est[['par']][2],
ll = - est[["value"]],
npars = 2)
return(result)
}
# TPL fitting
# ---------------------------------------
tplnorm <- function(expo, rate, xmin) {
maxit <- getOption("spatialwarnings.constants.maxit", default = 1e8L)
reltol <- getOption("spatialwarnings.constants.reltol", default = 1e-8)
a <- tplinfsum(expo, rate, xmin, maxit, reltol)
a
}
# P(x=k)
dtpl <- function(x, expo, rate, xmin, log = FALSE) {
const <- tplnorm(expo, rate, xmin)
if ( ! log ) {
ps <- x^(-expo) * exp(- x * rate) / const
} else {
ps <- - expo * log(x) - rate * x - log(const)
}
return( ifelse(x < xmin, NaN, ps) )
}
# P(x>=k)
iptpl <- function(x, expo, rate, xmin) {
const <- tplnorm(expo, rate, xmin)
# tplsum is vectorized over x
p_inf_to_k <- tplsum(expo, rate, x, xmin) / const
return( 1 - p_inf_to_k )
}
tpl_ll <- function(x, expo, rate, xmin, approximate = FALSE) {
x <- x[x>=xmin]
ll <- sum( dtpl(x, expo, rate, xmin, log = TRUE) )
if ( !is.finite(ll) ) {
ll <- sign(ll) * .Machine$double.xmax
}
return( ll )
}
#'@rdname pl_fit
#'@export
tpl_fit <- function(dat, xmin = 1) {
negll <- function(pars) {
- tpl_ll(dat, pars[1], pars[2], xmin)
}
# Initialize and find minimum
expo0 <- pl_fit(dat, xmin)[['plexpo']]
# Do a line search over the cutoff to find a minimum, starting from zero
# up to 100
is <- seq(0, 100, length = 128)
is <- 10^seq(-7, 2, l = 128)
lls <- unlist(lapply(is, function(i) {
negll(c(expo0, i))
}))
llmin <- min(lls[abs(lls) != .Machine$double.xmax & is.finite(lls)])
expmrate0 <- is[which(lls == llmin)]
# If multiple cutoffs produce the same ll, then use the lowest one
if ( length(expmrate0) > 1 ) {
expmrate0 <- expmrate0[1]
}
# Debug thing lls
# plot( log(spatialwarnings:::cumpsd(dat)) )
# vals <- ippl(dat, expo0, xmin = xmin)
# points(log(dat), log(vals), col = "red")
# vals <- iptpl(dat, expo0, rate = 0, xmin = xmin)
# points(log(dat), log(vals), col = "darkgreen")
# vals <- iptpl(dat, expo0, rate = expmrate0, xmin = xmin)
# points(log(dat), log(vals), col = "blue")
pars0 <- c(expo0, expmrate0)
est <- optim_safe(negll, pars0,
lower = c(TPL_EXPOMIN, TPL_RATEMIN),
upper = c(TPL_EXPOMAX, TPL_RATEMAX),
fit_on_logscale = FALSE)
# For very small cutoffs, we may want to fit on log scale to get a descent estimate,
# as otherwise there is too little variation in the cutoff for BFGS
if ( any(is.na(est[["par"]])) ) {
est <- optim_safe(negll, pars0,
lower = c(TPL_EXPOMIN, TPL_RATEMIN),
upper = c(TPL_EXPOMAX, TPL_RATEMAX),
fit_on_logscale = TRUE)
}
result <- list(type = 'tpl',
method = 'll',
plexpo = est[['par']][1],
cutoff = est[['par']][2],
ll = - est[["value"]],
npars = 2)
if ( warnNA() && is.na(result[["plexpo"]]) ) {
warning("Fitting of TPL returned NA/NaN")
}
if ( is.na(result[["plexpo"]]) ) {
return(result)
}
if ( warnbound() &&
any(abs(result[["plexpo"]] - c(TPL_EXPOMIN, TPL_EXPOMAX)) < 1e-8) ) {
warning("Estimated TPL exponent is on the boundary of the valid range")
}
# We do not warn for reaching minimum exponent because this is quite common when
# fitting truncated power laws to have a very small exponent
if ( warnbound() &&
any(abs(result[["cutoff"]] - c(TPL_RATEMIN, TPL_RATEMAX)) < 1e-8) ) {
warning("Estimated TPL cutoff is on the boundary of the valid range")
}
return(result)
}
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.