x <- rnorm(100) # simulate data
fit <- kde1d(x) # estimate density
-dkde1d(1000, fit) # evaluate density estimate
+dkde1d(0, fit) # evaluate density estimate
summary(fit) # information about the estimate
plot(fit) # plot the density estimate
curve(dnorm(x), add = TRUE, # add true density
@@ -117,7 +117,7 @@
Bounded data, log-linear
x <- rgamma(100, shape = 1) # simulate data
fit <- kde1d(x, xmin = 0, deg = 1) # estimate density
-dkde1d(1000, fit) # evaluate density estimate
+dkde1d(seq(0, 5, by = 1), fit) # evaluate density estimate
summary(fit) # information about the estimate
plot(fit) # plot the density estimate
curve(dgamma(x, shape = 1), # add true density
@@ -130,7 +130,7 @@
x <- rbinom(100, size = 5, prob = 0.5) # simulate data
x <- ordered(x, levels = 0:5) # declare as ordered
fit <- kde1d(x) # estimate density
-dkde1d(2, fit) # evaluate density estimate
+dkde1d(sort(unique(x)), fit) # evaluate density estimate
summary(fit) # information about the estimate
plot(fit) # plot the density estimate
points(ordered(0:5, 0:5), # add true density
diff --git a/docs/news/index.html b/docs/news/index.html
index 2d264db..14b9cd6 100644
--- a/docs/news/index.html
+++ b/docs/news/index.html
@@ -54,7 +54,7 @@
kde1d
- 0.1.2
+ 0.2.0
Changelog
Source: NEWS.md
+
+
NEW FEATURES
+
+- improved stability of density estimates near a boundary (#21).
+
+
BUG FIXES
+
+consistent behavior when dkde1d() and pkde1d() are called with non-ordered input although data are discrete (#19).
fixed bug in computation of kernel density estimates (#20).
adapt minimum bw allowed for discrete data to truncated Gaussian kernel (#20).
+
See a
Examples
- ## For reproducibility
-
set.seed(
0)
-
+
#> [1] 1e-12
summary(fit) # information about the estimate
#> kernel density estimate ('kde1d'), log-quadratic
+
dkde1d(
0,
fit)
# evaluate density estimate#> [1] 0.4123974
summary(fit) # information about the estimate
#> kernel density estimate ('kde1d'), log-quadratic
#> -----------------------------------------------------------------
-#> nobs = 100, bw = 0.37, loglik = -128.39, d.f. = 2.24
plot(fit) # plot the density estimate
curve(dnorm(x), add = TRUE, # add true density
+#> nobs = 500, bw = 0.29, loglik = -700.87, d.f. = 10.58
plot(fit) # plot the density estimate
curve(dnorm(x), add = TRUE, # add true density
col = "red")
#> [1] 482652.4
summary(fit) # information about the estimate
#> kernel density estimate ('kde1d'), log-linear with bounded support (xmin = 0)
+
dkde1d(
seq(
0,
5,
by = 1),
fit)
# evaluate density estimate#> [1] 1.144082921 0.318018791 0.108508811 0.036081984 0.009148931 0.001913149
summary(fit) # information about the estimate
#> kernel density estimate ('kde1d'), log-linear with bounded support (xmin = 0)
#> -----------------------------------------------------------------
-#> nobs = 100, bw = 0.46, loglik = -98.75, d.f. = 3.51
plot(fit) # plot the density estimate
curve(dgamma(x, shape = 1), # add true density
+#> nobs = 500, bw = 0.35, loglik = -471.12, d.f. = 39.2
plot(fit) # plot the density estimate
curve(dgamma(x, shape = 1), # add true density
add = TRUE, col = "red",
from = 1e-3)
#> [1] 0.332243
summary(fit) # information about the estimate
#> (jittered) kernel density estimate ('kde1d'), log-quadratic
+
dkde1d(
sort(
unique(
x)),
fit)
# evaluate density estimate#> [1] 0.04832389 0.15290500 0.28419148 0.30222901 0.18384615 0.02850447
summary(fit) # information about the estimate
#> (jittered) kernel density estimate ('kde1d'), log-quadratic
#> -----------------------------------------------------------------
-#> nobs = 100, bw = 0.5, loglik = -154.49, d.f. = 2.94
plot(fit) # plot the density estimate
points(ordered(0:5, 0:5), # add true density
+#> nobs = 500, bw = 0.36, loglik = -787.47, d.f. = 14.62
plot(fit) # plot the density estimate
points(ordered(0:5, 0:5), # add true density
dbinom(0:5, 5, 0.5), col = "red")
diff --git a/docs/reference/plot.kde1d-1.png b/docs/reference/plot.kde1d-1.png
index 966ed18..37bb3a2 100644
Binary files a/docs/reference/plot.kde1d-1.png and b/docs/reference/plot.kde1d-1.png differ
diff --git a/docs/reference/plot.kde1d-2.png b/docs/reference/plot.kde1d-2.png
index 5070e47..5804940 100644
Binary files a/docs/reference/plot.kde1d-2.png and b/docs/reference/plot.kde1d-2.png differ
diff --git a/docs/reference/plot.kde1d.html b/docs/reference/plot.kde1d.html
index f8b9d23..7bcf736 100644
--- a/docs/reference/plot.kde1d.html
+++ b/docs/reference/plot.kde1d.html
@@ -57,7 +57,7 @@
kde1d
- 0.1.2
+ 0.2.0
diff --git a/inst/include/lpdens.hpp b/inst/include/lpdens.hpp
index 95b856e..11cfdb0 100644
--- a/inst/include/lpdens.hpp
+++ b/inst/include/lpdens.hpp
@@ -138,7 +138,7 @@ inline Eigen::MatrixXd LPDens1d::fit_lp(const Eigen::VectorXd& x_ev,
Eigen::VectorXd kernels(x.size());
for (size_t k = 0; k < x_ev.size(); k++) {
// classical (local constant) kernel density estimate
- xx = x.array() - x_ev(k);
+ xx = (x.array() - x_ev(k)) / bw_;
kernels = kern_gauss(xx) / bw_;
f0 = kernels.mean();
res(k, 0) = f0;
@@ -149,13 +149,12 @@ inline Eigen::MatrixXd LPDens1d::fit_lp(const Eigen::VectorXd& x_ev,
f1 = xx.cwiseProduct(kernels).mean(); // first order derivative
b = f1 / f0;
- if (deg_ > 1) {
+ if (deg_ == 2) {
// more calculations for local quadratic
xx2 = xx.cwiseProduct(kernels) / (f0 * static_cast(n));
b *= std::pow(bw_, 2);
- s = 1.0 / (std::pow(bw_, 4) * xx.cwiseProduct(xx2).sum() -
- std::pow(b, 2));
- res(k, 0) *= std::sqrt(s) / bw_;
+ s = 1.0 / (std::pow(bw_, 4) * xx.transpose() * xx2 - std::pow(b, 2));
+ res(k, 0) *= bw_ * std::sqrt(s);
}
// final estimate
@@ -226,7 +225,7 @@ inline Eigen::VectorXd LPDens1d::boundary_transform(const Eigen::VectorXd& x,
if (!inverse) {
if (!std::isnan(xmin_) & !std::isnan(xmax_)) {
// two boundaries -> probit transform
- x_new = (x.array() - xmin_ + 5e-3) / (xmax_ - xmin_ + 1e-2);
+ x_new = (x.array() - xmin_ + 5e-5) / (xmax_ - xmin_ + 1e-4);
x_new = stats::qnorm(x_new);
} else if (!std::isnan(xmin_)) {
// left boundary -> log transform
@@ -240,8 +239,8 @@ inline Eigen::VectorXd LPDens1d::boundary_transform(const Eigen::VectorXd& x,
} else {
if (!std::isnan(xmin_) & !std::isnan(xmax_)) {
// two boundaries -> probit transform
- x_new = stats::pnorm(x).array() + xmin_ - 5e-3;
- x_new *= (xmax_ - xmin_ + 1e-2);
+ x_new = stats::pnorm(x).array() + xmin_ - 5e-5;
+ x_new *= (xmax_ - xmin_ + 1e-4);
} else if (!std::isnan(xmin_)) {
// left boundary -> log transform
x_new = x.array().exp() + xmin_ - 1e-3;
@@ -267,10 +266,10 @@ inline Eigen::VectorXd LPDens1d::boundary_correct(const Eigen::VectorXd& x,
Eigen::VectorXd corr_term(fhat.size());
if (!std::isnan(xmin_) & !std::isnan(xmax_)) {
// two boundaries -> probit transform
- corr_term = (x.array() - xmin_ + 5e-3) / (xmax_ - xmin_ + 1e-2);
+ corr_term = (x.array() - xmin_ + 5e-5) / (xmax_ - xmin_ + 1e-4);
corr_term = stats::dnorm(stats::qnorm(corr_term));
- corr_term /= (xmax_ - xmin_ + 1e-2);
- corr_term = 1.0 / corr_term.array().max(1e-4);
+ corr_term /= (xmax_ - xmin_ + 1e-4);
+ corr_term = 1.0 / corr_term.array().max(1e-6);
} else if (!std::isnan(xmin_)) {
// left boundary -> log transform
corr_term = 1.0 / (1e-3 + x.array() - xmin_);
@@ -330,6 +329,7 @@ inline Eigen::VectorXd LPDens1d::construct_grid_points(const Eigen::VectorXd& x)
//! @param grid_points the grid points.
inline Eigen::VectorXd LPDens1d::finalize_grid(Eigen::VectorXd& grid_points)
{
+ double range = grid_points.maxCoeff() - grid_points.minCoeff();
if (!std::isnan(xmin_))
grid_points(0) = xmin_;
if (!std::isnan(xmax_))
diff --git a/man/kde1d.Rd b/man/kde1d.Rd
index 9cee9bc..5261b45 100644
--- a/man/kde1d.Rd
+++ b/man/kde1d.Rd
@@ -40,22 +40,20 @@ log-transform is used if there is only one boundary (see, Geenens and Wang,
Discrete variables are handled via jittering (see, Nagler, 2018a, 2018b).
}
\examples{
-## For reproducibility
-set.seed(0)
## unbounded data
-x <- rnorm(100) # simulate data
+x <- rnorm(500) # simulate data
fit <- kde1d(x) # estimate density
-dkde1d(1000, fit) # evaluate density estimate
+dkde1d(0, fit) # evaluate density estimate
summary(fit) # information about the estimate
plot(fit) # plot the density estimate
curve(dnorm(x), add = TRUE, # add true density
col = "red")
## bounded data, log-linear
-x <- rgamma(100, shape = 1) # simulate data
+x <- rgamma(500, shape = 1) # simulate data
fit <- kde1d(x, xmin = 0, deg = 1) # estimate density
-dkde1d(1000, fit) # evaluate density estimate
+dkde1d(seq(0, 5, by = 1), fit) # evaluate density estimate
summary(fit) # information about the estimate
plot(fit) # plot the density estimate
curve(dgamma(x, shape = 1), # add true density
@@ -63,10 +61,10 @@ curve(dgamma(x, shape = 1), # add true density
from = 1e-3)
## discrete data
-x <- rbinom(100, size = 5, prob = 0.5) # simulate data
+x <- rbinom(500, size = 5, prob = 0.5) # simulate data
x <- ordered(x, levels = 0:5) # declare as ordered
fit <- kde1d(x) # estimate density
-dkde1d(2, fit) # evaluate density estimate
+dkde1d(sort(unique(x)), fit) # evaluate density estimate
summary(fit) # information about the estimate
plot(fit) # plot the density estimate
points(ordered(0:5, 0:5), # add true density
diff --git a/src/wrappers.cpp b/src/wrappers.cpp
index b003103..91c8a19 100644
--- a/src/wrappers.cpp
+++ b/src/wrappers.cpp
@@ -54,7 +54,17 @@ InterpolationGrid1d wrap_to_cpp(const Rcpp::List& R_object)
Eigen::VectorXd dkde1d_cpp(const Eigen::VectorXd& x,
const Rcpp::List& R_object)
{
- return wrap_to_cpp(R_object).interpolate(x);
+ Eigen::VectorXd fhat = wrap_to_cpp(R_object).interpolate(x);
+ double xmin = R_object["xmin"];
+ double xmax = R_object["xmax"];
+ if (!std::isnan(xmin)) {
+ fhat = (x.array() < xmin).select(Eigen::VectorXd::Zero(x.size()), fhat);
+ }
+ if (!std::isnan(xmax)) {
+ fhat = (x.array() > xmax).select(Eigen::VectorXd::Zero(x.size()), fhat);
+ }
+
+ return fhat;
}
//' computes the cdf of a kernel density estimate by numerical integration.
@@ -117,7 +127,7 @@ double select_bw_cpp(const Eigen::VectorXd& x,
bw *= mult;
if (discrete) {
- bw = std::max(bw, 0.5);
+ bw = std::max(bw, 0.5 / 5);
}
return(bw);
diff --git a/tests/testthat/test_kde1d.R b/tests/testthat/test_kde1d.R
index ae1a635..4451888 100644
--- a/tests/testthat/test_kde1d.R
+++ b/tests/testthat/test_kde1d.R
@@ -73,6 +73,15 @@ test_that("d/p/r/h functions work", {
expect_lte(max(pkde1d(sim, fit), 1), 1)
expect_that(all(qkde1d(u, fit) >= xmin), equals(TRUE))
expect_that(all(qkde1d(u, fit) <= xmax), equals(TRUE))
+ if (!is.nan(fit$xmin)) {
+ expect_equal(dkde1d(xmin - 1, fit), 0)
+ expect_equal(pkde1d(xmin - 1, fit), 0)
+ }
+
+ if (!is.nan(fit$xmax)) {
+ expect_equal(dkde1d(xmax + 1, fit), 0)
+ expect_equal(pkde1d(xmax + 1, fit), 1)
+ }
}
sims <- lapply(fits, function(x) rkde1d(n, x))
@@ -108,4 +117,22 @@ test_that("other generics work", {
}
lapply(fits, test_other_generics)
+})
+
+test_that("behavior for discrete data is consistent", {
+ n <- 1e3
+ x <- ordered(sample(5, n, TRUE), 1:5)
+ fit <- kde1d(x)
+ xx <- ordered(1:5, 1:5)
+ expect_equal(dkde1d(1:5, fit), dkde1d(xx, fit))
+ expect_equal(pkde1d(1:5, fit), pkde1d(xx, fit))
+ expect_true(all(is.na(dkde1d(c(0, 6), fit))))
+ expect_true(all(is.na(pkde1d(c(0, 6), fit))))
+ expect_true(all(rkde1d(n, fit) %in% x))
+})
+
+test_that("estimates for discrete data are reasonable", {
+ x <- ordered(sample(5, 1e5, TRUE), 1:5)
+ fit <- kde1d(x)
+ expect_true(all(abs(dkde1d(1:5, fit) - 0.2) < 0.1))
})
\ No newline at end of file
+dkde1d(0, fit) # evaluate density estimate (close to dnorm(0))#> [1] 0.4780549pkde1d(0, fit) # evaluate corresponding cdf (close to pnorm(0))#> [1] 0.5156368qkde1d(0.5, fit) # quantile function (close to qnorm(0))#> [1] -0.03249129hist(rkde1d(100, fit)) # simulate
