mable: Maximum Approximate Bernstein Likelihood Estimate of...
In mable: Maximum Approximate Bernstein/Beta Likelihood Estimation

mable

R Documentation

Maximum Approximate Bernstein Likelihood Estimate of Univariate or Multivariate Density Function

Description

Maximum approximate Bernstein/Beta likelihood estimation based on one-sample raw data with an optimal selected by the change-point method among m0:m1 or a preselected model degree m.

Usage

mable(
  x,
  M,
  interval = 0:1,
  IC = c("none", "aic", "hqic", "all"),
  vb = 0,
  controls = mable.ctrl(),
  progress = TRUE
)

Arguments

`x`	a (non-empty) numeric `n`-vector or `n x d` matrix or `data.frame` of `n` observations.
`M`	a positive integer or a vector `(m0, m1)`. If `M = m` or `m0 = m1 = m`, then `m` is a preselected degree. If `m0<m1` it specifies the set of consective candidate model degrees `m0:m1` for searching an optimal degree, where `m1-m0>3`.
`interval`	a vector containing the endpoints of supporting/truncation interval `c(a,b)`
`IC`	information criterion(s) in addition to Bayesian information criterion (BIC). Current choices are "aic" (Akaike information criterion) and/or "qhic" (Hannan–Quinn information criterion).
`vb`	code for vanishing boundary constraints, -1: f0(a)=0 only, 1: f0(b)=0 only, 2: both, 0: none (default).
`controls`	Object of class `mable.ctrl()` specifying iteration limit and the convergence criterion `eps`. Default is `mable.ctrl`. See Details.
`progress`	if TRUE a text progressbar is displayed

Details

Any continuous density function f on a known closed supporting interval [a,b] can be estimated by Bernstein polynomial f_m(x; p) = \sum_{i=0}^m p_i\beta_{mi}[(x-a)/(b-a)]/(b-a), where p = (p_0, \ldots, p_m), p_i \ge 0, \sum_{i=0}^m p_i = 1 and \beta_{mi}(u) = (m+1){m\choose i}u^i(1-x)^{m-i}, i = 0, 1, \ldots, m, is the beta density with shapes (i+1, m-i+1). For each m, the MABLE of the coefficients p, the mixture proportions, are obtained using EM algorithm. The EM iteration for each candidate m stops if either the total absolute change of the log likelihood and the coefficients of Bernstein polynomial is smaller than eps or the maximum number of iterations maxit is reached.

If m0<m1, an optimal model degree is selected as the change-point of the increments of log-likelihood, log likelihood ratios, for m \in \{m_0, m_0+1, \ldots, m_1\}. Alternatively, one can choose an optimal degree based on the BIC (Schwarz, 1978) which are evaluated at m \in \{m_0, m_0+1, \ldots, m_1\}. The search for optimal degree m is stoped if either m1 is reached with a warning or the test for change-point results in a p-value pval smaller than sig.level. The BIC for a given degree m is calculated as in Schwarz (1978) where the dimension of the model is d = \#\{i: \hat p_i\ge\epsilon, i = 0, \ldots, m\} - 1 and a default \epsilon is chosen as .Machine$double.eps.

If data show a clearly multimodal distribution by plotting the histogram for example, the model degree is usually large. The range M should be large enough to cover the optimal degree and the computation is time-consuming. In this case the iterative method of moment with an initial selected by a method of mode which is implemented by optimable can be used to reduce the computation time.

Value

A list with components

m the given or a selected degree by method of change-point
p the estimated vector of mixture proportions p = (p_0, \ldots, p_m) with the selected/given optimal degree m
mloglik the maximum log-likelihood at degree m
interval support/truncation interval (a,b)
convergence An integer code. 0 indicates successful completion (all the EM iterations are convergent and an optimal degree is successfully selected in M). Possible error codes are
- 1, indicates that the iteration limit maxit had been reached in at least one EM iteration;
- 2, the search did not finish before m1.
delta the convergence criterion delta value

and, if m0<m1,

