8
$\begingroup$

I have the data:

data0 = {{0, 0.047}, {40, 0.071}, {80, 0.103}, {120, 0.169}, {160, 0.292}, {200, 0.498}, {240, 0.777}, {280, 1.000}, {320, 1.227}};

I want to fit a logistic function to it, so:

F = NonlinearModelFit[data0, A/(1 + E^(-(t - t0)/b)), {{A, 1.2}, {t0, 200}, {b, 60}}, t];
F["ParameterTable"]

The fit doesn't look like bad: enter image description here

Now, we are given the corresponding standard deviations as: $0.003, 0.008, 0.010, 0.012, 0.023, 0.024, 0.061, 0.135, 0.221$. If we use these as weights, one can write:

F = NonlinearModelFit[data0, A/(1 + E^(-(t - t0)/b)), {{A, 1.2}, {t0, 200}, {b, 60}}, t, Weights -> 1/{0.003, 0.008, 0.010, 0.012, 0.023, 0.024, 0.061, 0.135, 0.221}];
F["ParameterTable"]

But now the fit becomes much worse: enter image description here

What is the most effective way to use standard deviation in order to obtain more informative and better fit?

Improve this question
$\endgroup$
3
  • $\begingroup$ Where did you get the idea that including weights will result in a "better looking" fit? If the weights are correct, then you'll (tend) to get better estimates of the parameters and more appropriate estimates of the associated standard errors. $\endgroup$
    – JimB
    Commented Jul 6, 2021 at 23:05
  • $\begingroup$ @JimB Thanks for your comment. Doesn't better estimates of the parameters lead to better fit? $\endgroup$
    – NoOne
    Commented Jul 6, 2021 at 23:08
  • $\begingroup$ It does (generally) result in a better fit. But you are stating that it is a worse fit by the way it looks to you. That is what is wrong. $\endgroup$
    – JimB
    Commented Jul 7, 2021 at 0:31

1 Answer 1

Reset to default
9
$\begingroup$

One needs to have a model (a predictive part and an error structure) to know what to do with the weights. If the weights are appropriate for the model, then one will tend to get better estimates of the parameters in the sense that the precision of the estimates will improve from not using weights. But one cannot judge the improvement just by looking at the overall change in fit from not using weights to using weights.

Consider the following model:

$$y_i={a \over{1+e^{-(t_i-t_0)/b}}}+w_i \epsilon_i$$

where $y_i$ is the $i$-th response, $t_i$ is the $i$-th predictor, $a$, $b$, and $t_0$ are constants to be estimated, the errors ($\epsilon_i$) are all independent, and $\epsilon_i \sim N(0,\sigma^2)$ for the $n$ observations with $i=1,2,\ldots,n$.

The Mathematica code for fitting such a model is as follows:

(* Data and weights *)
data = {{0, 0.047}, {40, 0.071}, {80, 0.103}, {120, 0.169}, {160, 0.292}, 
  {200, 0.498}, {240, 0.777}, {280, 1.000}, {320, 1.227}};
w = {0.003, 0.008, 0.010, 0.012, 0.023, 0.024, 0.061, 0.135, 0.221};

(* Fit *)
nlm = NonlinearModelFit[data, a/(1 + E^(-(t - t0)/b)), {{a, 1.2}, {t0, 200}, {b, 60}}, t, 
  Weights -> 1/w^2]
mle = nlm["BestFitParameters"]
σ = nlm["EstimatedVariance"]^0.5

