# tests/testthat/pli-R-v0.0.3-2007-07-25/weibull.R In spatialwarnings: Spatial Early Warning Signals of Ecosystem Degradation

```# Functions for estimation of a stretched exponential or Weibull distribution

# The use of a Weibull distribution to model values above a specified lower
# threshold can be done in one of two ways.  The first is simply to shift the
# the standard Weibull, i.e, to say that x-threshold ~ weibull.
# The other is to say that Pr(X|X>threshold) = weibull(X|X>threshold), i.e.,
# that the right tail follows the same functional form as the right tail of a
# Weibull, without necessarily having the same probability of being in the tail.
# These will be called the "shift" and "tail" methods respectively.
# The shift method is, computationally, infinitely easier, but not so suitable
# for our purposes.

# The basic R system provides dweibull (density), pweibull (cumulative
# distribution), qweibull (quantiles) and rweibull (random variate generation)

### Function for fitting:
# weibull.fit			Fit Weibull to data with choice of methods
### Distributional functions, per R standards:
# dweibull.tail			Tail-conditional probability density
# pweibull.tail			Tail-conditional cumulative distribution
### Backstage functions, not for users:
# weibull.fit.tail		Fit by maximizing tail-conditional likelihood
#				(default)
# weibull.fit.eqns		Fit by solving ML estimating equations for
#				shifted Weibull
# weibull.est.shape.inefficient Inefficient estimator of shape for shifted
#				Weibull, used to initialize fit.eqns
# weibull.est.scale.from.shape  MLE of scale given shape for shifted Weibull
# weibull.loglike.shift		Log-likelihood of shifted Weibull
# weibull.loglike.tail		Tail-conditional log-likelihood

# Fit Weibull to data
# Wrapper for functions implementing different ML methods
# Input: Data vector, lower threshold, method flag
# Output: List giving distribution type ("weibull"), parameters, log-likelihood
weibull.fit <- function(x, threshold = 0,method="tail") {
x <- x[x>=threshold]
switch(method,
tail = {
# Estimate parameters by direct maxmization of the tail-conditional
# log-likelihood
fit <- weibull.fit.tail(x,threshold)
},
eqns = {
# Estimate parameters by solving the ML estimating equations of a shifted
# Weibull
fit <- weibull.fit.eqns(x,threshold)
},
{
cat("Unknown method\n")
fit <- NA})
return(fit)
}

# Tail-conditional probability density function
# Returns NA on values < threshold
# Input: Data vector, distributional parameters, log flag
# Output: Vector of (log) probability densities
dweibull.tail <- function(x,shape,scale,threshold=0,log=FALSE) {
c <- pweibull(threshold,shape=shape,scale=scale,lower.tail=FALSE,log.p=log)
if (log) {
f <- function(y) {dweibull(y,shape,scale,log=TRUE) - c}
} else {
f <- function(y) {dweibull(y,shape,scale)/c}
}
d <- ifelse(x<threshold,NA,f(x))
return(d)
}

# Tail-conditional cumulative distribution function
# Returns NA on values < threshold
# Input: Data vector, distributional parameters, usual flags
# Output: Vector of (log) cumulative probabilities
pweibull.tail <- function(x,shape,scale,threshold=0,lower.tail=TRUE,
log.p=FALSE) {
c <- pweibull(threshold,shape,scale,lower.tail=FALSE)
c.log <- pweibull(threshold,shape,scale,lower.tail=FALSE,log.p=TRUE)
if ((!lower.tail)&&(!log.p)) {
f <- function(x) {pweibull(x,shape,scale,lower.tail=FALSE)/c}
}
if ((!lower.tail)&&(log.p)) {
f <- function(x) {pweibull(x,shape,scale,lower.tail=FALSE,log.p=TRUE) - c.log}
}
if ((lower.tail)&&(!log.p)) {
f <- function(x) {1 - pweibull(x,shape,scale,lower.tail=FALSE)/c}
}
if ((lower.tail)&&(log.p)) {
f <- function(x) {log(1 - pweibull(x,shape,scale,lower.tail=FALSE)/c)}
}
p <- ifelse(x<threshold,NA,f(x))
return(p)
}

# Fit Weibull to data by maximizing tail-conditional likelihood
# CONSTRAINTS: The shape and scale parameters must both be positive
# This will not necessarily give a _stretched_ exponential which would be
# a shape < 1
# Input: Data vector, lower threshold
# Output: List giving distribution type ("weibull"), parameters, log-likelihood
weibull.fit.tail <- function(x,threshold=0) {
# Use the whole data to produce initial estimates, via the estimating equation
initial_fit <- weibull.fit.eqns(x)
theta_0 <- c(initial_fit\$shape, initial_fit\$scale)
# Apply constraints: if we're outside the feasible set, default to a
# standardized exponential
if (theta_0 < 0) { theta_0 = 1 }
if (theta_0 < 0) { theta_0 = 1 }
# Now threshold and directly minimize the negative log likelihood
x <- x[x>= threshold]
n <- length(x)
negloglike <- function(theta) {
-weibull.loglike.tail(x,threshold,shape=theta,scale=theta)
}
ui <- rbind(c(1,0),c(0,1))
ci <- c(0,0)
est <- constrOptim(theta=theta_0, f=negloglike, grad=NULL, ui=ui, ci=ci)
fit <- list(type="weibull", shape=est\$par, scale=est\$par,
loglike = -est\$value, samples.over.threshold=n)
return(fit)
}

# Fit shifted Weibull to data by solving ML estimating equations
# Input: Data vector, lower threshold
# Output: List giving distribution type ("weibull"), parameters, log-likelihood
weibull.fit.eqns <- function(x, threshold = 0) {
# ML estimating equations of the shape and scale of the Weibull distribution,
# taken from Johnson and Kotz, ch. 20
# This assumes that x-threshold is Weibull-distributed
x <- x[x>=threshold]
x <- x-threshold
x.log <- log(x) # This will be needed repeatedly
n<-length(x) # So will this
h <- sum(x.log)/n # And this
scale_from_shape <- function(shape) {weibull.est.scale.from.shape(x,shape)}
# Note that the estimation of the shape parameter is only implicit,
# through a transcendental equation.
initial_estimates <- weibull.est.shape.inefficient(x)
shape <- initial_estimates
scale <- initial_estimates
map <- function(c) {(sum((x^c) * x.log)/sum(x^c) - h)^(-1)}
estimating_equation <- function(c) { (c - map(c))^2 }
shape <- nlm(f=estimating_equation,p=shape)\$estimate
scale <- scale_from_shape(shape)
loglike <- weibull.loglike.shift(x,threshold,shape,scale)
fit <- list(type="weibull", shape=shape, scale=scale, loglike=loglike, samples.over.threshold=n)
return(fit)
}

# Log-likelihood of shifted Weibull distribution
# Input: Data vector, parameters, lower threshold
# Output: real-valued log-likelihood
weibull.loglike.shift <- function(x, shape, scale, threshold = 0) {
# Assumes x - threshold is Weibull-distributed
# See Johnson and Kotz for more; they call the lower threshold xi_0
x <- x[x>=threshold]
x <- x-threshold
L <- sum(dweibull(x,shape,scale,log=TRUE))
return(L)
}

# Tail-conditional log-likelihood
# Input: Data vector, parameters, lower threshold
# Output: Real-valued log-likelihood
weibull.loglike.tail <- function(x, shape, scale, threshold = 0) {
# We want p(X=x|X>=threshold) = p(X=x)/Pr(X>=threshold)
x <- x[x>=threshold]
n <- length(x)
Like <- weibull.loglike.shift(x,shape,scale)
ThresholdProb <- pweibull(threshold,shape,scale,log=TRUE,lower.tail = FALSE)
L <- Like - n*ThresholdProb
return(L)
}

# Inefficient estimator of the shape parameter of a Weibull, plus scale
# Based on the moments of log of the data
# Input: Data vector, lower threshold
# Output: Real-valued estimates of shape and scale
weibull.est.shape.inefficient <- function(x,threshold=0) {
# The follow is not an efficient estimator of the shape, but serves to start
# the approximation process off.  (See Johnson and Katz, ch. 20, section 6.3,
# "Estimators Based on Distribution of log X".)
x <- x-threshold
shape <- (sqrt(6)/pi)*sd(log(x))
scale <- weibull.est.scale.from.shape(x,shape)
c(shape,scale)
}

# Maximum likelihood estimate of a Weibull scale parameter, given shape
# Input: Data vector, shape parameter
# Output: Real-valued scale parameter
weibull.est.scale.from.shape <- function(x,shape) {
# Given a value of the shape parameter, return the corresponding
# MLE of the scale parameter
n <- length(x)
scale <- ((1/n)*sum(x^shape))^(1/shape)
return(scale)
}
```

## Try the spatialwarnings package in your browser

Any scripts or data that you put into this service are public.

spatialwarnings documentation built on May 2, 2019, 5:16 p.m.