# Fitting in mathematica when dealing with non-gaussian noise and errors on the data points

How is it best to deal with fitting in Mathematica when the data has a non-normal or non-Gaussian noise on top of it? This could also be phrased as how to fit correctly in Mathematica when your data has asymmetrical error bars?

Here is a small example. Consider a simple function of the the form $$f(t)= A \sin(2\pi\nu_0 t) + \xi(t)$$ where $$\xi(t)$$ is a noise parameter drawn from a distribution of our choice. In Mathematica we set this up as

f[A_, \[Nu]0_, t_] := A Sin[2 \[Pi] \[Nu]0 t]

f0 = 10;
Amax = 1;
dt = 0.005;
t0 = 0;
tf = 0.4;


And the data is generated and noised with a Gaussian distribution by

\[Sigma]Gauss = 0.5;
\[Mu]Gauss = 0;

GaussianNoisedData =
Table[
{
t,
Around[f[Amax, f0, t] +  RandomVariate[NormalDistribution[\[Mu]Gauss, \[Sigma]Gauss]], \[Sigma]Gauss]
},
{t, t0, tf, dt}
];


Seeing as we know the distribution parameters of the Gaussian noise we can define the error bars and weights for fitting by the standard deviation, $$\sigma_\rm{G}$$, of the Normal (Gaussian) distribution

GaussianNoisedDataFit =
NonlinearModelFit[
GaussianNoisedData /. a_Around:>a["Value"],
f[A, \[Nu]0, t],
{{A, 1}, {\[Nu]0, 10}}, t,
Weights->1/ConstantArray[\[Sigma]Gauss, Length[NoiselessData]]^2,
VarianceEstimatorFunction->(1&)
]


Plotting the data and the fit gives us this:

For the case of Gaussian statistics and noise defining error bars and weights is easy. We can simply say $$1/\sigma^{2}$$ where $$\sigma$$ can be an error bar or a confidence interval. In this case I just use the standard deviation of the Gaussian noise I add to the data -- so a $$68.27\%$$ confidence interval.

As I understand it though, the fact that weights are usually $$1/\sigma^{2}$$ is a consequence of Gaussian statistics.

## What do we do for asymmetrical error bars and confidence intervals related to asymmetrical probability density function (PDF)?

For example lets replace the noise generation wit that of a Rayleigh distribution (I choose this as it is somewhat related to the Gaussian, it is asymmetric and has only one parameter)

\[Sigma]Rayleigh = 0.5;

RayleighNoisedData =
Table[
{
t,
Around[f[Amax, f0, t] +  RandomVariate[RayleighDistribution[\[Sigma]Rayleigh]], {0.336, 1.041}]
},
{t, t0, tf, dt}
];


Note that I have given it an error bar in the form of Around[f,{0.336, 1.041}]. I determine these values by finding the $$68\%$$ confidence interval around the mean of a Rayleigh distribution which is $$\sigma \sqrt{\pi/2}$$.

Given that I don't know how to deal with the error bars in the context of a weighted fit, I simply perform an unweighted fit

RayleighNoisedDataFit =
NonlinearModelFit[
RayleighNoisedData /. a_Around:>a["Value"],
f[A, \[Nu]0, t],
{{A, 1}, {\[Nu]0, 10}}, t
]


Plotting the data with the fit gives

Note the asymmetrical error bars.

## Summary and questions

When dealing with non-Gaussian statistics and errors, are there options in NonlinearModelFit to allow weighting with asymmetrical error bars and weights?

Can we give Mathematica a distribution and perform a fit based on a maximum likelihood estimation (MLE), if the noise or distribution of errors on data points is known?

Based on some input from a very good answer bellow, I think the example I gave above was probably quite bad for the question I was asking! For my example, as indicated in the answer below -- weighting is inappropriate!

It would be nice to see an example where data with asymmetric error bars/bands can be fitted using the error bars/bands as a weight. I haven't been able to find an example on how to do this in Mathematica!

