Non-linear least squares
Authors: Housam Binous and Ahmed Bellagi
Non-linear least squares
Authors: Housam Binous and Ahmed Bellagi
Authors: Housam Binous and Ahmed Bellagi
Hypothetical data generated by using the function a1 Log[x+a2]+Exp[a3/(x+1)]
with a1=2.5, a2=1.2 and a3=4 and a random number generator
with a1=2.5, a2=1.2 and a3=4 and a random number generator
In[]:=
data=Table[2.5Log[x+1.2]+Exp[4/(x+1)]+0.1Random[],{x,0,10,0.1}]
Out[]=
{55.081,38.6349,28.9293,22.7891,18.629,15.8154,13.6865,12.2046,11.0259,10.0894,9.41176,8.82455,8.44299,8.04677,7.68875,7.50915,7.25189,7.13385,6.95223,6.87158,6.76455,6.70323,6.63231,6.50514,6.48115,6.46556,6.40128,6.37902,6.4245,6.37849,6.39753,6.34447,6.32425,6.33216,6.33718,6.32202,6.34212,6.38844,6.35906,6.38297,6.36147,6.36078,6.37563,6.46567,6.45568,6.43673,6.45637,6.51883,6.48779,6.54771,6.60205,6.56339,6.56882,6.66136,6.58733,6.69594,6.71891,6.70361,6.72702,6.75537,6.76471,6.81119,6.77355,6.78887,6.83173,6.89195,6.8534,6.89345,6.96175,6.95701,6.91547,7.01019,7.0271,6.98039,7.0028,7.05476,7.08508,7.06488,7.08058,7.14334,7.17022,7.18499,7.19703,7.25126,7.1887,7.27691,7.24654,7.30542,7.31997,7.36794,7.39607,7.33484,7.40099,7.37356,7.46305,7.42686,7.49581,7.43474,7.5214,7.55278,7.51085}
In[]:=
xvalues=Table[x,{x,0,10,0.1}]
Out[]=
{0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.,4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5.,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.,6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.,7.1,7.2,7.3,7.4,7.5,7.6,7.7,7.8,7.9,8.,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8,9.9,10.}
In[]:=
a1=.;a2=.;a3=.;
In[]:=
dd=Table[{xvalues[[i]],data[[i]]},{i,1,Length[data]}]
Out[]=
{{0.,55.081},{0.1,38.6349},{0.2,28.9293},{0.3,22.7891},{0.4,18.629},{0.5,15.8154},{0.6,13.6865},{0.7,12.2046},{0.8,11.0259},{0.9,10.0894},{1.,9.41176},{1.1,8.82455},{1.2,8.44299},{1.3,8.04677},{1.4,7.68875},{1.5,7.50915},{1.6,7.25189},{1.7,7.13385},{1.8,6.95223},{1.9,6.87158},{2.,6.76455},{2.1,6.70323},{2.2,6.63231},{2.3,6.50514},{2.4,6.48115},{2.5,6.46556},{2.6,6.40128},{2.7,6.37902},{2.8,6.4245},{2.9,6.37849},{3.,6.39753},{3.1,6.34447},{3.2,6.32425},{3.3,6.33216},{3.4,6.33718},{3.5,6.32202},{3.6,6.34212},{3.7,6.38844},{3.8,6.35906},{3.9,6.38297},{4.,6.36147},{4.1,6.36078},{4.2,6.37563},{4.3,6.46567},{4.4,6.45568},{4.5,6.43673},{4.6,6.45637},{4.7,6.51883},{4.8,6.48779},{4.9,6.54771},{5.,6.60205},{5.1,6.56339},{5.2,6.56882},{5.3,6.66136},{5.4,6.58733},{5.5,6.69594},{5.6,6.71891},{5.7,6.70361},{5.8,6.72702},{5.9,6.75537},{6.,6.76471},{6.1,6.81119},{6.2,6.77355},{6.3,6.78887},{6.4,6.83173},{6.5,6.89195},{6.6,6.8534},{6.7,6.89345},{6.8,6.96175},{6.9,6.95701},{7.,6.91547},{7.1,7.01019},{7.2,7.0271},{7.3,6.98039},{7.4,7.0028},{7.5,7.05476},{7.6,7.08508},{7.7,7.06488},{7.8,7.08058},{7.9,7.14334},{8.,7.17022},{8.1,7.18499},{8.2,7.19703},{8.3,7.25126},{8.4,7.1887},{8.5,7.27691},{8.6,7.24654},{8.7,7.30542},{8.8,7.31997},{8.9,7.36794},{9.,7.39607},{9.1,7.33484},{9.2,7.40099},{9.3,7.37356},{9.4,7.46305},{9.5,7.42686},{9.6,7.49581},{9.7,7.43474},{9.8,7.5214},{9.9,7.55278},{10.,7.51085}}
finding the values of a1, a2 and a3 using the built-in command FindFit
In[]:=
FindFit[dd,a1Log[x+a2]+Exp[a3/(x+1)],{{a1,2.2},{a2,1.1},{a3,3.7}},x]
Out[]=
{a12.51839,a21.23937,a33.99895}
Using Reap and Sow we get the results of the first five iterations that FindFit has determined
In[]:=
sol=Reap[FindFit[dd,a1Log[x+a2]+Exp[a3/(x+1)],{{a1,2.2},{a2,1.1},{a3,3.7}},x,StepMonitorSow[{a1,a2,a3}]]]
Out[]=
{{a12.51839,a21.23937,a33.99895},{{{2.52889,1.15571,4.05127},{2.51888,1.23493,4.00028},{2.51839,1.23936,3.99895},{2.51839,1.23937,3.99895},{2.51839,1.23937,3.99895}}}}
In[]:=
z=Flatten[Transpose[sol[[2]]],1]
Out[]=
{{2.52889,1.15571,4.05127},{2.51888,1.23493,4.00028},{2.51839,1.23936,3.99895},{2.51839,1.23937,3.99895},{2.51839,1.23937,3.99895}}
In[]:=
s=Join[{{"Iteration",a1,a2,a3}},Transpose[Join[{{1,2,3,4,5}},Transpose[z]]]]
Out[]=
{{Iteration,a1,a2,a3},{1,2.52889,1.15571,4.05127},{2,2.51888,1.23493,4.00028},{3,2.51839,1.23936,3.99895},{4,2.51839,1.23937,3.99895},{5,2.51839,1.23937,3.99895}}
In[]:=
TableForm[s]
Out[]//TableForm=
Iteration | a1 | a2 | a3 |
1 | 2.52889 | 1.15571 | 4.05127 |
2 | 2.51888 | 1.23493 | 4.00028 |
3 | 2.51839 | 1.23936 | 3.99895 |
4 | 2.51839 | 1.23937 | 3.99895 |
5 | 2.51839 | 1.23937 | 3.99895 |
Let us do the same calculation using the Gauss–Newton algorithm for our non-linear least squares test problem.
The normal equations are written in matrix notation as: (J) Δβ = Δy
For details see: Non-linear least squares wikipedia
The normal equations are written in matrix notation as: (
T
J
T
J
For details see: Non-linear least squares wikipedia
We choose an initial guess to start the numerical procedure: a1=2.2;a2=1.1;a3=3.7;
In[]:=
a1=.;a2=.;a3=.;
In[]:=
f[a1_,a2_,a3_][x_]:=a1Log[x+a2]+Exp[a3/(x+1)]
In[]:=
jacobian1[a1_,a2_,a3_][x_]:=Evaluate[Outer[D,{f[a1,a2,a3][x]},{a1,a2,a3}]]
In[]:=
jacobian1[a1,a2,a3][x]
Out[]=
Log[a2+x],,
a1
a2+x
a3
1+x
1+x
First iteration
In[]:=
a1=2.2;a2=1.1;a3=3.7;
In[]:=
b=data-Map[f[a1,a2,a3],xvalues]
Out[]=
{14.4241,9.33976,6.52104,4.82748,3.68372,2.99825,2.41945,2.09638,1.80264,1.55424,1.41968,1.26642,1.23527,1.12445,1.00054,1.01408,0.916869,0.931874,0.861094,0.872841,0.842808,0.845536,0.827692,0.74426,0.756046,0.769371,0.728103,0.723731,0.782674,0.746244,0.77149,0.721679,0.702105,0.70837,0.70974,0.689141,0.70224,0.740177,0.70118,0.714379,0.681201,0.667984,0.669553,0.745631,0.721098,0.687072,0.691166,0.737683,0.690338,0.733634,0.771076,0.715275,0.70336,0.778381,0.686673,0.777476,0.782542,0.749242,0.754581,0.764819,0.75599,0.784285,0.728444,0.725565,0.750227,0.792258,0.73555,0.75746,0.807668,0.784871,0.725332,0.802095,0.801112,0.736568,0.741212,0.775469,0.788156,0.750401,0.748611,0.793967,0.803528,0.801048,0.795925,0.833076,0.753522,0.824824,0.777637,0.819777,0.817686,0.849093,0.860759,0.783142,0.833003,0.78937,0.862742,0.810533,0.863548,0.78663,0.857528,0.873231,0.815726}
In[]:=
J=Flatten[Map[jacobian1[a1,a2,a3],xvalues],1]
Out[]=
{{0.0953102,2.,40.4473},{0.182322,1.83333,26.2673},{0.262364,1.69231,18.1925},{0.336472,1.57143,13.2472},{0.405465,1.46667,10.0381},{0.470004,1.375,7.8554},{0.530628,1.29412,6.31228},{0.587787,1.22222,5.18538},{0.641854,1.15789,4.33954},{0.693147,1.1,3.68959},{0.741937,1.04762,3.17991},{0.788457,1.,2.7731},{0.832909,0.956522,2.44333},{0.875469,0.916667,2.1723},{0.916291,0.88,1.94682},{0.955511,0.846154,1.75718},{0.993252,0.814815,1.5961},{1.02962,0.785714,1.45808},{1.06471,0.758621,1.33885},{1.09861,0.733333,1.2351},{1.1314,0.709677,1.14422},{1.16315,0.6875,1.06412},{1.19392,0.666667,0.993123},{1.22378,0.647059,0.92987},{1.25276,0.628571,0.873244},{1.28093,0.611111,0.822325},{1.30833,0.594595,0.776347},{1.335,0.578947,0.734671},{1.36098,0.564103,0.696758},{1.38629,0.55,0.662153},{1.41099,0.536585,0.630467},{1.43508,0.52381,0.601368},{1.45862,0.511628,0.574571},{1.4816,0.5,0.549828},{1.50408,0.488889,0.526926},{1.52606,0.478261,0.505678},{1.54756,0.468085,0.485921},{1.56862,0.458333,0.467513},{1.58924,0.44898,0.450326},{1.60944,0.44,0.434251},{1.62924,0.431373,0.419187},{1.64866,0.423077,0.405048},{1.66771,0.415094,0.391754},{1.6864,0.407407,0.379237},{1.70475,0.4,0.367433},{1.72277,0.392857,0.356286},{1.74047,0.385965,0.345746},{1.75786,0.37931,0.335765},{1.77495,0.372881,0.326304},{1.79176,0.366667,0.317323},{1.80829,0.360656,0.30879},{1.82455,0.354839,0.300673},{1.84055,0.349206,0.292944},{1.8563,0.34375,0.285576},{1.8718,0.338462,0.278546},{1.88707,0.333333,0.271832},{1.90211,0.328358,0.265414},{1.91692,0.323529,0.259274},{1.93152,0.318841,0.253395},{1.94591,0.314286,0.247761},{1.96009,0.309859,0.242358},{1.97408,0.305556,0.237172},{1.98787,0.30137,0.232192},{2.00148,0.297297,0.227404},{2.0149,0.293333,0.2228},{2.02815,0.289474,0.218369},{2.04122,0.285714,0.214101},{2.05412,0.282051,0.209989},{2.06686,0.278481,0.206024},{2.07944,0.275,0.202198},{2.09186,0.271605,0.198505},{2.10413,0.268293,0.194938},{2.11626,0.26506,0.191491},{2.12823,0.261905,0.188158},{2.14007,0.258824,0.184934},{2.15176,0.255814,0.181814},{2.16332,0.252874,0.178792},{2.17475,0.25,0.175866},{2.18605,0.247191,0.173029},{2.19722,0.244444,0.170278},{2.20827,0.241758,0.16761},{2.2192,0.23913,0.165021},{2.23001,0.236559,0.162508},{2.24071,0.234043,0.160067},{2.25129,0.231579,0.157695},{2.26176,0.229167,0.15539},{2.27213,0.226804,0.153149},{2.28238,0.22449,0.150969},{2.29253,0.222222,0.148848},{2.30259,0.22,0.146783},{2.31254,0.217822,0.144773},{2.32239,0.215686,0.142816},{2.33214,0.213592,0.140909},{2.34181,0.211538,0.13905},{2.35138,0.209524,0.137238},{2.36085,0.207547,0.135471},{2.37024,0.205607,0.133748},{2.37955,0.203704,0.132067},{2.38876,0.201835,0.130425},{2.3979,0.2,0.128823},{2.40695,0.198198,0.127259}}
In[]:=
Length[J]
Out[]=
101
In[]:=
Transpose[J].b
Out[]=
{142.117,114.689,1154.61}
Second iteration
Third iteration
Plotting the hypothetical experimental data and the model