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.