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

Improve this question
$\endgroup$
1
  • $\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$ Commented Jul 27, 2015 at 16:58

3 Answers 3

Reset to default
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.

Improve this answer
$\endgroup$
4
  • $\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$ Commented Jul 27, 2015 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$ Commented Jul 27, 2015 at 22:05
  • 3
    $\begingroup$ @infinitezero: By the way, that is known as the Rayleigh distribution. $\endgroup$ Commented Jul 27, 2015 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$ Commented Jul 28, 2015 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.

Improve this answer
$\endgroup$
2
  • $\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$ Commented Jul 28, 2015 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$ Commented Jul 28, 2015 at 1:56
6
$\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

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.