To check on whether using weights tends to result in better estimates for the parameters one can perform simulations. Suppose we use the estimated coefficients from your data as the "true" model. (But I'm reducing the size of the error variance to a tenth of the above estimate to avoid lack-of-convergence issues which can be especially problematic when the weights are ignored.)

(* Number of data points *)
n = Length[data];
(* Predictors *)
x = data[[All, 1]];
(* Scale down the variance to avoid lack-of-convergence issues *)
σ = 0.1 nlm["EstimatedVariance"]^0.5
(* Number of simulations *)
nsim = 200;

(* Generate mean prediction and errors *)
meanPrediction = a/(1 + Exp[-(x - t0)/b]);
errors = RandomVariate[NormalDistribution[0, σ], {nsim, n}];

(* Generate nsim datasets *)
datasim = meanPrediction + meanPrediction w # & /@ errors /. mle;
datasim = Transpose[{x, #}] & /@ datasim;

(* Now perform fits with and without weights *)
initialValues = {{a, (a /. mle)}, {b, (b /. mle)}, {t0, (t0 /. mle)}};
withWeights = {a, b, t0} /. NonlinearModelFit[#, a/(1 + Exp[-(t - t0)/b]),
  initialValues, t, Weights -> 1/w^2, MaxIterations -> 5000]["BestFitParameters"] & /@ datasim;
withoutWeights = {a, b, t0} /. NonlinearModelFit[#, a/(1 + Exp[-(t - t0)/b]),
  initialValues, t, MaxIterations -> 10000]["BestFitParameters"] & /@ datasim;

Now plot a smooth histogram of the estimates of each of the parameters with and without using weights:

Show[SmoothHistogram[{withoutWeights[[All, 1]], withWeights[[All, 1]]},
  Automatic, "PDF", PlotRange -> {{0, 15}, {0, 2}},
  Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"a", "Density"},
  PlotLegends -> {"Without weights", "With weights"}],
 ListPlot[{{a /. mle, 0}}, PlotStyle -> {{Red, PointSize[0.02]}}],
 PlotRangeClipping -> False]

Show[SmoothHistogram[{withoutWeights[[All, 2]], withWeights[[All, 2]]},
  Automatic, "PDF", Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"b", "Density"},
  PlotLegends -> {"Without weights", "With weights"}],
 ListPlot[{{b /. mle, 0}}, PlotStyle -> {{Red, PointSize[0.02]}}],
 PlotRangeClipping -> False]

Show[SmoothHistogram[{withoutWeights[[All, 3]], withWeights[[All, 3]]},
  Automatic, "PDF", PlotRange -> {{250, 500}, {0, .08}},
  Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"t0", "Density"},
  PlotLegends -> {"Without weights", "With weights"}],
 ListPlot[{{t0 /. mle, 0}}, PlotStyle -> {{Red, PointSize[0.02]}}],
 PlotRangeClipping -> False]

Smooth histograms of parameter estimates with and without weights

The red dots are the true values. One can see a great improvement by using weights (at least in this case where we know how the data was generated as used an estimation approach that matched the model). Much "tighter" estimates using weights.

Improve this answer
$\endgroup$
12
  • $\begingroup$ Many thanks for your detailed answered! Please give me time to understand it. Just a first question: do you mean $t_i$ in the line under "Consider the following model", where you have written instead $x_i$? $\endgroup$
    – NoOne
    Commented Jul 7, 2021 at 8:20
  • $\begingroup$ Another question: In your model, $\epsilon_i$, is the error and $w_i$ is a coefficient. But this $w_i$ is not the $w$ in your code, which is the standard deviation. Am I correct? Because I found notation confusing. $\endgroup$
    – NoOne
    Commented Jul 7, 2021 at 8:40
  • $\begingroup$ And also, where is this $w_i$ used in your code? $\endgroup$
    – NoOne
    Commented Jul 7, 2021 at 8:58
  • $\begingroup$ This list of $w_i$ values is defined in code (as w) right after data is defined and used with NonlinearModelFit function for both the original data and the simulated data with Weights -> 1/w^2. Also the $w_i$ are known coefficients that you supplied. $\endgroup$
    – JimB
    Commented Jul 7, 2021 at 15:34
  • 1
    $\begingroup$ Not exactly. To be blunt: the fit is only "visually" worse to someone who doesn't understand the effect of weights. Note that the "worsening" of the visual fit is generally at the points where the standard deviation is higher. And the fit is "visually better" for the points with more weight (i.e., smaller standard deviations). $\endgroup$
    – JimB
    Commented Jul 7, 2021 at 15:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.