M the vector (m0, m1), where m1, if greater than m0, is the largest candidate when the search stoped
lk log-likelihoods evaluated at m \in \{m_0, \ldots, m_1\}
lr likelihood ratios for change-points evaluated at m \in \{m_0+1, \ldots, m_1\}
ic a list containing the selected information criterion(s)
pval the p-values of the change-point tests for choosing optimal model degree
chpts the change-points chosen with the given candidate model degrees

Note

Since the Bernstein polynomial model of degree m is nested in the model of degree m+1, the maximum likelihood is increasing in m. The change-point method is used to choose an optimal degree m. The degree can also be chosen by a method of moment and a method of mode which are implemented by function optimal().

Author(s)

Zhong Guan <zguan@iu.edu>

References

Guan, Z. (2016) Efficient and robust density estimation using Bernstein type polynomials. Journal of Nonparametric Statistics, 28(2):250-271. Wang, T. and Guan, Z.,(2019) Bernstein Polynomial Model for Nonparametric Multivariate Density, Statistics, Vol. 53, no. 2, 321-338

Examples


# Vaal Rive Flow Data
 data(Vaal.Flow)
 x<-Vaal.Flow$Flow
 res<-mable(x, M = c(2,100), interval = c(0, 3000), controls =
        mable.ctrl(sig.level = 1e-8, maxit = 2000, eps = 1.0e-9))
 op<-par(mfrow = c(1,2),lwd = 2)
 layout(rbind(c(1, 2), c(3, 3)))
 plot(res, which = "likelihood", cex = .5)
 plot(res, which = c("change-point"), lgd.x = "topright")
 hist(x, prob = TRUE, xlim = c(0,3000), ylim = c(0,.0022), breaks = 100*(0:30),
  main = "Histogram and Densities of the Annual Flow of Vaal River",
  border = "dark grey",lwd = 1,xlab = "x", ylab = "f(x)", col  = "light grey")
 lines(density(x, bw = "nrd0", adjust = 1), lty = 4, col = 4)
 lines(y<-seq(0, 3000, length = 100), dlnorm(y, mean(log(x)),
                   sqrt(var(log(x)))), lty = 2, col = 2)
 plot(res, which = "density", add = TRUE)
 legend("top", lty = c(1, 2, 4), col = c(1, 2, 4), bty = "n",
 c(expression(paste("MABLE: ",hat(f)[B])),
        expression(paste("Log-Normal: ",hat(f)[P])),
               expression(paste("KDE: ",hat(f)[K]))))
 par(op)


# Old Faithful Data
 library(mixtools)
 x<-faithful$eruptions
 a<-0; b<-7
 v<-seq(a, b,len = 512)
 mu<-c(2,4.5); sig<-c(1,1)
 pmix<-normalmixEM(x,.5, mu, sig)
 lam<-pmix$lambda; mu<-pmix$mu; sig<-pmix$sigma
 y1<-lam[1]*dnorm(v,mu[1], sig[1])+lam[2]*dnorm(v, mu[2], sig[2])
 res<-mable(x, M = c(2,300), interval = c(a,b), controls  =
        mable.ctrl(sig.level = 1e-8, maxit = 2000L, eps = 1.0e-7))
 op<-par(mfrow = c(1,2),lwd = 2)
 layout(rbind(c(1, 2), c(3, 3)))
 plot(res, which = "likelihood")
 plot(res, which = "change-point")
 hist(x, breaks = seq(0,7.5,len = 20), xlim = c(0,7), ylim = c(0,.7),
     prob  = TRUE,xlab = "t", ylab = "f(t)", col  = "light grey",
     main = "Histogram and Density of
               Duration of Eruptions of Old Faithful")
 lines(density(x, bw = "nrd0", adjust = 1), lty = 4, col = 4, lwd = 2)
 plot(res, which = "density", add = TRUE)
 lines(v, y1, lty = 2, col = 2, lwd = 2)
 legend("topright", lty = c(1,2,4), col = c(1,2,4), lwd = 2, bty = "n",
      c(expression(paste("MABLE: ",hat(f)[B](x))),
         expression(paste("Mixture: ",hat(f)[P](t))),
         expression(paste("KDE: ",hat(f)[K](t)))))
 par(op)

mable documentation built on Oct. 1, 2024, 9:06 a.m.