help mata mm_density()

Title

mm_density() -- Density estimation

Syntax

Initialize density estimation object

class mm_density scalar D
or
D = mm_density()

Setup

Provide data

D.data(X [, w, pw, sorted])

X is a real colvector containing data (missing values not allowed)

w is a real colvector containing weights (missing or negative values not allowed); default is w = 1 (no weights)

pw!=0 indicates that weights are sampling weights; default is pw = 0

sorted!=0 indicates that X is sorted and non-missing, that w is non-missing and
non-negative, and that X does not contain values that are outside of the support defined by
D.support(); default is sorted = 0; specify sorted!=0 to save some computer time
if you know that these conditions are true

Choose kernel

D.kernel(kernel [, adapt])

kernel is a string scalar specifying the kernel name; available
kernels are "epanechnikov",
"epan2", "biweight", "triweight",
"cosine", "gaussian", "parzen",
"rectangle", and "triangle"
(abbreviations allowed); default is "gaussian"; specifying ""
(empty string) selects the default

adapt is a real scalar specifying the number of stages of the adaptive
kernel density estimator; adapt may not be negative; default is adapt = 0
(non-adaptive kernel); specifying . (missing) selects the default

Set bandwidth or choose bandwidth selection method

D.bw(bw [, adjust, dpi, qui])

bw is a (strictly positive) real scalar specifying the value
of the bandwidth or a string scalar specifying the bandwidth selection
rule; available rules are "silverman" (optimal of Silverman),
"normalscale" (normal scale rule), "oversmoothed" (oversmoothed rule),
"sjpi" (Sheather-Jones solve-the-equation plug-in), "dpi"
(Sheather-Jones direct plug-in), and "isj" (bandwidth based on
diffusion estimator); default is "sjpi"; specifying "" (empty string)
or . (missing) selects the default

adjust is a real scalar specifying a multiplication factor
for the bandwidth; adjust must be strictly positive; default is
adjust = 1; specifying . (missing) selects the default

dpi is a real scalar specifying the number of stages of functional
estimation for the "dpi" rule; dpi may not be negative; default is
dpi = 2; specifying . (missing) selects the default

qui!=0 suppresses the warning message that is displayed if
SJPI or ISJ fails; default is qui = 0

Define support of data and choose boundary correction method

D.support(minmax [, method, rd])

minmax is a real vector specifying the lower and upper bounds
of the support of X; specify one value (lower boundary only) or
two values (lower and upper boundary); value . (missing) is interpreted as
(minus) infinity (i.e. unbounded); default is minmax = (.,.)
(unbounded support)

method is a string scalar specifying the boundary-correction
method to be used; available methods are "renormalization",
"reflection", and "linear correction" (abbreviations allowed); default
is "renormalization"; specifying ""
(empty string) selects the default

rd!=0 indicates that the data is to be interpreted as relative data (relative ranks); minmax defaults to
(0,1) in this case and automatic bandwidth selection is modified; values outside [0,1] are not allowed
in minmax if rd!=0; default is rd = 0

Set grid size of approximation estimator

D.n(n [, pad])

n is a real scalar specifying the grid size used by the fast approximation
estimator (and the bandwidth selectors); default
is n = 2^10 = 1024; specifying . (missing) selects the default

pad is a real scalar specifying the padding proportion of the approximation
grid; pad may not be negative; the default is pad = 0.1, which means that the grid will be padded
by 10 percent of the observed data range on each side of the data (unless the support of the
data has been restricted by D.support(), in which case the approximation grid
will span the defined support); specifying . (missing) selects the default

Retrieve settings

X = D.X() real colvector containing data

w = D.w() real colvector containing weights

nobs = D.nobs() real scalar containing number of observations (sum of weights)

pw = D.pw() real scalar containing pw flag

sorted = D.sorted() real scalar containing sorted flag

kernel = D.kernel() string scalar containing kernel name

adapt = D.adapt() real scalar containing number of stages of adaptive estimator

kh = D.kh() real scalar containing canonical bandwidth of kernel

bw = D.bw() string scalar containing name of bandwidth selection rule or real scalar containing user-provided bandwidth

adjust = D.adjust() real scalar containing bandwidth adjustment factor

dpi = D.dpi() real scalar containing number of DPI levels

minmax = D.support() real rowvector containing lower and upper bound of support

lb = D.lb() real scalar containing lower bound of support

ub = D.ub() real scalar containing upper bound of support

bc = D.bc() string scalar containing name of boundary-correction method

rd = D.rd() real scalar containing rd flag

