5
$\begingroup$

I have the following dataset (time vs. a quantity called $s$):

OPD = {{0, 0.0305}, {3, 0.14}, {6, 0.5}, {17, 4.02}, {19, 5.04}, {21, 5.52}, {24, 5.76}};

To the above, I have fitted a differential equation, $\dot{s} = \kappa (1 - s/\mu) s$, and obtained the parameters $s(0)$, $\kappa$, and $\mu$ as follows:

Bac = ParametricNDSolveValue[{s'[t] == kappa  (1 - s[t]/mu) s[t], s[0] == s0}, s, {t, 0, 24}, {s0, kappa, mu}];
J[s0_?NumericQ, kappa_?NumericQ, mu_?NumericQ] := Total@Map[(Bac[s0, kappa, mu][#[[1]]] - #[[2]])^2 &, OPD];
mini = NMinimize[{J[s0, kappa, mu], 0 < s0 < 2, 0 < kappa < 2, 0 < mu < 7}, {s0, kappa, mu}]

The fit looks like:

Show[ListPlot[OPD], Plot[Bac[s0, kappa, mu][t] /. mini[[2]], {t, 0, 24}, PlotRange -> All]]

enter image description here

I have just learned that Mathematica has an option called Weights, which can be used to give weights to the data points in order to take into account the uncertainties. It seems to me that, weights can be the standard deviations of the data points (?), which in my case are:

std = {0.0021, 0.0707, 0.0566, 0.5091, 0.2546, 0.8485, 0.3111};

I appreciate any help for how to include these standard deviations as weights to my data points in the course of fitting.

$\endgroup$
2
  • 1
    $\begingroup$ NonlinearModelFit[] has a Weight option that you may try. $\endgroup$ Commented Jan 29 at 13:18
  • 1
    $\begingroup$ BTW DSolve solves your DE analytically. DSolve[s'[t] == kappa (1 - s[t]/mu) s[t], s[t], t] gives your model as s[t] -> (E^(kappa t + mu C[1]) mu)/(-1 + E^(kappa t + mu C[1])) $\endgroup$ Commented Jan 29 at 13:26

2 Answers 2

5
$\begingroup$

Instead of NMinimize, use NonlinearModelFit which supports Weights. To define the fitting model, use ParametricNDSolveValue.

sol = ParametricNDSolveValue[{s'[t] == κ (1 - s[t]/μ) s[t],
    s[0] == s0}, s, {t, 0, 30}, {s0, κ, μ}]
NonlinearModelFit[OPD, sol[s0, κ, μ][t], {s0, κ, {μ, 6}}, t, Weights -> (1/std^2)]

Luckily, however, your differential equation has a well-known analytic solution, so you can use fit directly to the solution, provided by DSolveValue.

sol = DSolveValue[{s'[t] == κ (1 - s[t]/μ) s[t], s[0] == s0}, s[t], t]
(* (E^(t κ) s0 μ)/(-s0 + E^(t κ) s0 + μ) *)

OPD = {{0, 0.0305}, {3, 0.14}, {6, 0.5}, {17, 4.02}, {19, 5.04}, {21, 5.52}, {24, 5.76}};
std = {0.0021, 0.0707, 0.0566, 0.5091, 0.2546, 0.8485, 0.3111};
pts = MapThread[{#1[[1]], Around[#1[[2]], #2]} &, {OPD, std}];

fitNoWeights = NonlinearModelFit[OPD, sol, {s0, κ, {μ, 6}}, t];
fit = NonlinearModelFit[OPD, sol, {s0, κ, {μ, 6}}, t, Weights -> (1/std^2)];
Show[ListPlot[pts, PlotStyle -> Black], 
 Plot[{fitNoWeights[t], fit[t]}, {t, 0, 25}, 
  PlotLegends -> {"Without weights", "With weights"}]]

enter image description here

Note that when using weights, the fit will bind more tightly to the first three data points which have significantly smaller $\sigma$ than the last three points.

$\endgroup$
3
  • 1
    $\begingroup$ @Anovice, please see the edit. $\endgroup$
    – Domen
    Commented Jan 29 at 14:48
  • 1
    $\begingroup$ @Anovice, oh right, sorry, thanks! $\endgroup$
    – Domen
    Commented Jan 31 at 14:42
  • 1
    $\begingroup$ @Anovice, the order actually matters. "Positional" arguments come first, and named options (such as Weights->...) come at the end. $\endgroup$
    – Domen
    Commented Jan 31 at 17:19
4
$\begingroup$

Unfortunately, NMinimize does not have an option "Weights". Only fitting functions have it. However, you my multiply the residues by 1/std^2 to achieve the same effect like:

std = {0.0021, 0.0707, 0.0566, 0.5091, 0.2546, 0.8485, 0.3111};
OPD = {{0, 0.0305}, {3, 0.14}, {6, 0.5}, {17, 4.02}, {19, 5.04}, {21, 
    5.52}, {24, 5.76}};
Bac = ParametricNDSolveValue[{s'[t] == kappa   (1 - s[t]/mu)  s[t], 
    s[0] == s0}, s, {t, 0, 24}, {s0, kappa, mu}];
J[s0_?NumericQ, kappa_?NumericQ, mu_?NumericQ] := (1/std^2) . 
   Map[(Bac[s0, kappa, mu][#[[1]]] - #[[2]])^2 &, OPD];
mini = NMinimize[{J[s0, kappa, mu], 0 < s0 < 2, 0 < kappa < 2, 
   0 < mu < 7}, {s0, kappa, mu}]

![enter image description here

$\endgroup$
4
  • $\begingroup$ In a statistical sense, fitting using standard deviations is more appropriate. $\endgroup$ Commented Jan 29 at 10:42
  • $\begingroup$ @DanielHuber, I believe you have to divide by $\sigma^2$ (smaller $\sigma$ => larger contribution), namely (1/std^2) . Map[...]. Then, the fit also looks more convincing. $\endgroup$
    – Domen
    Commented Jan 29 at 12:03
  • 2
    $\begingroup$ I didn't check Daniel's code, so perhaps some parameter ranges need to be fixed. But still, when you are comparing by eye, you are not taking into account the error bars. Your fit is worse for the first three points, but they have very small uncertainties so the fit should be better there. $\endgroup$
    – Domen
    Commented Jan 29 at 12:37
  • $\begingroup$ @Domen You are 100% right. I corrected my mistake. Thank you for your hint. $\endgroup$ Commented Jan 29 at 13:00

Your Answer

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