kde {IsoplotR}

Create (a) kernel density estimate(s)

Description

Creates one or more kernel density estimates using a combination of the Botev (2010) bandwidth selector and the Abramson (1982) adaptive kernel bandwidth modifier.

Usage

kde(x, ...)

## Default S3 method:
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'other'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'UPb'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  type = 4,
  cutoff.76 = 1100,
  cutoff.disc = discfilter(),
  common.Pb = 0,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'detritals'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = FALSE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  ncol = NA,
  samebandwidth = TRUE,
  normalise = TRUE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'PbPb'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  common.Pb = 2,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'ArAr'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = FALSE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'KCa'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = FALSE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'ThPb'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = FALSE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'ThU'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [ka]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  Th0i = 0,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'ReOs'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = TRUE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'SmNd'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = TRUE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'RbSr'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = TRUE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'LuHf'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  i2i = TRUE,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'UThHe'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

## S3 method for class 'fissiontracks'
kde(
  x,
  from = NA,
  to = NA,
  bw = NA,
  adaptive = TRUE,
  log = FALSE,
  n = 512,
  plot = TRUE,
  rug = TRUE,
  xlab = "age [Ma]",
  ylab = "",
  kde.col = rgb(1, 0, 1, 0.6),
  hist.col = rgb(0, 1, 0, 0.2),
  show.hist = TRUE,
  bty = "n",
  binwidth = NA,
  hide = NULL,
  nmodes = 0,
  sigdig = 2,
  ...
)

Arguments

x

a vector of numbers OR an object of class UPb, PbPb, ThPb, ArAr, KCa, ReOs, SmNd, RbSr, UThHe, fissiontracks, ThU or detrital

...

optional arguments to be passed on to R's density function.

from

minimum age of the time axis. If NULL, this is set automatically

to

maximum age of the time axis. If NULL, this is set automatically

bw

the bandwidth of the KDE. If NULL, bw will be calculated automatically using the algorithm by Botev et al. (2010).

adaptive

logical flag controlling if the adaptive KDE modifier of Abramson (1982) is used

log

transform the ages to a log scale if TRUE

n

horizontal resolution (i.e., the number of segments) of the density estimate.

plot

show the KDE as a plot

rug

add a rug plot?

xlab

the x-axis label

ylab

the y-axis label

kde.col

the fill colour of the KDE specified as a four element vector of r, g, b, alpha values

hist.col

the fill colour of the histogram specified as a four element vector of r, g, b, alpha values

show.hist

logical flag indicating whether a histogram should be added to the KDE

bty

change to "o", "l", "7", "c", "u", or "]" if you want to draw a box around the plot

binwidth

scalar width of the histogram bins, in Myr if log = FALSE, or as a fractional value if log = TRUE. Sturges' Rule (\log_2[n]+1, where n is the number of data points) is used if binwidth = NA

hide

vector with indices of aliquots that should be removed from the plot.

nmodes

label the nmodes most prominent modes of the distribution. Change to 'all' to label all the modes.

sigdig

the number of significant digits to which the modes should be labelled. Only used if nmodes is a positive integer or 'all'.

type

scalar indicating whether to plot the ^{207}Pb/^{235}U age (type=1), the ^{206}Pb/^{238}U age (type=2), the ^{207}Pb/^{206}Pb age (type=3), the ^{207}Pb/^{206}Pb-^{206}Pb/^{238}U age (type=4), the concordia_age (type=5), or the ^{208}Pb/^{232}Th age (type=6).

cutoff.76

the age (in Ma) below which the ^{206}Pb/^{238}U and above which the ^{207}Pb/^{206}Pb age is used. This parameter is only used if type=4.

cutoff.disc

discordance cutoff filter. This is an object of class discfilter.

common.Pb

common lead correction:

0: none

1: use the Pb-composition stored in

settings('iratio','Pb207Pb206') (if x has class UPb and x$format<4);

settings('iratio','Pb206Pb204') and settings('iratio','Pb207Pb204') (if x has class PbPb or x has class UPb and 3<x$format<7); or

settings('iratio','Pb206Pb208') and settings('iratio','Pb207Pb208') (if x has class UPb and x$format=7,8).

