How can I approximate a pdf knowing the estimated CDF in R? I have an estimate of a CDF in R (nonparametric) and I need to compare this distribution to another one by Kullback-Leibler. In order to do so, I need to find the pdf of this random variable. What is the best way to estimate the pdf in this case? thank you.
 A: I have two different techniques for estimating the pdf.  The first technique just uses linear interpolation of the quantiles.  This seems to provide a decent estimate, but it is a little "choppy" depending on the underlying distribution. The second technique has an additional smoothing step which is accomplished via Monte Carlo simulation.  This seems to improve the estimate most of the time (not all the time though.  You are going to have to tweak with it a little more to see what you like, I attached R code below).  
Both techniques can be shown to be consistent estimators of the density, in that that will converge onto the true density as the sample size tends to infinity.  
Disclaimer:  I derived these myself last night and I have no references for you (not to say someone somewhere hasn't done something similar before).  I cannot guarantee anything about these techniques but my simulations suggest that they work well providing certain conditions are met (which I discuss more below).   
The Problem
To restate the problem; let $X \sim f(X)$. Your goal is to estimate $f(X)$. This would be easily accomplished with KDE if a sample $x_1,..,x_n \stackrel{iid}{\sim} f(X)$ was available.  
The dilemma here is that you do not observe a sample of data generated by $f(X)$, rather you observe another sample $y_1,...,y_n \stackrel{iid}{\sim} g(Y)$ and the quantiles $F(y_1),...F(y_n)$. For convenient notation also let $y_1,y_2,..,y_n$ be an ordered vector such that $y_1 < y_2 <...<y_n$.
Linear Interpolation
An enticing approach to estimating $f(x)$ would be to estimate the cdf, $F(x)=\int_{-\infty}^x f(t)dt$ and use the relationship  $\frac{\partial F(x)}{\partial x} = f(x)$.
One simple way to estimate the cdf is to use linear interpolation between known quantiles, i.e
$$
\hat F(x) = F(y_{k(x)}) + \bigg(\frac{x-y_{k(x)}}{y_{k(x)+1}-y_{k(x)}}\bigg)[F(y_{k(x)+1})-F(y_{k(x)})]
$$
where $k(x) \in \{1,2,...,n\}$ is defined such that $y_{k(x)} \le x < y_{k(x)+1}$.
Differentiating the above gives;
$$
\hat f(x) = \frac{F(y_{k(x)+1})-F(y_{k(x)})}{y_{k(x)+1}-y_{k(x)}}
$$
Smoothing and Monte Carlo
$\hat f(x)$ is a mixture of uniform distributions that steps up and down with changes in $k(x)$.  This can create a "choppy" binned curve, similar to a histogram in appearance.  Smoothing $\hat f(x)$ may increase it's accuracy.  However it would be challenging to use things like moving averages and local regressions because the density estimates should be positive and the law of total probability should be maintained.
The solution I have here is to use the probability integral transform (PIT) to simulate from $\hat F$ and let KDE do it's magic on the resulting sample.
Formally, I draw $u_1,...,u_m \sim Unif(0,1)$ then I caclulate
$$
\hat x_j =\hat F ^{-1}(u_j) =  \frac{u_j-\hat F(y_{k(\hat x_j)})}{\hat F(y_{k(\hat x_j)+1})-\hat F(y_{k(\hat x_j)})}(y_{k(\hat x_j)+1}-y_{k(\hat x_j)})+y_{k(\hat x_j)}
$$
Then I use the kernel density estimator (KDE) 
$$
\tilde {f}_h(x) = \frac{1}{mh} \sum_{j=1}^m K\Big(\frac{x-\hat x_j}{h}\Big)
$$
The KDE can be thought of a smoothing technique.
Caveats and Examples
Let $f(X) \equiv N(0,1)$ and draw $y_1,...,y_{260}$ from a shifted gamma distribution s.t. $z_i \sim \mathrm{Gamma}(1,.3)$ and $y_i =z_i-7$.  In the below plot the blue line is the pdf of $y_i$, the red line is the true pdf of $X$ (standard normal) and the black line is the estimate using the above techniques.

In this case the smoothed estimator seems more accurate.  Notice though, how each estimator is most accurate in the right tail of the distribution and that they loose accuracy on the left.   This is because the $y_i$'s become more dispersed as you move left (i.e.the $y_i$'s are distributed more densely on the right). 
In general, if the distribution of the $y_i$'s has very low density over an area where the density of $f(X)$ is high, the estimators will perform poorly.
Below is a more extreme case where $y_i \sim N(1,.6)$

Here you see that the smoothed estimator is completely off.  The linear interpolated estimator is actually accurate, but it doesn't estimate the entire distribution.
The problem above goes away if we just increase the variance of the $y_i$'s using $y_i \sim N(1,2)$; 

R Code
Play around with this
rm(list=ls())

n=260
x_pars=c(0,1);
y_pars=c(1,.3,7);

#observed data follows  gamma-shifted distribution
y=rgamma(n,y_pars[1],y_pars[2])-y_pars[3]
sy=sort(y)
#quantiles come from the standard normal distribution
p=pnorm(sy,x_pars[1],x_pars[2])

############################
# Linear Interpolation
############################
w=abs(diff(sy)) #interval width
dp=abs(diff(p)) #differenced quantiles
fx=dp/w         #estimated hight of pdf 
windows()
par(mfrow=c(2,1))
plot(sy[-n],fx,typ='l',main="Linear Interpolation",
      ylab="density",xlab="y",xlim=c(-3,3),ylim=c(0,.5)) #estimated pdf
points(sy,dnorm(sy,x_pars[1],x_pars[2]),type='l',
       col=2,lty=2) #actual pdf
points(sy,dgamma(sy+y_pars[3],y_pars[1],y_pars[2]),
       type='l',col=4,lty=4) #pdf of observations
legend("topright",legend=c("Actual PDF","Estimated PDF","Obs. PDF"),
       lty=c(2,1,4),col=c(2,1,4))

############################
# Monte Carlo Smoothing
############################
precision=100 #if this is to high the KDE may have difficulty + longer run-time
m=n*precision


u=sort(runif(m))

ind=which(u>p[1] & u<=p[n]) #throw these out for now (assume  0 density on either side of the observations)
m2=length(ind)
u=u[ind];

#simulate from linear interp. cdf using probability integral transform
z=numeric(m2)
st=1
for(i in 1:m2){
  for(j in 2:n){
    if(u[i]<=p[j]){
      st=j;
      interp=(u[i]-p[j-1])/dp[j-1]
      z[i]=sy[j-1]+w[j-1]*interp
      break
    }
  }
}


plot(density(z),type='l',main="Linear Interpolation + MC Smoothing",
     ylab="density",xlab="y",xlim=c(-3,3),ylim=c(0,.5)) #estimated pdf
points(sy,dnorm(sy,x_pars[1],x_pars[2]),type='l',
       col=2,lty=2) #actual pdf
points(sy,dgamma(sy+y_pars[3],y_pars[1],y_pars[2]),
       type='l',col=4,lty=4) #pdf of observations
legend("topright",legend=c("Actual PDF","Estimated PDF","Obs. PDF"),
       lty=c(2,1,4),col=c(2,1,4))