• I think this would be a better question for CrossValidated as there are many, many existing statistical approaches. There's no need to re-invent things here. Also, error bars are usually to show sampling/measurement error on individual sample points. If you're interested in estimating the parameters of a function, then it is some sort of error band about the curve that would be of interest. And as you seem to suggest, asymmetrical error bars - or asymmetrical error bands - and many times the way to go.
– JimB
Dec 28, 2019 at 6:19
• @JimB ". If you're interested in estimating the parameters of a function, then it is some sort of error band about the curve that would be of interest. And as you seem to suggest, asymmetrical error bars - or asymmetrical error bands " Yes!! exactly this is what I am interested in!! My question is really about how to weight a fit in Mathematica when dealing with asymmetrical error bars/bands.
– Q.P.
Dec 28, 2019 at 15:28
• I forgot to mention another issue with your example: All of the weights are equal which eliminates the need (if there really is a need) to have an analysis that accounts for the weights. Only if the "weights" vary among the observations might weighting be useful.
– JimB
Dec 28, 2019 at 16:53
• @JimB Yes that did occur to me actually -- its a pretty useless example all round. I'm surprised I haven't received any close requests on it actually. But still, if you are able to show a small example of weighting with an asymmetric error bar or error band, that would be extremely helpful to me and probably others interested in this type of fit.
– Q.P.
Dec 28, 2019 at 16:58
• If you have a particular example in mind, then I'd be happy to help. In the meantime: (1) If you know the model and the form of the error distribution (whether the errors are identical or not), then I'd use maximum likelihood or a Bayesian approach. (2) If you know the form of the model but not the form of the error distribution, then robust regression techniques might do what you want. (3) If you don't know the form of the model or the form of the error distribution, then you might want to use quantile regression or gams or nonparametric regression. So many techniques, so little time.
– JimB
Dec 28, 2019 at 22:36

When you have just a single additive error term with a known distributional form, you should avoid thinking about "weights" (except for some distributions where using weights and least squares gives you good starting values for troublesome datasets and models).

For the example you give using maximum likelihood is pretty straightforward. If you have some knowledge about the parameters of interest that can be characterized with a probability distribution, then a Bayesian approach would be your best bet.

Here is a maximum likelihood approach:

(* Define function *)
f[A_, ν0_, t_] := A Sin[2 π ν0 t]

(* Generate some data *)
f0 = 10;
Amax = 1;
dt = 0.005;
t0 = 0;
tf = 0.4;
σRayleigh = 4;
SeedRandom[12345];
RayleighNoisedData = Table[{t, f[Amax, f0, t] +
RandomVariate[RayleighDistribution[σRayleigh]]}, {t, t0, tf, dt}];

(* Find log of the likelihood and maximize it *)
logL = LogLikelihood[RayleighDistribution[σ], (#[[2]] - f[A, ν0, #[[1]]]) & /@ RayleighNoisedData];
sol = FindMaximum[{logL[[1, 1, 1]], σ > 0}, {{A, Amax}, {ν0, f0}, {σ, σRayleigh}}]
(* {-176.659, {A -> 0.790322, ν0 -> 10.3217, σ -> 3.83458}} *)

(* Plot results *)
Show[ListPlot[RayleighNoisedData],
Plot[{f[A, ν0, t] /. {A -> Amax, ν0 -> f0},
f[A, ν0, t] /. sol[[2]]}, {t, t0, tf},
PlotLegends -> {"True", "Estimated"}], PlotRange -> All]


In this case there's just no need for symmetric or asymmetric error bands on the data points. If there is sampling/measurement error in addition to the random deviations from the underlying curve, then you need to use software that accounts for multiple error sources (i.e., mixed models). Mathematica does not yet offer mixed model functions (although one can obtain estimates of mixed model parameters with Mathematica, there are no higher level functions equivalent to NonlinearModelFit to do so).

(I've left off error bars (symmetric and/or asymmetric) because I don't see that they help in interpretation in this case. Also, I've not included essential confidence bands around the estimated curve because I've run out of time tonight.)

• +1. A clear, concise and sound answer as usual, JimB.
– ciao
Dec 28, 2019 at 7:55
• Thanks for the answer here @JimB! I think my example was not very good for the question I was asking -- I guess what I am really asking is: How when dealing with an asymmetrical error-bar or error-band on some data can we use this as a weight in a fit? I can't find any documentation. I have made a modification to my question, in light of your explanations here.
– Q.P.
Dec 28, 2019 at 14:47

I found a very interesting paper that discusses dealing with asymmetrical error bars and how to use them as weights. I'm including it here as it may be useful for others in the community, but I will now accept @JimB's answer as he provided a very nice example MLE fitting and, as usual, some excellent comments on the original question.

• Thanks for the link. I must say the terminology used is not at all standard statistical terminology and there are no non-Gaussian distributions mentioned (except for a reference to Poisson which is discrete rather than continuous so I don't know why that's mentioned). It will take a bit of studying on my part.
– JimB
Jan 6, 2020 at 16:53
• @JimB I'd be curious to get your opinions once you've studied up on it.
– Q.P.
Jan 6, 2020 at 20:12