5
$\begingroup$

I've been trying to find a solution to what I think is a common problem. I have a list of Data like

data = {{0.05, 15}, {0.15, 51}, {0.25, 64}, {0.35, 107}, {0.45, 113}, {0.55, 162}, {0.65, 163}, {0.75, 167}, {0.85, 182}, {0.95, 187}, {1.05, 165}, {1.15, 168}, {1.25, 151}, {1.35, 143}, {1.45, 121}, {1.55, 130}, {1.65, 109}, {1.75, 91}, {1.85, 91}, {1.95, 63}, {2.05, 48}, {2.15, 34}, {2.25, 29}, {2.35, 24}, {2.45, 14}, {2.55, 11}, {2.65, 6}, {2.75, 6}, {2.85, 9}, {2.95, 4}, {3.05, 4}, {3.15, 2}, {3.25, 0}, {3.35, 1}, {3.45, 5}, {3.55, 3}, {3.65, 2}, {3.75, 1}, {3.85, 2}, {3.95, 1}, {4.05, 0}, {4.15, 0}}

All I want is a Gaussian Distribution to be put over my ListPlot. I tried this, but the result was not even close to satisfying.

As I generate the data myself, I can change the format of data, if needed.

What I found and tried so far:

FindFit[data, Gaussian[x, ampl, x0, sigma], x, {ampl, x0, sigma}]

There I don't fulfill the criteria for the FindFit function as my number of coordinates is not equal to the number of variables.

model = PDF[NormalDistribution[μ, σ]]
sol = FindFit[data, model[x], {μ, σ}, x]

Please enlighten me

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey Jul 27 '15 at 16:58
9
$\begingroup$

Two problems: first, Gaussian is not a built-in function in Mathematica; you need to use the form involving PDF and NormalDistribution instead. Second, the arguments of FindFit need to be in the order data, model, parameters, variable; you have the last two switched. The correct version would be as follows:

model[x_] = ampl Evaluate[PDF[NormalDistribution[x0, sigma], x]];
fit = FindFit[data, model[x], {ampl, x0, sigma}, x]

(* {ampl -> 274.765, x0 -> 1.02404, sigma -> 0.614853} *)

Show[ListPlot[data], Plot[model[x] /. fit, {x, 0, 4.2}, PlotStyle -> Red]]

enter image description here

As you can see, the fit is not great at small $x$, largely because you have no data for negative $x$. You might want to think about whether this data is actually well-approximated by a Gaussian, or whether another distribution would be more appropriate.

$\endgroup$
  • $\begingroup$ To find out whether that data is well-approximated by a Gaussian or by another one is actually part of the problem. Thank you very much for your answer. I managed to test for a Cauchy-Distribution myself now. :) $\endgroup$ – infinitezero Jul 27 '15 at 17:17
  • $\begingroup$ Just in case someone is wondering. The function that describes this very good is b x e^(-a x^2) $\endgroup$ – infinitezero Jul 27 '15 at 22:05
  • 3
    $\begingroup$ @infinitezero: By the way, that is known as the Rayleigh distribution. $\endgroup$ – Rahul Jul 27 '15 at 23:30
  • $\begingroup$ I actually did a fit to $b x^n \exp (-a (x-x_0)^2 ) $ right after I posted my answer, and found $n \approx 1.003$ or something like that. $\endgroup$ – Michael Seifert Jul 28 '15 at 2:45
6
$\begingroup$

If you are trying to estimate the best fitting distribution you might also try something along these lines. I typically prefer to work with distributions through the powerful built-in framework. I've chosen a RayleighDistribution* since your data does not contain values below zero.

dist = EstimatedDistribution[WeightedData@@Transpose[data], RayleighDistribution[s]]

We can get a good estimate of the amplitude by taking the least squares estimate for a constant c that best scales the PDF to fit the original data.

x = PDF[dist, data[[All, 1]]];
y = data[[All, 2]];
c = 1/x.x x.y;

Now we can plot the rescaled PDF against the data.

Show[ListPlot[data], Plot[PDF[dist, x]*c, {x, 0, 4}], PlotRange -> All]

enter image description here

