# How to fit a list of points to a certain form given error bands?

I am trying to fit a list of points to the form F[0]/(1 - a*q^2 + b*q^4), where F[0] is the value at zero of a function, which can be considered an arbitrary constant that is greater than 0 for the sake of this question. Given three lists of points (one list that represents the mean function and two other that represent its upper and lower error bands), is there a way to find a fit to the mean list that will also produce uncertainties in a & b, from the form above? The three lists are as follows:

mean = {0.2500865832159757, 0.26529296193700797, 0.28194119360213926,
0.3002289030061626, 0.32039061666124463, 0.34270675020630514,
0.3675153635706632, 0.3952277382834709, 0.42634931999342446,
0.4615083323370265, 0.5014955856944889, 0.5473210060843121,
0.6002958069020696, 0.6621552007690517, 0.7352475215375661,
0.8228368038041404}

lower = {0.24459370549279338, 0.25748792043843843, 0.2715510835036353,
0.2869434924605513, 0.3038559622306384, 0.3225174644373506,
0.3432051944218582, 0.3662580214919165, 0.3920947371974243,
0.42123924031319887, 0.45435597032047914, 0.49230085694950415,
0.5361964272474037, 0.5875457504361626, 0.6484112063214175,
0.7217063441697118}

upper = {0.24313521571573743, 0.2595505081790414, 0.27759285870235867,
0.2974862871076955, 0.31949591470633976, 0.3439376891236028,
0.3711910031452933, 0.40171527094969367, 0.43607200112298017,
0.47495463915065683, 0.5192296094758698, 0.5699938651870915,
0.6286573955934119, 0.6970646053682822, 0.7776783645963563,
0.8738693713780176}


I've been looking specifically at trying to find the absolute error of the mean function for each point, then using that to find errors in a & b.

NonlinearModelFit can do this, but you have to translate your error bands into variances first. Assuming that your error bands are $$\pm1\sigma$$ (one standard deviation), the variances are

variance = ((upper - lower)/2)^2;


The nonlinear fit is then

f = NonlinearModelFit[mean,
F0/(1 - a*q^2 + b*q^4),
{{F0, 0.2}, {a, 0}, {b, 0}}, q,
VarianceEstimatorFunction -> (1 &), Weights -> 1/variance];


The best-fit parameters and their standard errors are

f["ParameterTable"]


(you can extract the numbers with f["ParameterTableEntries"]).

Let's compare a plot of the data points to the fit:

P = ListPlot[Around @@@ Transpose[{mean, Sqrt[variance]}], PlotStyle -> Black];
Q = Plot[Evaluate@f["MeanPredictionBands"], {q, 1, Length[mean]},
PlotStyle -> Blue, Filling -> {1 -> {2}}];
Show[P, Q]