I was wondering if I could get some help on this problem.

By simulating 10,000 times a stochastic system that yields a positive number at the end, I collect 10,000 random positive numbers. This is what I call a simulated stationary distribution.

I want to use this distribution to calculate an expectation with a value function associated with the same domain. I thought there could be two approaches: (1) find the best-fitting parametric function or (ii) make the numbers coarser and classify them into subintervals.

I’m trying to follow the second approach. Then the problem would be, for instance, putting 0.13333321 and 0.13312119 into the interval of [0.133,0.134) and counting how many rows are put in that interval to calculate the probability density function.

Would there be an existing package which does this? Or could anyone suggest how to achieve this goal efficiently? Any suggestion would be tremendously helpful.