Estimating a probability distribution and sampling from it in Julia

Question

I am trying to use Julia to estimate a continuous univariate distribution using N observed data points (stored as an array of Float64 numbers), and then sample from this estimated distribution. I have no prior knowledge restricting attention to some family of distributions.

I was thinking of using the KernelDensity package to estimate the distribution, but I'm not sure how to sample from the resulting output.

Any help/tips would be much appreciated.

Answer 1

Without any restrictions on the estimated distribution, a natural candidate would be the empirical distribution function (see Wikipedia). For this distribution there are very nice theorems about convergence to actual distribution (see Dvoretzky–Kiefer–Wolfowitz inequality).

With this choice, sampling is especially simple. If dataset is a list of current samples, then dataset[rand(1:length(dataset),sample_size)] is a set of new samples from the empirical distribution. With the Distributions package, it could be more readable, like so:

using Distributions
new_sample = sample(dataset,sample_size)

Finally, Kernel density estimation is also good, but might need a parameter to be chosen (the kernel and its width). This shows a preference for a certain family of distributions. Sampling from a kernel distribution is surprisingly similar to sampling from the empirical distribution: 1. choose a sample from the empirical distributions; 2. perturb each sample using a sample from the kernal function.

For example, if the kernel function is a Normal distribution of width w, then the perturbed sample could be calculated as:

new_sample = dataset[rand(1:length(dataset),sample_size)]+w*randn(sample_size)