2: use the isochron intercept as the initial Pb-composition

3: use the Stacey-Kramers two-stage model to infer the initial Pb-composition (only valid if x has class UPb).

ncol

scalar value indicating the number of columns over which the KDEs should be divided.

samebandwidth

logical flag indicating whether the same bandwidth should be used for all samples. If samebandwidth = TRUE and bw = NULL, then the function will use the median bandwidth of all the samples.

normalise

logical flag indicating whether or not the KDEs should all integrate to the same value.

i2i

‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘common’) ^{40}Ar/^{36}Ar, ^{40}Ca/^{44}Ca, ^{207}Pb/^{204}Pb, ^{87}Sr/^{86}Sr, ^{143}Nd/^{144}Nd, ^{187}Os/^{188}Os, ^{230}Th/^{232}Th, ^{176}Hf/^{177}Hf or ^{204}Pb/^{208}Pb ratio from an isochron fit. Setting i2i to FALSE uses the default values stored in settings('iratio',...).

Th0i

initial ^{230}Th correction.

0: no correction

1: project the data along an isochron fit

2: if x$format is 1 or 2, correct the data using the measured present day ^{230}Th/^{238}U, ^{232}Th/^{238}U and ^{234}U/^{238}U activity ratios in the detritus. If x$format is 3 or 4, correct the data using the measured ^{238}U/^{232}Th activity ratio of the whole rock, as stored in x by the read.data() function.

3: correct the data using an assumed initial ^{230}Th/^{232}Th-ratio for the detritus (only relevant if x$format is 1 or 2).

Details

Given a set of n age estimates \{t_1, t_2, ..., t_n\}, histograms and KDEs are probability density estimators that display age distributions by smoothing. Histograms do this by grouping the data into a number of regularly spaced bins. Alternatively, kernel density estimates (KDEs; Vermeesch, 2012) smooth data by applying a (Gaussian) kernel:

KDE(t) = \sum_{i=1}^{n}N(t|\mu=t_i,\sigma=h[t])/n

where N(t|\mu,\sigma) is the probability of observing a value t under a Normal distribution with mean \mu and standard deviation \sigma. h[t] is the smoothing parameter or ‘bandwidth’ of the kernel density estimate, which may or may not depend on the age t. If h[t] depends on t, then KDE(t) is known as an ‘adaptive’ KDE. The default bandwidth used by IsoplotR is calculated using the algorithm of Botev et al. (2010) and modulated by the adaptive smoothing approach of Abramson (1982). The rationale behind adaptive kernel density estimation is to use a narrower bandwidth near the peaks of the sampling distribution (where the ordered dates are closely spaced in time), and a wider bandwidth in the distribution's sparsely sampled troughs. Thus, the resolution of the density estimate is optimised according to data availability.

Value

If x has class UPb, PbPb, ArAr, KCa, ReOs, SmNd, RbSr, UThHe, fissiontracks or ThU, returns an object of class KDE, i.e. a list containing the following items:

x

horizontal plot coordinates

y

vertical plot coordinates

bw

the base bandwidth of the density estimate

ages

the data values from the input to the kde function

log

copied from the input

modes

a two-column matrix with the x and y values of the nmodes most prominent modes. Only returned if nmodes is a positive integer or 'all'.

or, if x has class =detritals, an object of class KDEs, i.e. a list containing the following items:

kdes

a named list with objects of class KDE

from

the beginning of the common time scale

to

the end of the common time scale

themax

the maximum probability density of all the KDEs

xlabel

the x-axis label to be used by plot.KDEs(...)

References

Abramson, I.S., 1982. On bandwidth variation in kernel estimates-a square root law. The annals of Statistics, pp.1217-1223.

Botev, Z. I., J. F. Grotowski, and D. P. Kroese. "Kernel density estimation via diffusion." The Annals of Statistics 38.5 (2010): 2916-2957.

Vermeesch, P., 2012. On the visualisation of detrital age distributions. Chemical Geology, 312, pp.190-194.

See Also

radialplot, cad

Examples

kde(examples$UPb)

dev.new()
kde(examples$FT1,log=TRUE)

dev.new()
kde(examples$DZ,from=1,to=3000,kernel="epanechnikov")