0
$\begingroup$

The data is given below:

{{-2.9, 1}, {-2.7, 0}, {-2.5, 0}, {-2.3, 2}, {-2.1, 2}, {-1.9,
3}, {-1.7, 5}, {-1.5, 7}, {-1.3, 3}, {-1.1, 11}, {-0.9, 7}, {-0.7, 
3}, {-0.5, 14}, {-0.3, 9}, {-0.1, 24}, {0.1, 17}, {0.3, 26}, {0.5, 
11}, {0.7, 14}, {0.9, 11}, {1.1, 9}, {1.3, 5}, {1.5, 2}, {1.7, 
5}, {1.9, 3}, {2.1, 3}, {2.3, 1}, {2.5, 1}, {2.7, 1}, {2.9, 0}}

I have used the formula of NormalDistribution with $\mu$ and $\sigma$ as the parameters of the fitting through NonlinearModelFit and the result is follows

gaussian dist

The data should satisfy NormalDistribution in a much more decent form as can be seen by eye estimation and moreover the data is taken from some nuclear experiments satisfying the distribution. I guess there might be some other ways to fit these kind of statistical data points.

Any guidance regarding any fitting issues will be very helpful to resolve my problem. Thanks for all valuable suggestions in advance !

$\endgroup$
  • 1
    $\begingroup$ NonlinearModelFit[] seems inappropriate here; have a look at EstimatedDistribution[] or FindDistributionParameters[]. $\endgroup$ – J. M. is away Aug 10 '17 at 18:14
  • $\begingroup$ I assume that the integer values are counts of observations? Why do you need to fit any particular distribution? Is it to compare parameter estimates from other datasets? $\endgroup$ – JimB Aug 10 '17 at 19:40
  • $\begingroup$ You cannot fit a probability density to your data since the area under your data is not 1. You want to fit to a gaussian pulse as shown by David Stork. $\endgroup$ – Bob Hanlon Aug 10 '17 at 19:45
  • 1
    $\begingroup$ @BobHanlon, if the OP's data are the midpoints and frequency counts of a histogram, then a regression approach is completely inappropriate. (I'm not as diplomatic as J.M.) $\endgroup$ – JimB Aug 10 '17 at 19:51
3
$\begingroup$

If one wants to fit a curve that just happens to be of the form of a standard probability density function (with an additive and/or multiplicative constant included), then a regression approach makes sense. But there is no probabilistic interpretation that is induced by such a fit.

But when one has a random sample from that probability distribution, performing a regression makes no sense.

The data presented for the above question is almost certainly pairs of midpoints and the associated frequency of occurrence from a sample from some distribution. The question is about fitting a normal distribution from such data. (The goodness of the fit is another issue.)

If one had the raw data and was fitting a normal distribution, then calculating the sample mean and the sample standard deviation would be a reasonable way to estimate the underlying parameters. When only the binned data is available there are two typical ways to perform the estimation for a normal distribution.

(* Get estimates of the mean and standard deviation using the binned data *)
μ0 = data[[All, 1]].data[[All, 2]]/Total[data[[All, 2]]]
(* 0.045999999999999944 *)
σ0 = (((data[[All, 1]]^2).data[[All, 2]])/Total[data[[All, 2]]] - μ0^2)^0.5
(* 1.0017404853553638 *)

The maximum likelihood estimates of $\mu$ and $\sigma$ based on the binned data are found by maximizing the likelihood (or equivalently the log of the likelihood):

$$\text{Likelihood}=\prod_{i=1}^n \left(\Phi\left(\frac{x_i+\frac{w}{2}-\mu}{2}\right)-\Phi\left(\frac{x_i-\frac{w}{2}-\mu}{2}\right)\right)^{f_i}$$

where $x_i$ is the $i$-th bin midpoint, $f_i$ is the frequency count of the $i$-th bin, $n$ is the number of bins, $w$ is the bin width, and $\Phi(z)$ the standard normal cumulative distribution function.

(* Get maximum likelihood estimates *)
Φ[z_NumericQ] := CDF[NormalDistribution[0, 1], z]
w = 0.2; (* Width of histogram bin *)

(* Log of the likelihood *)
logL = Sum[data[[i, 2]] Log[Φ[(data[[i, 1]] + w/2 - μ)/σ] - Φ[(data[[i, 1]] - w/2 - μ)/σ]],
 {i, Length[data]}];

(* Find values of the parameters that maximize the log of the likelihood *)
FindMaximum[{logL, σ > 0}, {{μ, μ0}, {σ, σ0}}]
{-606.0230945881103`, {μ -> 0.045999689972915744, σ ->  1.0000755520127822}}

For this data there is very little difference between the estimation procedures. A next step would be to assess the quality of the fit.

$\endgroup$
  • 2
    $\begingroup$ Here's a slightly compacted implementation: μ0 = data[[All, 1]].Normalize[data[[All, 2]], Total]; σ0 = Sqrt[((data[[All, 1]]^2).Normalize[data[[All, 2]], Total] - μ0^2)]; FindArgMax[{data[[All, 2]].Log[Erf[(# - μ - w/2)/(σ Sqrt[2]), (# - μ + w/2)/(σ Sqrt[2])]/2 &[data[[All, 1]]]], σ > 0}, {{μ, μ0}, {σ, σ0}}], where dot products were used for efficiency. Note as well the use of two-argument Erf[], which is much more numerically stable than subtracting two instances of Erf[]. $\endgroup$ – J. M. is away Aug 11 '17 at 4:21
2
$\begingroup$

Here is the fit to your data, based on a normal distribution:

FindFit[data, 
 a PDF[NormalDistribution[μ, σ], x], {a, μ, σ}, x]

(*

{a -> 37.2923, μ -> 0.134454, σ -> 0.834692}

*)

Show[
  ListPlot[data, PlotStyle -> Red],
  Plot[37.2923 PDF[NormalDistribution[0.134454, 0.834692], x], 
  {x, -3, 3}]
    ]

enter image description here

$\endgroup$
1
$\begingroup$

For what it's worth, here's a quick and easy approximation to the answer given by @Jim Baldwin:

EstimatedDistribution[
  Flatten[
    Cases[data, {x_, n_} :> ConstantArray[x, n]]
    ],
  NormalDistribution[m, s]
  ]
(* NormalDistribution[0.046, 1.00174] *)
$\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.