0
$\begingroup$

I have written a small program that performs a Monte Carlo simulation to calculate the uncertainty of a quantity. I store all the results of this quantity in a list and from this list I generate a histogram.

h = 2*(Quantile[listR13A, 3/4] - Quantile[listR13A, 1/4])/nr^(1/3);
Show[Histogram[listR13A, {h}, "Probability"], PlotStyle -> Thick]

nr is the number of trials, listR13A is the above list where I store my results. h is the width of the bins. I wanted my own binning routine so that I know what is happening. Now would like to fit a Gaussian function to this histogram and obtain sigma and mu. How can I do this?

$\endgroup$
  • $\begingroup$ Use the raw data rather than the censored data placed into histogram bins to estimate the parameters of a normal distribution. $\endgroup$ – JimB Mar 11 at 16:00
  • $\begingroup$ okay, how do i get them? Sorry, I am a total beginner. BinLists, right? $\endgroup$ – Zorg Mar 11 at 16:01
  • $\begingroup$ By the raw data I mean the individual values of the quantity you're simulating. All you need is the sample mean and sample variance to estimate the population mean and population variance. That can be done by just collecting the sum of the simulated values, the sum of the squares of the individual values, and the total number of simulations. But if the number of simulations is not large (say less than 10,000,000), then you could keep all of the raw values and calculate the sample statistics at the end. $\endgroup$ – JimB Mar 11 at 16:13
  • $\begingroup$ okay, BinCounts and BinLists was not the right choice either. So what is it? $\endgroup$ – Zorg Mar 11 at 16:13
  • $\begingroup$ @JimB okay this is also a good approach, I will compare the results, if I mange to get want I like. but thank you! $\endgroup$ – Zorg Mar 11 at 16:18
4
$\begingroup$

Instead of fitting a function to the histogram (an estimate of the PDF), it is generally better to fit a distribution (not a function) to the raw, unbinned data.

Use FindDistributionParameters or EstimatedDistribution.


Example:

data = RandomVariate[NormalDistribution[], 200];

EstimatedDistribution[data, NormalDistribution[mean, sigma]]

You could do this, but DON'T. It is not reliable (it depends on the binning you use) and it does not give you a proper normalized PDF.

{bins, vals} = HistogramList[data, Automatic, "PDF"]

NonlinearModelFit[{MovingAverage[bins, 2], vals}\[Transpose], 
 c*Exp[-b*(x - m)^2], {b, c, m}, x]
| improve this answer | |
$\endgroup$
  • $\begingroup$ So what is the better option now? I did not get it. $\endgroup$ – Zorg Mar 11 at 17:09
  • $\begingroup$ @Zorg I showed two methods, and told you to use one, and not the other. Do you have any specific suggestions for improving the phrasing? $\endgroup$ – Szabolcs Mar 11 at 17:11
  • 1
    $\begingroup$ @Zorg If you're asking about FindDistributionParameters vs EstimatedDistribution: they do the same thing, they just return the result differently. $\endgroup$ – Szabolcs Mar 11 at 17:23
  • $\begingroup$ Okay, 'EstimatedDistribution' seems to be the best option for what I want. But I have encountered a different problem. I want to plot the histogram and the PDF in the same plot, and I like to have the histogram with the "Probability" option, but how can I normalize the PDF? $\endgroup$ – Zorg Mar 13 at 11:01

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.