# Chi-square minimization error

I'm trying to fit some data into a model using the standard procedure of minimizing the chi-square function using this code:

ClearAll["Global*"]
Needs["ErrorBarPlots"];
a0 = 1;
b0 = 2;

sol[a_, b_] :=
NDSolve[{x'[t] == -a*x[t] + Sin[t], y'[t] == a*x[t] - b*y[t]^2,
z'[t] == b*y[t] - Cos[t], x[0] == 1, y[0] == 0, z[0] == 0}, {x, y,
z}, {t, 0, 10}, MaxSteps -> Infinity];

F[t_] := x[t] + y[t]^2 + z[t]^3;

wsol[a_,
b_] := (wsol[a, b] =
NDSolve[{w'[t] - w[t] - ((F[t] /. sol[a, b])[[1]]) == 0,
w[0] == 0}, w, {t, 0, 10}]);

model[t_, a_, b_] := (w[t] /. wsol[a, b])[[1]] // Chop;

SeedRandom[1264645];
rangoT = Range[0, 9, 0.5];
ndat = Length[rangoT];

DATA = Table[{t,
model[t, a0, b0] + RandomVariate[NormalDistribution[0, .1]],
RandomReal[{.05, .1}]}, {t, rangoT}];

chi2[a_, b_] :=
Sum[((DATA[[i, 2]] - model[DATA[[i, 1]], a, b])/DATA[[i, 3]])^2, {i,
1, ndat}]

FindMinimum[chi2[a, b], {w, 0.9}, {b, 1.1}]


But when I run it I get error messages from the FindMinimum command related to NDSolve and ReplaceAll. Everything runs just fine until using that command. I was wondering why's not working and what should be changed.

• What error message do you get? Do you still get an answer? Dec 3, 2023 at 2:27
• The evaluation of DATA doesn't work. NDSolve throws copious errors at that point. Is that not the case for you? Dec 3, 2023 at 2:33
• @MarcoB that's my case... But everything appears when using the FindMinimum. Errors related to NDSolve and ReplaceAll... If you evaluate chi2 function, everything is ok Dec 3, 2023 at 3:47
• 1. So something is already wrong in the line DATA=..., please correct the description. 2. Use something rather than t in the Table will fix the issue. (Don't forget t is already used in those NDSolves and you haven't localize them! ) 3. {w, 0.9} in last line is obviously wrong. 4. You need _?NumericQ, see this post for more info: mathematica.stackexchange.com/a/26037/1871 Dec 3, 2023 at 4:28

1. solve for all variables simultaneously instead of solving for $$(x,y,z)$$ first and $$w$$ second,
2. use NonlinearModelFit instead of writing your own least-squares solver, and
3. going all the way up to $$t_{\text{max}}=10$$ seems a bit much and the differential equation seems to diverge.
tmax = 6;
sol = ParametricNDSolveValue[{
x'[t] == -a*x[t] + Sin[t],
y'[t] == a*x[t] - b*y[t]^2,
z'[t] == b*y[t] - Cos[t],
w'[t] - w[t] == x[t] + y[t]^2 + z[t]^3,
x[0] == 1, y[0] == 0, z[0] == 0, w[0] == 0},
{x, y, z, w}, {t, 0, tmax}, {a, b},
MaxSteps -> Infinity];
model[t_?NumericQ, a_?NumericQ, b_?NumericQ] := sol[a, b][[4]][t]

a0 = 1;
b0 = 2;
SeedRandom[1264645];
rangoT = Range[0, tmax, 0.5];
ndat = Length[rangoT];

DATA = Table[{t,
model[t, a0, b0] + RandomVariate[NormalDistribution[0, .1]],
RandomReal[{.05, .1}]}, {t, rangoT}];

result = NonlinearModelFit[DATA[[All, 1 ;; 2]], model[t, a, b],
{{a, a0}, {b, b0}}, t,
VarianceEstimatorFunction -> (1 &),
Weights -> 1/DATA[[All, 3]]^2]

result["ParameterTable"]


$$\begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline a & 1.0013 & 0.00103671 & 965.843 & \text{1.8409735010100922\grave{ }*{}^{\wedge}-28} \\ b & 1.9994 & 0.000490316 & 4077.78 & \text{2.423162678987459\grave{ }*{}^{\wedge}-35} \\ \end{array}$$

## update: analytic $$x(t)$$

Note that the differential equation for $$x(t)$$ can be solved analytically, which improves the stability and accuracy of the result:

DSolve[{x'[t] == -a*x[t] + Sin[t], x[0] == 1}, x[t], t] // FullSimplify
(*    {{x[t] -> ((2 + a^2) E^(-a t) - Cos[t] + a Sin[t])/(1 + a^2)}}    *)


Therefore we can define the problem more accurately as

x[a_, t_] = ((2 + a^2) E^(-a t) - Cos[t] + a Sin[t])/(1 + a^2);
sol = ParametricNDSolveValue[{
y'[t] == a*x[a, t] - b*y[t]^2,
z'[t] == b*y[t] - Cos[t],
w'[t] - w[t] == x[a, t] + y[t]^2 + z[t]^3,
y[0] == 0, z[0] == 0, w[0] == 0},
{y, z, w}, {t, 0, tmax}, {a, b}, MaxSteps -> Infinity];
model[t_?NumericQ, a_?NumericQ, b_?NumericQ] := sol[a, b][[3]][t]

• Wow! Thanks... The use of ParametricNDSolveValue and NonlinearModelFit solved those issues. What I was wondering is that if could plot the contour plots for sigma amplitudes using the results of NonlinearModelFit and the definition of chi-square function used above. Dec 3, 2023 at 23:50
• @Syn1110 Graphics[result["ParameterConfidenceRegion"], Frame -> True, FrameLabel -> {a, b}] shows the confidence ellipse of the parameters. Please have a look at result["Properties"] to see all the goodies you can get out of the fit result. Dec 4, 2023 at 7:43
• Thanks a lot! That only shows confidence region at 1sigma amplitude right? Is there a way to obtain 2sigma and 3sigma? Dec 5, 2023 at 0:30

The problem is caused by "NDSolve". It should localize the dummy variable "t", but fails to do so. Consider:

t = 0;
NDSolve[{f'[t] == f[t], f[0] == 1}, f, {t, 0, 10}]


This gives the error message:

I think this is a bug.

• No, it's not. NDSolve simply doesn't have attribute like HoldAll that adjusts the evaluation order, this is by design. Dec 3, 2023 at 9:26
• @xzczd Do you see any reason why this should be by design? Dec 3, 2023 at 9:30
• Otherwise things like this will happen on Solve/DSolve/NDSolve: opt=PlotStyle->Red; Plot[x,{x,0,1},opt] Dec 3, 2023 at 9:51
• Actually in early versions, even domain = {x, 0, 1}; Plot[x, domain] is enough to cause trouble: i.stack.imgur.com/LpOl8.png This is all because of the HoldAll. If Solve, etc. behaves similarly, it'll be rather inconvenient. (Just imagine the world that code like eq = {a + b == 0, a - 2 b == 3}; vars = {a, b}; Solve[eq, vars] results in failure! ) Dec 3, 2023 at 11:03