2
$\begingroup$

I feel like this should be a lot more straight-forward than it is... I have some 1D data that I can plot (data available here), and I want to fit a Gaussian to the peak such that I can extract a mean and standard deviation. There is a baseline background that you can see, and an obvious peak. I don't care if the peak really matches a different distribution, I just want to fit a Gaussian.

So... I've tried:

data = Import["MyData.txt","Table"];
data = data[[All,1]];
h=DistributionFitTest[data, Automatic, "HypothesisTestData"];
h["FittedDistribution"]
h["TestDataTable"]
Show[Histogram[data, Automatic, "PDF"], Plot[PDF[h["FittedDistribution"], x], {x, Min[g46Nqin], Max[g46Nqin]}, PlotStyle -> Thick, 

PlotRange -> Full]]

Which yields: NormalDistribution[46.4739, 4.56424]

enter image description here

This is obviously not a great fit. The peak is a little to the right, and the standard deviation is much tighter than shown.

So I tried making my own fit which includes the baseline:

nlm = NonlinearModelFit[data, {(\[CapitalAlpha]/(\[Sigma] Sqrt[2 Pi])) Exp[-(((x - \[Mu])/\[Sigma])^2)/2] + \[CapitalBeta], \[Mu] > 30, \[Sigma] > 0}, {{\[Sigma], 3.}, {\[Mu], 46}, \[CapitalAlpha], \[CapitalBeta]}, x]
Show[Show[Histogram[data, Automatic, "PDF"], Plot[nlm[x], {x, Min[data], Max[data]}]]]

And this just returns absolute crap. I've toyed with methods, starting parameter guesses, constraints, etc... Nothing helps.

Any options or assistance with this is highly valued! Thank you!

$\endgroup$
  • 1
    $\begingroup$ Not trying to be insulting but the question displays some very unconventional ideas about statistics. Mean[data] and StandardDeviation[data] ((Length[data] - 1)/Length[data])^0.5 gets you exactly the values that h["FittedDistribution"]. All fine so far. From the fit and histogram one sees that the distribution is far from normal. This is because your data isn't normally distributed and not because there's something wrong with the estimation procedures. $\endgroup$ – JimB Sep 29 '17 at 15:26
  • $\begingroup$ Fitting with NonlinearModelFit is confusing regression with fitting a distribution. Regression has nothing to do with fitting curve to a histogram. You want something like skd = SmoothKernelDistribution[data]; Show[Histogram[data, Automatic, "PDF"], Plot[PDF[skd, x], {x, Min[data], Max[data]}]]. $\endgroup$ – JimB Sep 29 '17 at 15:33
  • $\begingroup$ No offense taken, I am still learning a lot about stats (never had a formal class). Using an SKD, I was able to plot 2D data and fit to those points, so thanks! $\endgroup$ – Matt Stein Sep 29 '17 at 16:08
1
$\begingroup$

The problem is it is not NormalDistribution. So it's trying its best to find a possible fit and your data are messing it up. To find a best model you can do:

data=Import["https://drive.google.com/uc?export=download&id=0B-No0nABM5N8YnJEWDZLRDdwYXc"];
d1=data[[All,1]];
his1=Histogram[d1,Automatic,"Probability",PlotTheme->"Detailed"];
dis1=FindDistribution[d1,PerformanceGoal->"Quality"]
pl1=Plot[PDF[dis1,x],Prepend[MinMax[d1],x],PlotRange->All,PlotTheme->"Detailed"];
Show[his1,pl1]

Which gives not Gaussian:

Out[]= MixtureDistribution[{0.492401, 0.507599},
{UniformDistribution[{34.9459, 55.9893}], 
InverseGaussianDistribution[47.3888, 38032.8]}] 

enter image description here

I think it gave you a hint seen also by naked eye that UniformDistribution correction is needed. You can use it to force to fit a Gaussian:

his1=Histogram[d1,Automatic,"Probability",PlotTheme->"Detailed"];
dis1=FindDistribution[d1,TargetFunctions->{UniformDistribution,NormalDistribution},PerformanceGoal->"Quality"]
pl1=Plot[PDF[dis1,x],Prepend[MinMax[d1],x],PlotRange->All,PlotTheme->"Detailed"];
Show[his1,pl1]

Now you get:

Out[]= MixtureDistribution[{0.407219, 0.592781}, 
{UniformDistribution[{35.0021, 55.9962}], 
NormalDistribution[47.6887, 2.08879]}]

