FindDistributionParamters vs. NonlinearModelFit for NormalDistribution

Question

I have the following data:

dat["Sag"] ={1.277, 1.913, 1.48, 1.472, 1.787, 1.736, 1.896, 1.715, 1.131}

I want to fit a normal Distribution to it and get the according mean prediction bands. Letter one I can only get by using NonlinearModelFit. Both things are not a problem. Only when comparing the fit results between FindDistributionParameters and NonlinearModelFit I get way different results. Here the code:

\[ScriptCapitalH] = 
 DistributionFitTest[dat["Sag"], NormalDistribution[\[Mu], \[Sigma]], 
  "HypothesisTestData"];
\[ScriptCapitalH]["FittedDistributionParameters"]
(*{\[Mu]->1.60078,\[Sigma]->0.259774}*)
pdf[rad_] = PDF[\[ScriptCapitalH]["FittedDistribution"], rad];(*getting the PDF*)
(*creating data for fitting*)
hlist = HistogramList[dat["Sag"], 7, PDF];
fDat = Table[{(hlist[[1, n + 1]] + hlist[[1, n]])/2, 
    hlist[[2, n]]}, {n, Length[hlist[[2, All]]]}];
(*fitting the data*)
fitFkt = NonlinearModelFit[fDat, 
  PDF[NormalDistribution[\[Mu], \[Sigma]], sag], {\[Mu], \[Sigma]}, 
  sag];
fitFkt["BestFitParameters"]
(*{\[Mu] -> 1.71268, \[Sigma] -> 0.320319}->Wondering why that is different*)
(*Mean prediction bands*)
cl[0.95] = fitFkt["MeanPredictionBands", ConfidenceLevel -> 0.95][[2]];
cl[0.9] = fitFkt["MeanPredictionBands", ConfidenceLevel -> 0.9][[2]];
cl[0.5] = fitFkt["MeanPredictionBands", ConfidenceLevel -> 0.5][[2]];
(*Plotting fits*)
pdfPlot = 
 Plot[{pdf[sag], fitFkt[sag], cl[0.5], cl[0.9], cl[0.95]}, {sag, 0, 
   3}, Filling -> {2 -> {3}, 3 -> {4}, 4 -> {5}}, 
  PlotLegends -> "Expressions", 
  PlotLabel -> 
   "Comparison of fitting result by FindDistributionParameters
(pdf(sag) ) and NonlinearModelFit (fktFit(sag))\nincluding mean
prediction bands for confidence level (0.5, 0.9, 0.95) "]
(*Median und Standardabweichung*)
med["Sag"] = Median[dat["Sag"]]
std["Sag"] = StandardDeviation[dat["Sag"]]

The blue curve and the orange curve are the fitting results for FindDistributionParameters and NonlinearModelFit, respectively. As one can see the do not match at all. Interestingly, the NonlinearModelFit µ value fits better the Median. In opposition, MiniTab gets also the fit parameter values of FindDistributionParameters.

Do I something wrong with NonlinearModelFit?

My goal is to get a better idea about the uncertainty of the fit, as only 9 samples have been used. I tried that with the mean prediction bands. however, the fit doesn't match with FindDistributionParameters.

Is there an alternative way to get the uncertainty of the fit?

Using an estimation for the uncertainty of sigma:

HoldForm[\[CapitalDelta]\[Sigma]/\[Sigma] = 1/Sqrt[
   2 (N - 1)]] // TraditionalForm
1/Sqrt[2 (N - 1)] std["Sag"] + std["Sag"] /. N -> 9

In this way I would get Dsigma+sigma that includes the number of samples and would in that way account for the low number of samples in the fit result.

Letter one is what I am locking for.

Answer 1

In general I would say that fitting the PDF of a normal distribution to a histogram of your data using least square estimation (which is what NonlinearModelFit uses by default) is a pretty poor method for fitting a distribution to data. For one, least square estimation assumes that the fit residuals have a constant variance over the x-axis. This is assumption is most certainly not met when you try to fit a PDF to a histogram.

Furthermore, if you believe that there are no outliers in your data and that the data really comes from a normal distribution, then the correct way to estimate the distribution is just to calculate the mean and standard deviation. That's pretty much a defining characteristic of the normal distribution.

If you want to get a measure in the uncertainty of your fit parameters, I recommend a Bayesian analysis using conjugate priors. Or if you don't like Bayesian analysis you can just compute the confidence intervals of the mean and variance. For the mean, you can use the new function MeanAround for a quick shortcut:

In[120]:= MeanAround[dat["Sag"]]

1.60 +/- 0.09

edit You may find the following answer useful:

How do I calculate Standard Deviation confidence intervals in Mathematica?

Answer 2

@SjoerdSmit has given you the answer to obtain confidence intervals for the two parameters of the normal distribution. (I also agree that using NonlinearModelFit is a poor method to estimate the density except that I would go further and state that it is usually nonsensical to do so.)

You really shouldn't expect much with just 9 data points (no matter how expensive it was to obtain that data). Using a bootstrap approach one can obtain a "95% confidence band" for the probability density.

(* Data *)
x = {1.277, 1.913, 1.48, 1.472, 1.787, 1.736, 1.896, 1.715, 1.131};

(* Get mean and standard deviation of data *)
xbar = Mean[x];
xstdev = StandardDeviation[x];

(* Perform 1,000 bootstraps on a grid of 101 values *)
nBoot = 1000;
nGrid = 100;
(* Values to estimate density *)
z = xbar - 3 xstdev + 6 xstdev Range[0, nGrid]/nGrid;
density = ConstantArray[0, {nBoot, nGrid + 1}];
SeedRandom[12345];
Do[xBoot = RandomChoice[x, Length[x]];
 xBootMean = Mean[xBoot];
 xBootStdev = StandardDeviation[xBoot];
 density[[i]] = PDF[NormalDistribution[xBootMean, xBootStdev], #] & /@ z,
 {i, nBoot}]

(* Plot 95% confidence band *)
lower = Transpose[{z, Quantile[density[[All, #]], 0.025] & /@ Range[nGrid + 1]}];
upper = Transpose[{z, Quantile[density[[All, #]], 0.975] & /@ Range[nGrid + 1]}];
estimate = Transpose[{z, PDF[NormalDistribution[xbar, xstdev], #] & /@ z}];
ListPlot[{estimate, lower, upper}, Joined -> True, PlotStyle -> {Black, LightGray, LightGray}]

This gives you an indication that 9 sample points doesn't pin down the density very well.