getting a density function from a sample of 500

Question

I'm pretty new to statistics so please excuse me if the answer is obvious.

The scenario is the following: I am using mcmc to sample from a posterior distribution of a parameter. I then need to estimate a quantity that depends on that parameter by integrating with respect to the posterior distribution of the parameter.

My question is, how can I get a density function from a sample of 500 points? I know I can get the cdf and the density is the derivative of it but I am unsure about how to do all of this in R.

Thanks!

EDIT:

One approach I was thinking about is to discretize the distribution and do a sum instead of an integral, but I'm wondering if there is a better way.

Answer 1

You need to get something like

$$ E[f(x)] = \int f(x) \, dx $$

MCMC was designed to approximate such integrals. By the law of large numbers

$$ E[f(x)]\approx \frac{1}{N} \sum_{i=1}^N f(x) $$

So the proper way is just to plug-in the MCMC samples from the posterior distribution to your function $f$ and calculate arithmetic mean (or median, quantiles, or whatever else you want).

As a sidenote, whatever you use MCMC for, if you want it to work, you need large samples, so the number of MCMC samples should be rather thousands, or tens of thousands, then hundreds. With 500 samples the estimate wouldn't be reliable for most of the cases. Moreover, usually the number of effective samples (e.g. measured by $n_{eff}$ statistic), i.e. the number of "independent" samples, would usually be much smaller then the number of MCMC draws, so your 500 samples can possibly have "value" of much smaller number of samples.