0
$\begingroup$

I want to Calculating (stepwise) cumulative probability from discontinuous density in R. The density was estimated based non-parametric method (density() function in R), and the command returned a class including x (coordinates of the points where the density is estimated), y (estimated density values based on x) and the selected bandwidth.

What I want is to construct a (stepwise) cumulative probability, which means, start from the lower bond of x (d.s$x[1]) and update the x with the bandwidth every time (Note the sample size of x is several times larger than the original sample itself), until it's larger than the upper bound. Finally, I can obtain a vector of estimated probability at the edge of band, and the probabilities of points within each interval are assumed to be identical. Furthermore, I can simulate a much more artificial sample and assign value to them according to the vector of probabilities and their own draws (between 0 and 1).

In my own case, I have an original sample of household's financial asset (130 observations). I use following code to conduct the estimation

library(foreign)
finasset <- read.dta("H:/CHFS/finasset2.dta")
asset = finasset$logasset
income = finasset$income



## MAIN
h.s  = sd(asset)*(4/(3*length(asset)))^.2 #Silverman’s rule of thumb,eq(10.25)
h.SJ = bw.SJ(asset)
h.t  = 3*sd(asset)*((1/(2*sqrt(pi)))/(35*length(asset)))^.2
d.s  = density(asset, bw=h.s,  kernel="gaussian")
d.SJ = density(asset, bw=h.SJ, kernel="gaussian")
d.t  = density(asset, bw=h.t,  kernel="gaussian")


## OUTPUT
h.s       # SILVERMAN'S BANDWIDTH
h.SJ      # SHEATHER-JONES BANDWIDTH (Slight variation in MASS package)
h.t       # TERRELL'S MAXIMAL SMOOTHING BANDWIDTH

(max(d.s$x)-min(d.s$x))/h.s shows the total supports of sample space can contain roughly 21 band. But there are 512 points for x, which means the coordinates are denser. I just want the distance between tow supports are near the bandwidth, or even larger. How to do it? By interpolation? or I can directly use the estimated x and y? And how to code it up in R (I'm relatively new to R)

One more question is, what should I do if I want to exclude those observations with value of zero, and assign the probability to them separately?


Thank you @Tim for your timely and helpful answer. Meanwhile, you write code in a compact way, this, however, makes it a little bit difficult.

For example, I'm messed up towards the part behind function(i): pnorm(grid, x[i], h)/length(x) is supposed to be normal distribution in each sub-interval, scaled by simple size(?) According to the syntax of pnorm, grid should be the real data(each asset observation in my case), x[i] should be mean, but how to decide it? And how to decide the length(or bandwidth) of each sub-interval, if it's not automatically selected by the non-parametric method?

Next, what about the next 3 arguments, numeric(length(grid))), 1, sum sum should mean add all up distributions in small intervals, what about the other two components? Can you elabrate it a little bit more.

$\endgroup$
1
$\begingroup$

Your question is not totally clear for me, but there is a simple answer to your main question. Kernel density estimate is in fact a mixture distribution

$$ f(x) = \sum_{i=1}^n p_i \, f_i(x; x_i, h) $$

where $p_i$ are mixing proportions, all equal to $1/n$ and $f_i$ distributions are your kernels, each parametrized by mean $x_i$ and bandwidth $h$. Cumulative distribution function of a mixture distribution is

$$ F(x) = \sum_{i=1}^n p_i \, F_i(x; x_i, h) $$

so with Gaussian kernel it is simply

$$ F(x) = \frac{1}{n} \sum_{i=1}^n \Phi(\tfrac{x-x_i}{h}) $$

where $\Phi$ is a standard normal cdf.

You do not need to use density function for it, instead you can simply use the direct algorithm:

# your data
x <- c(rnorm(50), rnorm(120, 1.5, 0.5), rnorm(50, -1, 0.1))
# points you want to evaluate your cdf on
grid <- seq(-6, 6, by = 0.01)
# bandwidth
h <- 0.1

p <- numeric(length(grid))
for (i in 1:length(x))
  p <- p + pnorm(grid, x[i], h)
p <- p/length(x)

This can be achieved in R also with more compact code:

apply(vapply(1:length(x), function(i) pnorm(grid, x[i], h)/length(x), numeric(length(grid))), 1, sum)

Below you can see the cdf of such distribution plotted against empirical cumulative distribution function of the simulated data (as above).

cdf of kernel density

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your helpful answer. However, since I'm new to R, your compact coding style is a little bit difficult to me. I add some doubts at the end of my original question, can you take a look. $\endgroup$ – zlqs1985 Dec 31 '16 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.