| mable.decon | R Documentation |
Maximum approximate Bernstein/Beta likelihood estimation in additive density deconvolution model with a known error density.
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
)
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 |
interval |
a finite vector |
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 |
controls |
Object of class |
progress |
if |
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).
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\}
Zhong Guan <zguan@iu.edu>
Guan, Z., (2019) Fast Nonparametric Maximum Likelihood Density Deconvolution Using Bernstein Polynomials, Statistica Sinica, doi:10.5705/ss.202018.0173
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.