# Non linear fit/Find distribution parameters for large gradient data;

searching around I've found that this is a common problem with fitting but I haven't found a work-able solution.

I have a symmetric data set that I've manipulated so that it is a probability distribution, ie. integral sums to 1, but the peak occurs between -0.002 and 0.002. Using FindDistributionParameters returns an error but I can fit it with a normal distribution using:

data = ToExpression@ImportString[Import["https://pastebin.com/raw/cVExSjYq"], "Text"]

{\[Mu], \[Sigma]} = (NonlinearModelFit[
data,
1/Sqrt[2*\[Pi]*\[Sigma]^2]*E^(-(x - \[Mu])^2/(
2*\[Sigma]^2)), {\[Mu], \[Sigma]}, x][
"BestFitParameters"][[All, 2]]) /. {x_, y_} -> {x, Abs[y]}


outputs:

{1.24156*10^-20, 0.000250996}


The distribution doesn't closely match, hence I'm trying to use other distributions but fitting it with

FindDistributionParameters[data,StudentTDistribution[0, \[Sigma], \[Nu]], {{\[Sigma], 0.00025}, {\[Nu], 0.5}}]


returns the error:

"One or more data points are not in support of the process or distribution StudentTDistribution[0,\[Sigma],\[Nu]]."


I get the same error using WeibullDistribution too and trying to manually find the NonlinearModelFit just returns {1,1,1} for {alpha,beta,mu}.

I've tried removing zero values using Select[data,#[[2]]!=0&] but that doesn't help.

Plotted data with StudentT[0,0.00025,0.5] in blue and NormalDistribution[0,0.00025] in red is shown here:

Any recommendations?

• What are you trying to do with the obtained fit? – Anton Antonov Aug 2 '17 at 0:45
• @AntonAntonov I'm trying to get a reasonable comparison metric for various trials that are returning similar distributions. I also want to generate a RandomVariate[] to see how a virtual data set responds some other analysis. – Andrew Stewart Aug 2 '17 at 19:09
• How sure are you that the distribution/function is defined at x=0? Or that it's continuously differentiable? It looks to me to have the form Abs[a x^-b], or something similar, ie: undefined at x=0. Maybe a GammaDistribution (flipped over the vertical axis)? – aardvark2012 Aug 3 '17 at 1:16
• Part of your manipulations seems to have made the manipulated data perfectly symmetric. Was that on purpose? Was the data duplicated by reflecting it around zero? – JimB Aug 3 '17 at 2:34
• @aardvark2012 Yes the problem set is defined at 0 and should be continuously differentiable. This data set here was created by mirroring about x=0, but since it is finite at x=0 I didn't find ExponentialDistribution[] or something similar returning appropriate values. I'm not extremely well versed in statistics but I'll take a look at the GammaDistribution. Thanks! – Andrew Stewart Aug 3 '17 at 18:43

This, too, is an extended comment (and maybe just an repeat of @aardvark2012's and @AntonAntonov's comments).

The data set you present is a relative frequency distribution (and the actual sample size is lost). The integral doesn't sum to 1 as there is no integral to sum or integrate. What you do have is 0.0001*Total[data[[All,2]]] equaling 1.0.

Using "least squares" to fit a probability distribution is usually not a recommended approach as it usually makes no sense to do so. One of the many reasons is that the error variance is not constant: the tail areas that drop to zero would have to have very small error variances and that is not one of the assumptions in the least squares fitting process you've used.

(Note that for any real values of μ and positive real values of σ the integral Integrate[E^(-(x - μ)^2/(2*σ^2))/Sqrt[2*π*σ^2],{x,-∞,∞} will always equal 1.)

In short, you are mixing up regression and fitting probability distributions.

You have a regression situation if you want to predict the data[[All,2]] values from the data[[All,1]] values and believe that the curve form is the same as a probability density function. This, however, does not impart any probabilistic properties to the fit.

If the raw data consists of a random sample from some probability distribution, then you are better off to use that raw data either with FindDistributionParameters if you really know the form of the probability density function. If you don't know the form and have a large enough sample size, then using SmoothKernelDistribution should be considered.

• Thanks for the response. Yes, the data I presented is heavily modified from what was measured. I have the original frequency histogram but it has irregular intervals, (corrected using a method from Scott & Scott, 2008) but I wasn't able to run that dataset through FindDistributionParameters either. Returned similar error. I'll have to read up on why the least squares is inappropriate. I don't really follow the error variance issue. I'll keep reading. Thanks! – Andrew Stewart Aug 2 '17 at 19:14
• Here is a link to the Scott & Scott article. One really needs the original data which gets you the associated sample size. – JimB Aug 3 '17 at 2:04
• the original sample size for this set of data is about 12500 samples. I'm not an expert in stats but I've read enough to know that many of the fit tests break down with >5000 samples. I would put up the file but it gets awkward working with files this large. Does this change anything of great significance? – Andrew Stewart Aug 3 '17 at 18:49

Not an answer, but an extended comment.

The first point I'd make is that the data slot in FindDistributionParameters is not the same as in NonlinearModelFit -- it's just a list of outcomes, not 2D points (which explains the error message).

Second, there are a lot of data points lying on the x-axis. And all of them carry the same weight as the few that give the distribution its shape. This may distort the results returned by NonlinearModelFit.

Third, your data is very spikey, and I don't think either of the models you've tried are, by themselves, spikey enough to get a good fit.

Here's a slightly modified version of the StudentTDistribution:

model[x_, \[Alpha]_, \[Sigma]_, \[Nu]_] := \[Alpha] (\[Nu]/(\[Nu] +
x^2/\[Sigma]^2))^((1 + \[Nu])/2)/(
Sqrt[\[Nu]] \[Sigma] Beta[\[Nu]/2, 1/2])

Manipulate[
Show[ListPlot[data],
Plot[model[x, \[Alpha], \[Sigma], \[Nu]], {x, -0.003, 0.003}, PlotRange -> All],
PlotRange -> {{-0.003, 0.003}, {0, 2000}}],
{{\[Alpha], 1.36201}, 10^-5, 3}, {{\[Sigma], 0.000149}, 10^-5, 0.0005},
{{\[Nu], 0.284}, 10^-5, 1}
]


Unfortunately, FindFit and NonlinearModelFit both seem pretty unhelpful using this model, even with these initial values for the parameters.

I'll try something else if I get a chance.

• Thanks. I'll take a look and try and quickly scheme something to quantify the closeness. After reading @JimBaldwin's comment though not sure I'm on the right line of thinking. I'll keep reading. Thanks! – Andrew Stewart Aug 2 '17 at 19:08