enter image description here

which is pretty good. BTW another function to keep on the radar is FindDistributionParameters - you have to really look through those docs article examples to understand best practices.

$\endgroup$
  • $\begingroup$ Absolutely brilliant! Thank you! Is there also a way to output the parameter tables on this such that I can get error on the parameters given? $\endgroup$ – Matt Stein Sep 29 '17 at 16:07
  • $\begingroup$ @MattStein The parameters are given in Out[]= - for example thsese are parameters NormalDistribution[47.6887, 2.08879] according to this definition PDF[NormalDistribution[μ, σ], x] (evaluate) or this article. $\endgroup$ – Vitaliy Kaurov Sep 29 '17 at 16:12
  • $\begingroup$ Yes I see those, but I was hoping it also had errors that were not immediately show. For example, the NonlinearModelFit function won't show errors immediately until you print out a property table... I am having trouble figuring out if that's possible with this method. $\endgroup$ – Matt Stein Sep 29 '17 at 16:22
  • $\begingroup$ @MattStein as I said already at the end of the post look at the details section and all examples in the DOCS to see related functionality and how to get it work. I am out of time now so might not respond quickly. $\endgroup$ – Vitaliy Kaurov Sep 29 '17 at 16:27
  • $\begingroup$ I've combed over that pretty extensively now, and it just doesn't appear to be there. I know that NonlinearModelFit gives that option, but there doesn't seem to be an equivalent for this method. Thank you for your time though and the excellent response. You answered the central issue. $\endgroup$ – Matt Stein Sep 29 '17 at 16:44
1
$\begingroup$

Update

I added a way to obtain approximate standard errors for the location of the peak in the probability density function at the end of answer.

If you want to estimate the mean and standard deviation from a random sample of just about any distribution, the sample mean and standard deviation is usually what you want:

Mean[data]
StandardDeviation[data]

But if the objective is to estimate where the maximum occurs of the (unknown) probability density function, then it depends on what you really know. If the distribution is Gaussian, then the best estimate where the peak occurs is the sample mean. If you "know" that the distribution is symmetrical then the sample mean or median will likely suffice.

But when you don't know either the distribution is Gaussian or symmetrical, then it's a bit more complicated (but not much more complicated if you have lots of data).

One can fit a nonparametric kernel density estimator and then find the location of the maximum.

skd = SmoothKernelDistribution[data];
sol = FindMaximum[PDF[skd, x], {{x, Mean[data]}}]
(* {0.13075329164814645, {x\[Rule]47.825414957389846}} *)

Show[Histogram[data, Automatic, "PDF"],
 Plot[PDF[skd, x], {x, Min[data], Max[data]}],
 ListPlot[{{x, 0}, {x, sol[[1]]}} /. sol[[2]],
  PlotStyle -> {Thick, Red}, Joined -> True]]

Histogram and density estimator and peak

Approximate confidence intervals for the location of the peak can be found by using a bootstrap procedure. ("All estimates should be accompanied by a measure of precision.") Here is an implementation that constructs standard errors for both the nonparametric kernel density approach and the mixture distribution approach suggested by @VitaliyKaurov:

nsim = 1000;
n = Length[data];
xmaxKernel = ConstantArray[0, nsim];
xmaxMixture = ConstantArray[0, nsim];
Do[
 (* Bootstrap sample *)
 bootstrap = RandomChoice[data, n];

 (* Nonparametric density estimator *)
 skd = SmoothKernelDistribution[bootstrap];
 Quiet[sol = FindMaximum[PDF[skd, x], {{x, Mean[data]}}]];
 xmaxKernel[[isim]] = x /. sol[[2]];

 (* Mixture approach *)
 dis1 = FindDistribution[bootstrap,
   TargetFunctions -> {UniformDistribution, NormalDistribution},
   PerformanceGoal -> "Quality"];
 Quiet[sol = FindMaximum[PDF[dis1, x], {{x, Mean[data]}}]];
 xmaxMixture[[isim]] = x /. sol[[2]],

 {isim, nsim}]

(* Standard deviations of the two estimation methods *)
sdKernel = StandardDeviation[xmaxKernel]
(* 0.08128476166850837 *)
sdMixture = StandardDeviation[xmaxMixture]
(* 0.24461412569292096 *)

The standard error for the mixture distribution approach is about 3 times the size of the nonparametric kernel approach.

$\endgroup$

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.