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

```# Poisson distribution for the right/uppper tail, and its estimation

# R by default provides dpois (probability mass function), ppois (cumulative
# probability function), qpois (quantile function) and rpois (random variate
# generation)

### Function for fitting:
# pois.tail.fit		Fit Poisson to tail of data by numerically maximizing
#			likelihood
### Distributional function, per R standard:
# dpois.tail		Probability mass function, conditional on being in tail
### Backstage function, not intended for users:
# pois.tail.loglike	Calcualte log-likelihood

# Probability mass function, conditional on being in tail
# Outputs "NA" on data points below cut-off
# Input: Data vector, Poisson rate, lower cut-off, log-flag
# Output: Vector of (log) probabilities
dpois.tail <- function(x, rate, threshold=0, log=FALSE) {
C <- ppois(threshold-1,rate,lower.tail=FALSE,log.p=log)
if (log) {
f <- function(x) {-C + dpois(x,rate,log=TRUE)}
} else {
f <- function(x) {dpois(x,rate)/C}
}
d <- ifelse(x<threshold, NA, f(x))
return(d)
}

# Calculate tail-conditional log likelihood
# Input: Data vector, Poisson rate, lower cut-off
# Output: Log likelihood
pois.tail.loglike <- function(x, rate, threshold=0) {
n <- length(x)
JointProb <- sum(dpois(x,rate,log=TRUE))
ProbOverThreshold <- ppois(threshold-1,rate,lower.tail=FALSE,log.p=TRUE)
return(JointProb - n*ProbOverThreshold)
}

# Fit Poisson distribution to right/upper tail of data, by numerically
# maximizing log likelihood
# Input: Data vector, lower cut-off
# Output: List giving distribution type, estimated rate, information about
#	  fit and data
pois.tail.fit <- function(x, threshold=0) {
# Do a simple moment-based fit first on the whole data as a pump-primer
rate_0 <- mean(x)
# discard data points below threshold
x <- x[x>=threshold]
rate_1 <- mean(x)
n <- length(x)
negloglike <- function(rate) { -pois.tail.loglike(x, rate, threshold) }
est <- nlm(f=negloglike, p=rate_0)
fit <- list(type="pois.tail", rate = est\$estimate, threshold=threshold,
loglike=-est\$minimum, samples.over.threshold=n, full.mean=rate_0,
mean.over.threshold=rate_1)
return(fit)
}
```

## 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.