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) 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
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
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.