logConDiscrCI: Compute pointwise confidence bands for the log-concave MLE of... In logcondiscr: Estimate a Log-Concave Probability Mass Function from Discrete i.i.d. Observations

Description

Compute pointwise confidence bands for the log-concave maximum likelihood estimate of a log-concave probability mass function based on the limiting theory developed in Balabdaoui et al (2012).

Usage

 1 2 logConDiscrCI(dat, conf.level = 0.95, type = c("MLE", "all")[1], B = 1000, output = TRUE, seed = 2011)

Arguments

 dat Data to compute MLE and confidence band for. conf.level The confidence level to be used. type To compute confidence bands one theoretically needs to know the knots of the true PMF. For type MLE the knots of the MLE will be used instead and for type all all observations will be considered knots. The latter is conservative and gives pointwise confidence intervals that are based on standard errors from a Normal approximation (the latter comes from the asymptotic theory in Balabdaoui et al, 2011). B Number of samples to be drawn to compute resampling quantiles. output If TRUE, progress of computations is output. seed Optional seed.

Details

The pointwise confidence bands are based on the limiting theory in Balabdaoui et al (2011).

Value

A list with the following components:

 MLE The estimated MLE (simply the output list of the function logConDiscrMLE applied to dat). emp A dataframe containing two columns: the unique sorted observations and the empirical PMF. CIs The computed confidence intervals for each x \in \{min(dat), ..., max(dat)\}.

Note

Values outside [0, 1] will be clipped. As a consequence, coverage may be higher than 1 - α.

Author(s)

Kaspar Rufibach (maintainer) [email protected]
http://www.kasparrufibach.ch
Fadoua Balabdaoui [email protected]
Hanna Jankowski [email protected]
http://www.math.yorku.ca/~hkj
Kathrin Weyermann

References

Balabdaoui, F., Jankowski, H., Rufibach, K., and Pavlides, M. (2013). Maximum likelihood estimation and confidence bands for a discrete log-concave distribution. J. R. Stat. Soc. Ser. B Stat. Methodol., 75(4), 769–790.

Weyermann, K. (2007). An Active Set Algorithm for Log-Concave Discrete Distributions. MSc thesis, University of Bern (Supervisor: Lutz Duembgen).

Examples

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 # ------------------------------------------------------------- # compute MLE and confidence bands for a random sample from a # Poisson distribution # ------------------------------------------------------------- set.seed(2011) x <- sort(rpois(n = 100, lambda = 2)) mle <- logConDiscrMLE(x) psi <- mle$psi # compute confidence bands CIs <- logConDiscrCI(x, type = "MLE", output = TRUE, seed = 20062011)$CIs # plot estimated PMF and log of estimate par(mfrow = c(1, 2), las = 1) true <- dpois(0:20, lambda = 2) plot(mle$x, exp(psi), type = "b", col = 2, xlim = c(min(x), max(x) + 1), xlab = "x", ylim = c(0, max(exp(psi), true, CIs[, 3])), ylab = "PMF", main = "Estimate MLE from a Poisson", pch = 19) legend("topright", c("truth", "MLE", "confidence bands"), col = c(4, 2, 2), lty = c(1, 1, 2), pch = c(0, 19, NA), bty = "n") # add true PMF lines(0:20, true, type = "l", col = 4) # add confidence bands matlines(CIs[, 1], CIs[, 2:3], type = "l", lty = 2, col = 2) # log-density plot(mle$x, psi, type = "p", col = 2, xlim = c(min(x), max(x) + 1), xlab = "x", ylab = "PMF", main = "Estimate MLE from a Poisson", pch = 19, ylim = c(-6, log(max(exp(psi), true, CIs[, 3])))) lines(0:20, log(true), type = "l", col = 4) # add confidence bands matlines(CIs[, 1], log(CIs[, 2:3]), type = "l", lty = 2, col = 2) # ------------------------------------------------------------- # compute confidence bands when only estimate (not original data) # are available (as a an example we simply use the estimator from # above) # ------------------------------------------------------------- x.est <- 0:6 est <- c(0.09, 0.30, 0.30, 0.19, 0.09, 0.02, 0.01) # generate original data (up to given precision) x <- rep(0:6, times = 100 * est)

logcondiscr documentation built on May 30, 2017, 7:11 a.m.