* I had initially chosen a WeibullDistribution but @Rahul expertly identified Rayleigh which is more parsimonious.

$\endgroup$
  • $\begingroup$ I didn't know about EstimatedDistribution, thanks for that! By the way, given that this is a probability distribution, the scaling factor ought to just be the total number of observations times the spacing between bins, i.e. c = 0.1 Total@data[[All,2]]. That gives $264.9$ compared to $266.675$ from your method; I don't have a good explanation for the difference. $\endgroup$ – Rahul Jul 28 '15 at 1:24
  • $\begingroup$ @Rahul they make different assumptions about what the data represent. In your case we assume that the data points are left endpoints on the bins of a histogram (a perfectly reasonable assumption in this case). In the case of least squares there is a presumption that the data contain noise. $\endgroup$ – Andy Ross Jul 28 '15 at 1:56
3
$\begingroup$

If the above data is binned from a random sample from a probability distribution, then using regression methods (FindFit, NonlinearModelFit, etc.) is at best bad form. For one thing, none of the regression assumptions hold (constant variance about the curve, etc.)

When one only has binned data and wants to estimate the parameters of a specified probability distribution, the maximum likelihood estimates for the parameters can be obtained by appropriately acknowledging that the data is "censored". In other words, the exact value of an observation is not known but rather only that the value is between two values.

Assuming a Rayleigh distribution suggested by @Rahul and implemented in the answer by @AndyRoss, label the desired cumulative distribution function $F$.

The probability of an observation being found in a bin with center $x_i$ and bin width $\Delta$ is given by

$$F\left(x_i+\frac{\Delta }{2}\right)-F\left(x_i-\frac{\Delta }{2}\right)$$

As we have $c_i$ observations in the $i$-th bin the likelihood function $$Likelihood = \prod _{i=1}^n \left(F\left(x_i+\frac{\Delta }{2}\right)-F\left(x_i-\frac{\Delta }{2}\right)\right)^{c_i}$$

with the log of the likelihood being

$$\log{Likelihood}=\sum _{i=1}^n c_i \log \left(F\left(x_i+\frac{\Delta }{2}\right)-F\left(x_i-\frac{\Delta }{2}\right)\right)$$

We can find the value of $\sigma$ that maximizes the log of the likelihood...

(* Log of the likelihood *)
binWidth = 0.1
logL = Total[#[[2]] Log[CDF[RayleighDistribution[σ], #[[1]] + binWidth/2] - 
        CDF[RayleighDistribution[σ], #[[1]] - binWidth/2]] & /@ data];

(* Starting value for σ using the sample mean *)
σ0 = 
 Sqrt[2/π] data[[All, 1]].data[[All, 2]]/Total[data[[All, 2]]]
(* 0.909224 *)

(* Find maximum likelihood estimate *)
sol = FindMaximum[{logL, σ > 0}, {{σ, σ0}}]
(* {-8345.23, {σ -> 0.90983}} *)

And thou shalt not produce an estimate without an estimate of precision:

(* Estimated standard error using the Delta Method *)
σse = ((-1/D[logL, {{σ}, 2}])[[1, 1]] /. sol[[2]])^0.5
(* 0.00884761 *)

Here are the fits using both a probability density scale and a count scale:

Show[ListPlot[
  Transpose[{data[[All, 1]], 
    data[[All, 2]]/(binWidth Total[data[[All, 2]]])}],
  Frame -> True, FrameLabel -> {"x", "Probability density"}], 
 Plot[PDF[RayleighDistribution[σ /. sol[[2]]], x], {x, 
   Min[data[[All, 1]]] - binWidth/2, 
   Max[data[[All, 1]]] + binWidth/2}]]

Fit using probability density scale

Show[ListPlot[Transpose[{data[[All, 1]], data[[All, 2]]}],
  Frame -> True, FrameLabel -> {"x", "Frequency count"}], 
 Plot[binWidth Total[data[[All, 2]]] PDF[
    RayleighDistribution[σ /. sol[[2]]], x], {x, 
   Min[data[[All, 1]]] - binWidth/2, 
   Max[data[[All, 1]]] + binWidth/2}]]

Fit using count scale

$\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.