densityFolded {mcmcOutput} | R Documentation |
Parameters are often constrained to be greater than zero (eg, standard deviation) or within the range (0, 1) (eg, probabilities), but the function density
often returns non-zero densities outside these ranges. Simple truncation does not work, as the area under the curve is < 1. The function densityFolded
attempts to identify these constraints and gives an appropriate density.
If x
is a matrix, detection of constraints and selection of bandwidth is applied to the pooled values, but a separate density curve is fitted to each column of the matrix.
densityFolded(x, bw = "nrd0", adjust = 1, from=NA, to=NA, ...)
x |
a numeric vector or matrix from which the estimate is to be computed; missing values not allowed. |
bw |
the smoothing bandwidth to be used; see |
adjust |
the bandwidth used is actually |
from , to |
the lower and upper ends of the grid at which the density is to be estimated; if NA, range will cover the values in x; ignored and replaced with 0 or 1 if a constraint is detected. |
... |
other arguments passed to |
Returns a list containing the following components:
the n coordinates of the points where the density is estimated.
a vector or matrix with the estimated density values.
the bandwidth used.
the sample size after elimination of missing values.
the call which produced the result.
the deparsed name of the x argument.
If y
is a vector, the output will have class density
.
Mike Meredith
require(graphics)
oldpar <- par(mfrow=2:1)
x1 <- rnorm(1e4) # no constraint on x1
plot(density(x1))
plot(densityFolded(x1)) # no difference
x2 <- abs(rnorm(1e4)) # x2 >= 0, with mode at 0
plot(density(x2)) # density > 0 when x2 < 0, mode around 0.2
abline(v=0, col='grey')
plot(densityFolded(x2)) # mode plotted correctly
abline(v=0, col='grey')
x3 <- rbeta(1e4, 1.5, 1.5) # 0 <= x3 <= 1
plot(density(x3)) # density > 0 when x2 < 0 and x2 > 1
abline(v=0:1, col='grey')
plot(densityFolded(x3))
abline(v=0:1, col='grey')
x4 <- rbeta(1e4, 1.5, 0.9) # 0 <= x4 <= 1, with mode at 1
plot(density(x4)) # mode appears to be around 0.95
abline(v=0:1, col='grey')
plot(densityFolded(x4)) # mode plotted correctly
abline(v=0:1, col='grey')
# Try with a matrix
x5 <- cbind(rbeta(1e4, 2,2), rbeta(1e4, 2,3), rbeta(1e4, 3,2))
plot(density(x5))
tmp <- densityFolded(x5)
with(tmp, matplot(x, y, type='l'))
par(oldpar)