Compute bivariate kernel density estimate using five parameter Gaussian kernels which can also use non equally spaced and adaptive bandwidths
Usage
KernSur(x, y, xgridsize=100, ygridsize=100, correlation, xbandwidth,
ybandwidth, range.x, range.y, na.rm=FALSE)
Arguments
x
vector of x values
y
vector of y values
xgridsize
integer for number of ordinates at which to calculate the smoothed estimate: default=100
ygridsize
integer for number of ordinates at which to calculate the smoothed estimate: default=100
correlation
x,y correlation, or vector of local correlations: default=cor(x,y)
xbandwidth
value of x window width, or vector of local window widths: default=dpik(x)
ybandwidth
value of y window width, or vector of local window widths: default=dpik(y)
range.x
total range of the estimate in the x dimension, or a vector giving the x ordinates: default=range +- 1.5 * mean bandwidth
range.y
total range of the estimate in the y dimension, or a vector giving the y ordinates: default=range +- 1.5 * mean bandwidth
na.rm
NA behaviour: TRUE drops cases with NA's, FALSE stops function with a warning if NA's are detected: default=FALSE
Value
returns two vectors and a matrix:
xords
vector of ordinates at which the density has been estimated in the x dimension
yords
vector of ordinates at which the density has been estimated in the y dimension
zden
matrix of density for f(x,y) with dimensions xgridsize, ygridsize
Acknowledgements
Written in collaboration with A.M.Pollard <mark.pollard@rlaha.ox.ac.uk> with the financial support of the Natural Environment Research Council (NERC) grant GR3/11395
Note
Slow code suitable for visualisation and display of correlated p.d.f, where highly generalised k.p.d.fs are needed - bkde2D is much faster when uncorrelated, uniformly grided, single bandwidth, k.p.d.fs are required.
This function doesn't use bins as such, it calculates the density at a set of points in each dimension. These points can be thought of as 'bin centres' but in reality they're not.
From version 1.00 onwards a number of improvements have been made: NA's are now handled semi-convincingly by dropping if required. A multi-element vector of bandwidths associated with each case can be sent for either dimension, so it is possible to accept the default, give a fixed bandwidth, or a bandwidth associated with each case. A multi-element vector of correlations can be sent, rather than a single correlation.
It should be noted that if a vector is sent for correlation, or either bandwidth, they must be of the same length as the data vectors. Furthermore, vectors which approximate the bin centres, can be sent rather than the extreme limits in the range; which means that the points at which the density is to be calculated need not be uniformly spaced.
Unlike KernSec this function does not yet support local bandwidths.
If the default bandwidth is to be used there must be at least five unique values for in the x and y vectors. If not the function will return an error. If you don't have five unique values in the vector then send a value, or vector for bandwidth
The number of ordinates defaults to the length of range.x if range.x is a vector of ordinates, otherwise it is xgridsize, or 100 if that isn't specified.
Finally, the various modes of sending parameters can be mixed, ie: the extremes of the range can be sent to define the range for x, but a multi-element vector could be sent to define the ordinates in the y dimension, or, a vector could be sent to describe the bandwidth for each case in the x direction, and a single-element vector defines all bandwidths in the y.
Version 1.1-0 has a bugfix in that it now outputs the magnetude of the density function at the specified bi-variate points, not an approximation to the volumes.
Lucy, D. Aykroyd, R.G. & Pollard, A.M.(2002) Non-parametric calibration for age estimation. Applied Statistics51(2): 183-196
See Also
KernSecperdensityhistbkdebkde2Ddpik
Examples
x <- c(2,4,6,8,10) # make up some x-y data
y <- x
# calculate and plot a surface with zero correlation based on above data
op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=0,
xbandwidth=1, ybandwidth=1, range.x=c(0,13), range.y=c(0,13))
image(op$xords, op$yords, op$zden, col=terrain.colors(100), axes=TRUE)
contour(op$xords, op$yords, op$zden, add=TRUE)
box()
# re-calculate and re-plot the above using a 0.8 correlation
op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=0.8,
xbandwidth=1, ybandwidth=1, range.x=c(0,13), range.y=c(0,13))
image(op$xords, op$yords, op$zden, col=terrain.colors(100), axes=TRUE)
contour(op$xords, op$yords, op$zden, add=TRUE)
box()
# calculate and plot a surface of the above data with an ascending
# correlation and bandwidths and a vector of equally spaced ordinates
bands <- c(1,1.1,1.2,1.3,1.0)
cors <- c(0,-0.2,-0.4,-0.6, -0.7)
rnge.x <- seq(from=0, to=13, length=100)
op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=cors,
xbandwidth=bands, ybandwidth=bands, range.x=rnge.x, range.y=c(0,13))
image(op$xords, op$yords, op$zden, col=terrain.colors(100), axes=TRUE)
contour(op$xords, op$yords, op$zden, add=TRUE)
box()
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(GenKern)
Loading required package: KernSmooth
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GenKern/KernSur.Rd_%03d_medium.png", width=480, height=480)
> ### Name: KernSur
> ### Title: Bivariate kernel density estimation
> ### Aliases: KernSur
> ### Keywords: distribution smooth
>
> ### ** Examples
>
> x <- c(2,4,6,8,10) # make up some x-y data
> y <- x
>
> # calculate and plot a surface with zero correlation based on above data
> op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=0,
+ xbandwidth=1, ybandwidth=1, range.x=c(0,13), range.y=c(0,13))
> image(op$xords, op$yords, op$zden, col=terrain.colors(100), axes=TRUE)
> contour(op$xords, op$yords, op$zden, add=TRUE)
> box()
>
> # re-calculate and re-plot the above using a 0.8 correlation
> op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=0.8,
+ xbandwidth=1, ybandwidth=1, range.x=c(0,13), range.y=c(0,13))
> image(op$xords, op$yords, op$zden, col=terrain.colors(100), axes=TRUE)
> contour(op$xords, op$yords, op$zden, add=TRUE)
> box()
>
> # calculate and plot a surface of the above data with an ascending
> # correlation and bandwidths and a vector of equally spaced ordinates
> bands <- c(1,1.1,1.2,1.3,1.0)
> cors <- c(0,-0.2,-0.4,-0.6, -0.7)
> rnge.x <- seq(from=0, to=13, length=100)
>
> op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=cors,
+ xbandwidth=bands, ybandwidth=bands, range.x=rnge.x, range.y=c(0,13))
> image(op$xords, op$yords, op$zden, col=terrain.colors(100), axes=TRUE)
> contour(op$xords, op$yords, op$zden, add=TRUE)
> box()
>
>
>
>
>
>
> dev.off()
null device
1
>