# Calculating (stepwise) cumulative probability from density in R

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. ## 1 Answer 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).

• 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. Dec 31 '16 at 3:39