I am trying to estimate Probability Density Function (PDF) of a big amount of data ($1e^6$ , $1e^7$, and higher) coming from Mote Carlo (MC) simulation.

My objective is to estimate the PDF (e.g. with Kernel Density estimation), and then extract the maximum point of this function that shall correspond to the mode of the population.

I am getting quite stuck in calculating this PDF (both with MATLAB using ksdensity function, and with R through fMultivar package and desnity2d function), because of the size of the population. Both MATLAB and R return error in allocating vectors of many GBs, and I suppose this is not the right way. I am saving the output population in .csv file.

To put you more in context.

The output of each MC iteration, namely $\bar{y}$, is the solution of an Optimal Control Problem where a 2D(bivariate) normal distrbution stochastic variable is inserted, i.e. $\bar{x}$. Therefore, the PDF of $\bar{x}$ is $\mathcal{N}\sim(\mathbf{\mu,\Sigma})$ (where $\mu$ and $\Sigma$ are the vector of expected value and diagonal covariance matrix of $\bar{x}$, respectively) and $\bar{y}$ is also a bidimensional variable with unknown distribution in principle . So my doubt is:

Once I run a MC simulation using numerous samples of $\bar{x}$ extracted from bivariate normal PDF, how can I calculate the PDF of $\bar{y}$ ?

I hope my question was clear, many thank you in advance for you useful support.


So you have a large iid sample $(y_{1,i},y_{2,i})$ from an unknown distribution and would like to estimate the pdf of the underlying distribution. Methodologically simple and computationally robust is a 2D-Histogram (see Matlab's histcounts2). You can run this in batches loading from a csv of any size. Since you have loads of data, you might want to run a recursive procedure, where you refine your histogram grid to zoom in on areas of largest density.

For example you start with a 20-by-20 grid and find one of the 400 resulting cells contains the largest amount of your points. So you know you your mode will probably be in this cell. Now you run the histogram procedure again, but only for the points in the cell. Say with a ten-by-ten grid. And so on. You should not make the grid too fine, since your density estimate (number of points in the cell divided by its area) will become more volatile as the cells contain lesser and lesser points.

Except for this, you will face similar challenges as for smooth density estimates. Your estimate might depend a lot on the way you set up your grid. For example ff your pairs are highly correlated you might also want to try a grid aligned with the principal components instead of x- and y-axis.

So you need to experiment to ensure robustness. It also helps a lot to develop some a-priori idea on what your distribution should look like. (e.g. Is it unimodal?) If you have a clear idea, you should also try a parametric fit. How to do this fit depends on your data of course. If it makes sense, start with a normal density and then try to distort it, if necessary, to improve the fit. From such a fit you can read off the mode and other properties easily and it serves to validate the non-parametric histogram approach.

  • $\begingroup$ "Since you have loads of data, you might want to run a recursive procedure, where you refine your histogram grid to zoom in on areas of largest density. "How can I do that ? What do you mean for ' iteratively refine your histogram, could you be more precise please ? Regarding the kind of output data I may suppose they underly on bivariate normal distribution, but it is a strong assumption. So you suggest to 'fit' the distribution with parametric estimate, validating the histcount approach and exctract my mode ? $\endgroup$ – thisisstan_. Jun 25 '19 at 19:32
  • $\begingroup$ See my edits, Whether the histcount or the fitting works better depends of course on the data and your ingenuity to come up with parametric families to fit. $\endgroup$ – g g Jun 25 '19 at 20:46

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.