4

First off, I'm not entirely sure if this is the correct place to be posting this, as perhaps it should go in a more statistics-focussed forum. However, as I'm planning to implement this with R, I figured it would be best to post it here. Please apologise if I'm wrong.

So, what I'm trying to do is the following. I want to simulate data for a total of 250.000 observations, assigning a continuous (non-integer) value in line with a kernel density estimate derived from empirical data (discrete), with original values ranging from -5 to +5. Here's a plot of the distribution I want to use.

enter image description here

It's quite essential to me that I don't simulate the new data based on the discrete probabilities, but rather the continuous ones as it's really important that a value can be say 2.89 rather than 3 or 2. So new values would be assigned based on the probabilities depicted in the plot. The most frequent value in the simulated data would be somewhere around +2, whereas values around -4 and +5 would be rather rare.

I have done quite a bit of reading on simulating data in R and about how kernel density estimates work, but I'm really not moving forward at all. So my question basically entails two steps - how do I even simulate the data (1) and furthermore, how do I simulate the data using this particular probability distribution (2)?

Thanks in advance, I hope you guys can help me out with this.

nikUoM
  • 639
  • 1
  • 8
  • 18
  • If I understand your question, you could probably use Monte Carlo simulations. Since your distribution is known, you can use random sampling from that distribution to create your simulated data. Does that sound like what you're trying to do? – Lloyd Christmas Sep 26 '16 at 14:23
  • Hi @LloydChristmas, thanks for your reply. The thing is, the distribution is known, but I don't want to simulate the data using the actual, discrete distribution but rather the continuous one depicted in the kernel. If that's what you're referring to, this would be awesome. How can I do it? – nikUoM Sep 26 '16 at 14:40
  • 1
    There is some interesting discussion here: http://stats.stackexchange.com/questions/30303/how-to-simulate-data-that-satisfy-specific-constraints-such-as-having-specific-m – Lloyd Christmas Sep 26 '16 at 14:57
  • I was going to recommend a "rejection method", but that is covered in the link above. Also, @gung provided a solution that might work for you as well. – Lloyd Christmas Sep 26 '16 at 14:58
  • Hi. I've read through @gung's approach, and I don't quite understand most of it but think it might work. I'll read through the paper he/she recommended and see if I can produce something meaningful using this approach. Thank you very much! – nikUoM Sep 26 '16 at 15:34

1 Answers1

8

With your underlying discrete data, create a kernel density estimate on as fine a grid as you wish (i.e., as "close to continuous" as needed for your application (within the limits of machine precision and computing time, of course)). Then sample from that kernel density, using the density values to ensure that more probable values of your distribution are more likely to be sampled. For example:

Fake data, just to have something to work with in this example:

set.seed(4396)
dat = round(rnorm(1000,100,10))

Create kernel density estimate. Increase n if you want the density estimated on a finer grid of points:

dens = density(dat, n=2^14)

In this case, the density is estimated on a grid of 2^14 points, with distance mean(diff(dens$x))=0.0045 between each point.

Now, sample from the kernel density estimate: We sample the x-values of the density estimate, and set prob equal to the y-values (densities) of the density estimate, so that more probable x-values will be more likely to be sampled:

kern.samp = sample(dens$x, 250000, replace=TRUE, prob=dens$y)

Compare dens (the density estimate of our original data) (black line), with the density of kern.samp (red):

plot(dens, lwd=2)
lines(density(kern.samp), col="red",lwd=2)

enter image description here

With the method above, you can create a finer and finer grid for the density estimate, but you'll still be limited to density values at grid points used for the density estimate (i.e., the values of dens$x). However, if you really need to be able to get the density for any data value, you can create an approximation function. In this case, you would still create the density estimate--at whatever bandwidth and grid size necessary to capture the structure of the data--and then create a function that interpolates the density between the grid points. For example:

dens = density(dat, n=2^14)

dens.func = approxfun(dens)

x = c(72.4588, 86.94, 101.1058301)

dens.func(x)

[1] 0.001689885 0.017292405 0.040875436

You can use this to obtain the density distribution at any x value (rather than just at the grid points used by the density function), and then use the output of dens.func as the prob argument to sample.

eipi10
  • 91,525
  • 24
  • 209
  • 285