# Estimating probability density function of big amount of data coming from MC simulations

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.