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

Questions
, , ,
1

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.

2

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]]

fit.PNG
fit.PNG370×248 8.69 KB

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.

3

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!

4

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.

5

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]

9839fit.PNG
9839fit.PNG381×488 38.6 KB