Message Boards Message Boards

[?] Fit noisy data as a Least squares problem?

Posted 4 years ago

I want to fit a set of noisy data, obtained from the equation fm with added noise, to the same equation (with y=0; t=1).

 Dif = 0.0000001; K0 = 0.5; z = 0; Intens = 25000;
    angulo1 = 45*Pi/180;
    lsup1 = Cos[angulo1]; linf\[Alpha] = 0.0;
    m1 = -Tan[angulo1];

fm = T145[x_, y_, t_, b_, l_, d_] := 
  Intens/(2*\[Pi]*K0)*
   NIntegrate[
    Sqrt[1 + m1^2]*
     Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m1*\[Alpha] \
- d))^2]/Sqrt[
        Dif*4*t]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - \
(m1*\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha],
      l*lsup1}]

Here data with added noise and then fitting:

Block[{y = 0, t = 1, b = 0.001, l = 0.001, d = 0.0001}, 
  data = Table[{x, 
     fm[x, y, t, b, l, d] + 
      Random[NormalDistribution[0, 0.1]]}, {x, -0.0015, 0.002, 
     0.000015}]];

model = Intens/(2*\[Pi]*K0)*
   NIntegrate[
    Sqrt[1 + m1^2]*
     Erfc[Sqrt[(x - \[Alpha])^2 + (-\[Beta])^2 + (z - (m1*\[Alpha] - 
             d))^2]/Sqrt[
        Dif*4*1]]/(Sqrt[(x - \[Alpha])^2 + (-\[Beta])^2 + (z - (m1*\
\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha], 
     l*lsup1}];

fit = FindFit[data, model, {b, l, d}, x, 
  Method -> "LevenbergMarquardt"]

But the results are {b->1.,l->1.,d->1.}, that are patently wrong. Where is the problem? Thanks to all.

POSTED BY: Lorenzo Fuggiano
4 Replies

Not sure how your code actually produced a result since it produced some errors for me.

NIntegrate only works when the boundaries are numerical. Therefore when you use it as a fit function you have to be sure the function only evaluates if the input parameters are numerical.

You can use the same function for your data generation and your fitting.

Also with minimization problems the initial fit value should be close to the expected values since you are using local minimization "LevenbergMarquardt". Therefore give initial values close to the expected values.

fm[x_, y_, t_, b_, l_, d_] /; 
  NumberQ[x] && NumberQ[b] && NumberQ[l] && NumberQ[d] := 
 Intens/(2*\[Pi]*K0)*
  NIntegrate[
   Sqrt[1 + m1^2]*
    Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m1*\[Alpha] -
               d))^2]/
       Sqrt[Dif*4*
         t]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m1*\
\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha], 
    l*lsup1}]

Block[{y = 0, t = 1, b = 0.0015, l = 0.0011, d = 0.00013}, 
  data = Table[{x, 
     fm[x, y, t, b, l, d] + 
      Random[NormalDistribution[0, 0.1]]}, {x, -0.0015, 0.002, 
     0.000015}]
  ];

fit = FindFit[data, 
  fm[x, 0, 1, b, l, d], {{b, 10.^-3}, {l, 10.^-3}, {d, 10.^-4}}, x, 
  Method -> "LevenbergMarquardt"]

Show[ListLinePlot[data, PlotStyle -> Gray, PlotRange -> All], 
 Plot[fm[x, 0, 1, b, l, d] /. fit, {x, -0.0015, 0.002}, 
  PlotStyle -> Red]]

enter image description here

Or switch to a global optimizer "NMinimize" which will probably also do the job but is much slower and has the risk of going to regions where your function is ill defined.

POSTED BY: Martijn Froeling

Thank you very much. It solved my problem. So adding the command I'm saying that the function must be evaluated only if the input parameters are numeric. Thanks!

POSTED BY: Lorenzo Fuggiano
Posted 4 years ago

Do you think it would be possible to obtain a more robust system to the variations of the initial parameters? I mean that I would like the fitting results (b, l, d) to be as similar as possible even by varying the initial search values. Thanks and sorry for the further inconvenience.

POSTED BY: Updating Name

Well I think you only have to find the correct order of magnitude of the fit parameters. Which I think you should know beforehand.

If not you could do just a brute force evaluation of different orders of magnitude of the initial parameters and select the best values for your initial guess.

Block[{y = 0, t = 1, b = 0.0015, l = 0.0011, d = 0.00013}, 
  data = Table[{x, 
     fm[x, y, t, b, l, d] + 
      Random[NormalDistribution[0, 0.1]]}, {x, -0.0015, 0.002, 
     0.000015}]];

bi = li = di = {10.^-2, 10.^-3, 10.^-4, 10.^-5}
Block[{y = 0, t = 1},
  all = SortBy[Flatten[Table[
      datai = 
       Table[{x, 
         fm[x, y, t, b, l, d] + 
          Random[NormalDistribution[0, 0.1]]}, {x, -0.0015, 0.002, 
         0.000015}];
      {RootMeanSquare[data[[All, 2]] - datai[[All, 2]]], 
       datai, {b, l, d}},
      {b, bi}, {l, li}, {d, di}], 2], First]
  ];
{lsq, datai, {bi, li, di}} = First[all];

{lsq, {bi, li, di}}
ListLinePlot[all[[All, 2]], PlotRange -> All]
ListLinePlot[{data, datai}, PlotRange -> All]

enter image description here

POSTED BY: Martijn Froeling
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract