[R] Density

[William Dunlap]

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-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of li li
> Sent: Thursday, August 09, 2012 6:55 AM
> To: dcarlson at 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 <dcarlson at 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-bounces at r-project.org [mailto:r-help-bounces at 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 at 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 at 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.