logcondiscr (version 1.0.6)

logConDiscrCI: Compute pointwise confidence bands for the log-concave MLE of a PMF

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

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.

Value

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)$\}$.

Details

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

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

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