n = D.n() real scalar containing grid size of approximation estimator

pad = D.pad() real scalar containing padding proportion of approximation grid

Technical remarks

o Changing any of the above settings will clear all results that may already exist in D.

o To preserve memory, D.data() does not store a copy of the
data in D; it only stores pointers to X and w. Modifying X
or w after D.data() has been applied will change the data used
by mm_density(), but it will not clear existing results in
D. This may lead to inconsistent results. It is best not to modify
X and w as long as D is active.

Estimation

Compute density estimate

d = D.d(at [, exact])
or
d = D.d(k, from, to [, exact])

at is a real vector specifying custom values at which the density is to be estimated

k is a real scalar specifying the size of the grid over which the density is to be estimated
(regular grid with k points)

from is a real scalar specifying the lowest value of the
grid; specify from >= . to determine the value automatically. Let xmin
be the minimum of the data, let tau be a padding value depending on
bandwidth and kernel, and let lb be the lower bound of the support as set
by D.support() (or minus infinity if unbounded). The lowest value of
the grid is then determined as follows:

if from < .: from

else if from = .z: xmin

else if from = .y: max(lb, xmin-tau)

else if bounded: lb

else: xmin-tau

to is a real scalar specifying the largest value of the
grid; specify to >= . to determine the value automatically (analogously to
from)

exact!=0 requests the exact density estimator to be used; default is
exact = 0 (approximation estimator)

Retrieve results after estimation

d = D.d() real colvector containing density estimate fitted last

at = D.at() real colvector containing evaluation points of density estimate

Self-threading functions

h = D.h() real scalar containing bandwidth

D = D.D() real colvector containing full-grid approximation estimator

AT = D.AT() real colvector containing evaluation grid of approximation estimator

W = D.W() real colvector containing grid counts of approximation estimator

L = D.L() real colvector containing local bandwidth factors of approximation estimator

l = D.l() real colvector containing local bandwidth factors of exact estimator

The self-threading functions can be used irrespective of whether D.d() has been applied or
not; results will be the same.

Other functions

Obtain levels or first derivatives of kernel function using current settings
(including boundary correction)

levels = D.K(X, at, h)
derivatives = D.Kd(X, at, h)

X is a real colvector containing the points at which the levels or derivatives will be evaluated

at is a real colvector containing kernel locations (midpoints)

h is a real colvector containing kernel bandwidths

X, at, and h are required to be r-conformable
(see [M-6] glossary ).

Description

mm_density() is a class-based system for univariate kernel density
estimation.

Examples

Density estimate with default settings:

. mata:
: x = rnormal(1000, 1, 0, 1)
: S = mm_density()
: S.data(x)
: S.d(13, -3, 3), S.at()
: end

Display full grid approximation estimator:

. mata:
: x = rnormal(50000,1, 0, 1) \ rnormal(50000,1, 4, 1.5)
: S = mm_density()
: S.data(x)
: mm_plot((S.D(), S.AT()), "line")
: end

Compute bandwidth (without estimating the density):

. mata:
: x = rnormal(1000,1, 0, 1) \ rnormal(1000,1, 4, 1.5)
: S = mm_density()
: S.data(x)
: S.bw("sj"); S.h() // Sheather-Jones solve-the-equation
: S.bw("dpi"); S.h() // Sheather-Jones direct plug-in
: S.bw("isj"); S.h() // diffusion estimator bandwidth
: S.bw("silver"); S.h() // optimal of Silverman
: S.bw("over"); S.h() // oversmooothed rule
: S.bw("normal"); S.h() // normal scale rule
: end

Source code

mm_density.mata

Methods and Formulas

For a brief overview of the methods used by mm_density() see
Jann (2007). A book-length
treatment can be found in Wand and Jones (1995). Not covered in these
references is the improved Sheather-Jones bandwidth selector based on diffusion
methods; see Botev et al. (2010).

References

Botev, Z.I., J.F. Grotowski, and D.P. Kroese (2010). Kernel density
estimation via diffusion. Annals of Statistics
38(5): 2916-2957. DOI: 10.1214/10-AOS799.

Jann, B. (2007). Univariate kernel density
estimation. DOI: 10.7892/boris.69421.

Wand, M.P., M.C. Jones (1995). Kernel Smoothing. London: Chapman and Hall.

Author

Ben Jann, University of Bern, ben.jann@soz.unibe.ch

Also see

Online: help for
mm_ddens(),
mm_ecdf(),
mm_histogram(),
kdensity,
histogram,
moremata