You can use approx(d, xout=) (or spline()) on the output of density() to get density estimates at the points in xout. E.g., > X <- c(rnorm(20, -1, 1), rgamma(50,1/2)) > d <- density(X) > plot(d) > points(approx(d, xout=-2:2))
Or you could use the functions in package:logspline to fit the density function and evaluate it where you wish > z <- logspline(X) > points(-2:2, dlogspline(-2:2, fit=z), col="red") (The vector '-2:2' about could be any set of numbers, including your original data points.) Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com > -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On > Behalf > Of li li > Sent: Thursday, August 09, 2012 6:55 AM > To: dcarl...@tamu.edu > Cc: r-help > Subject: Re: [R] Density > > Hi David, > Thanks a lot for the reply. > I might not have stated the problem clearly. Let me try again. > > Given a set of observations X, I want to find out the estimated density > values for the observations X? > > I believe that the values "x" returned from "density" function is not the > observations > that are fed into the function and the returned "y" values are estimated > density values for "x". > Below are in the R manual > > x > > the n coordinates of the points where the density is estimated. > y > > the estimated density values. These will be non-negative, but can be zero. > > We can also check this using the code below. > > X <- rnorm(100) > density(X)-> den0 > den0 > X[1:10] > (den0$x)[1:10] > (den0$y)[1:10] > round(dnorm((den0$x)[1:10]), 6) > round(dnorm(X[1:10]), 6) > > Thank you. > Hannah > > > 2012/8/8 David L Carlson <dcarl...@tamu.edu> > > > The numbers are there, they just aren't listed by the default print method > > for density. When you type the object name, den0, R runs > > print.density(den0) which provides summary statistics. You need to look at > > the manual page (?density) and pay close attention to the section labeled > > "Value" which provides information about what values are returned by the > > function. > > > > str(den0) > > den0$x > > den0$y > > plot(den0$x, den0$y, typ="l") > > > > ---------------------------------------------- > > David L Carlson > > Associate Professor of Anthropology > > Texas A&M University > > College Station, TX 77843-4352 > > > > > > > -----Original Message----- > > > From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- > > > project.org] On Behalf Of li li > > > Sent: Wednesday, August 08, 2012 9:03 PM > > > To: r-help > > > Subject: [R] Density > > > > > > Dear all, > > > Given a set of observations X, I want to evaluate the kernel density > > > estimator > > > at these observed values. If I do the following, I do not get the those > > > estimated values directly. > > > Can anyone familiar with this give an idea on how to find out the > > > estimated > > > density values for X? > > > > > > > X <- rnorm(100) > > > > > > > density(X)-> den0 > > > > > > > den0 > > > > > > > > > Call: > > > > > > density.default(x = X) > > > > > > > > > Data: X (100 obs.); Bandwidth 'bw' = 0.354 > > > > > > > > > x y > > > > > > Min. :-3.2254 Min. :0.0002658 > > > > > > 1st Qu.:-1.6988 1st Qu.:0.0359114 > > > > > > Median :-0.1721 Median :0.1438772 > > > > > > Mean :-0.1721 Mean :0.1635887 > > > > > > 3rd Qu.: 1.3545 3rd Qu.:0.2866889 > > > > > > Max. : 2.8812 Max. :0.3776935 > > > > > > > > > I did write the code for the kernel density > > > > > > estimator myself. So once I find a proper bandwidth, > > > > > > I can use the following function. However, it would be nicer > > > > > > if there is a more direct way. > > > > > > > > > > fhat <- function(x, X){ > > > > > > + h <- density(X, bw="SJ")$bw > > > > > > + n <- length(X) > > > > > > + 1/(n*h)*sum(dnorm((x-X)/h))} > > > > > > > > > > > > > > est <- numeric(length(X)) > > > > > > > for (i in 1:length(X)){est[i] <- fhat(x=X[i], X=X)} > > > > > > > > > > > > > > est > > > > > > > > > > > > > > > Thanks in advance. > > > > > > Hannah > > > > > > [[alternative HTML version deleted]] > > > > > > ______________________________________________ > > > R-help@r-project.org mailing list > > > https://stat.ethz.ch/mailman/listinfo/r-help > > > PLEASE do read the posting guide > > > http://www.R-project.org/posting-<http://www.r- > project.org/posting-> > > > guide.html > > > and provide commented, minimal, self-contained, reproducible code. > > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
