mable.decon {mable}

Mable deconvolution with a known error density

Description

Maximum approximate Bernstein/Beta likelihood estimation in additive density deconvolution model with a known error density.

Usage

mable.decon(
  y,
  gn = NULL,
  ...,
  M,
  interval = c(0, 1),
  IC = c("none", "aic", "hqic", "all"),
  vanished = TRUE,
  controls = mable.ctrl(maxit.em = 1e+05, eps.em = 1e-05, maxit.nt = 100, eps.nt = 1e-10),
  progress = TRUE
)

Arguments

y

vector of observed data values

gn

error density function if known, default is NULL if unknown

...

additional arguments to be passed to gn

M

a vector (m0, m1) specifies the set of consective candidate model degrees, M = m0:m1. If gn is unknown then M a 2 x 2 matrix whose rows (m0,m1) and (k0,k1) specify lower and upper bounds for degrees m and k, respectively.

interval

a finite vector (a,b), the endpoints of supporting/truncation interval if gn is known. Otherwise, it is a 2 x 2 matrix whose rows (a,b) and (a1,b1) specify supporting/truncation intervals of X and \epsilon, respectively. See Details.

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

vanished

logical whether the unknown error density vanishes at both end-points of [a1,b1]

controls

Object of class mable.ctrl() specifying iteration limit and other control options. Default is mable.ctrl.

progress

if TRUE a text progressbar is displayed

Details

Consider the additive measurement error model Y = X + \epsilon, where X has an unknown distribution F on a known support [a,b], \epsilon has a known or unknown distribution G, and X and \epsilon are independent. We want to estimate density f = F' based on independent observations, y_i = x_i + \epsilon_i, i = 1, \ldots, n, of Y. We approximate f by a Bernstein polynomial model on [a,b]. If g=G' is unknown on a known support [a1,b1], then we approximate g by a Bernstein polynomial model on [a1,b1], a1<0<b1. We assume E(\epsilon)=0. AIC and BIC methods are used to select model degrees (m,k).

Value

A mable class object with components, if g is known,

• M the vector (m0, m1), where m1 is the last candidate degree when the search stoped

• m the selected optimal degree m

• p the estimate of p = (p_0, ..., p_m), the coefficients of Bernstein polynomial of degree m

• 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\}

• convergence An integer code. 0 indicates an optimal degree is successfully selected in M. 1 indicates that the search stoped at m1.

• 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

if g is unknown,

• M the 2 x 2 matrix with rows (m0, m1) and (k0,k1)

• nu_aic the selected optimal degrees (m,k) using AIC method

• p_aic the estimate of p = (p_0, ..., p_m), the coefficients of Bernstein polynomial model for f of degree m as in nu_aic

• q_aic the estimate of q = (q_0, ..., q_k), the coefficients of Bernstein polynomial model for g of degree k as in nu_aic

• nu_bic the selected optimal degrees (m,k) using BIC method

• p_bic the estimate of p = (p_0, ..., p_m), the coefficients of Bernstein polynomial model for f of degree m as in nu_bic

• q_bic the estimate of q = (q_0, ..., q_k), the coefficients of Bernstein polynomial model for g of degree k as in nu_bic

• lk matrix of log-likelihoods evaluated at m \in \{m_0, \ldots, m_1\} and k \in \{k_0, \ldots, k_1\}

• aic a matrix containing the Akaike information criterion(s) at m \in \{m_0, \ldots, m_1\} and k \in \{k_0, \ldots, k_1\}

• bic a matrix containing the Bayesian information criterion(s) at m \in \{m_0, \ldots, m_1\} and k \in \{k_0, \ldots, k_1\}

Author(s)

Zhong Guan <[e-mail]>

References

Guan, Z., (2019) Fast Nonparametric Maximum Likelihood Density Deconvolution Using Bernstein Polynomials, Statistica Sinica, doi:10.5705/ss.202018.0173

Examples

 # A simulated normal dataset
 set.seed(123)
 mu<-1; sig<-2; a<-mu-sig*5; b<-mu+sig*5;
 gn<-function(x) dnorm(x, 0, 1)
 n<-50;
 x<-rnorm(n, mu, sig); e<-rnorm(n); y<-x+e;
 res<-mable.decon(y, gn, interval = c(a, b), M = c(5, 50))
 op<-par(mfrow = c(2, 2),lwd = 2)
 plot(res, which="likelihood")
 plot(res, which="change-point", lgd.x="topright")
 plot(xx<-seq(a, b, length=100), yy<-dnorm(xx, mu, sig), type="l", xlab="x",
     ylab="Density", ylim=c(0, max(yy)*1.1))
 plot(res, which="density", types=c(2,3), colors=c(2,3))
 # kernel density based on pure data
 lines(density(x), lty=4, col=4)
 legend("topright", bty="n", lty=1:4, col=1:4,
 c(expression(f), expression(hat(f)[cp]), expression(hat(f)[bic]), expression(tilde(f)[K])))
 plot(xx, yy<-pnorm(xx, mu, sig), type="l", xlab="x", ylab="Distribution Function")
 plot(res, which="cumulative",  types=c(2,3), colors=c(2,3))
 legend("bottomright", bty="n", lty=1:3, col=1:3,
     c(expression(F), expression(hat(F)[cp]), expression(hat(F)[bic])))
 par(op)