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

Question

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?

---

Addition

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!

Answer 1

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.)

Answer 2

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.