Fitting different data sets with different models using FindFit

Questions
, ,
1

Referring to a fitting made previously thanks to the help of a member:
Fit noisy data as a Least squares problem?

I would like to ask for help to fit 2 noisy data sets to 2 models, referring to the same equation with different conditions.
I have an equation T1:

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

I create noisy data for y=0 and x=0:

Block[{y = 0, t = 1, b = 0.001, l = 0.001, d = 0.0001, 
   m = -Tan[45*Pi/180], Intens = 25000}, 
  datax = Table[{x, 
     T1[x, y, t, b, l, d, m, Intens] + 
      Random[NormalDistribution[0, 0.1]]}, {x, -0.002, 0.002, 
     0.000015}]];
Block[{x = 0, t = 1, b = 0.001, l = 0.001, d = 0.0001, 
   m = -Tan[45*Pi/180], Intens = 25000}, 
  datay = Table[{y, 
     T1[x, y, t, b, l, d, m, Intens] + 
      Random[NormalDistribution[0, 0.1]]}, {y, -0.002, 0.002, 
     0.000015}]];

and after, all I want to do, is to fit the two noisy data sets to the equation above for y=0 and x=0, in order to find the best fitting parameters {b,l,d,m,Intens}.

I thought about using Join as follow (t=1):

alldata = Join[datax, datay];
  fit=FindFit[alldata,T1[x,y,1,b,l,d,m,Intens],{b,l,d,m,Intens},{x,y},Method->"LevenbergMarquardt"]

but nothing. Could anyone help me? Thanks.

2

If I reformulate with Integrate instead of NIntegrate it does not give an error:

FindFit[alldata,
 0.3183098861837907` Intens Integrate[
   Sqrt[1 + m^2]*
    Erfc[Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\[Alpha] - 
              d))^2]/Sqrt[
        Dif*4*1]]/(Sqrt[(x - \[Alpha])^2 + (y - \[Beta])^2 + (z - (m*\
\[Alpha] - d))^2]), {\[Beta], -b/2, b/2}, {\[Alpha], linf\[Alpha], 
    l*(1/(1 + m^2))}],
 {b, l, d, m, Intens}, {x, y}, Method -> "LevenbergMarquardt"]

but it takes such a long time that I did not wait. The numerical evaluation of the double integral for so many different values of the parameters may be very onerous.

3

That is interesting.

I had a look at your data ( datax and datay ) and tried this:

ff[x_, a_, b_, c_] := a  Exp[-b (x - c)^2]

and

sola =  FindFit[datax, ff[x, a, b, c], {{a, 5.85}, {b, 5.1 10^6}, {c, 0.00015}}, x]

(* {a -> 5.54423, b -> 4.52534*10^6, c -> 0.000158701} *)

Then

Plot[Evaluate[ff[x, a, b, c] /. sola], {x, -0.002, 0.002},
 Epilog -> Point /@ datax, PlotRange -> {0, 7}]

Similarly

solb =  FindFit[datay,   ff[x, a, b, c], {{a, 5.05}, {b, 2.56 10^6}, {c, -0.00001}}, x]
(* {a -> 5.54337, b -> 2.59886*10^6, c -> 1.4393*10^-7} *)
Plot[Evaluate[ff[x, a, b, c] /. solb], {x, -0.002, 0.002},
 Epilog -> Point /@ datay, PlotRange -> {0, 7}]

These are the results from FindFit, and I believe it. Anyhow, one could get the impression (optically) that there is somethig better for datax, but note, this is an impression. Play around with a, b and c:

Manipulate[
 Plot[ff[x, a, b, c], {x, -0.002, 0.002}, PlotRange -> {0, 7},
  Epilog -> Point /@ datax
  ],
 {a, 5, 7}, {b, .5 10^6., 10^7}, {c, -.002, .002}]
4

I tried to perform this calculation, but I don’t get any value from the parameters. I also get an error:

FindFit::fitc: Number of coordinates (1) is not equal to the number of variables (2).
5

Thank you for the reply. Anyway, I don’t understand well. ff in this case what does it refer to?

6

Summary: I want to fit the noisy data to the model I found them with (T1). The parameters are the same for both the fit to be made {b,l,d,m,Intens}. For datax the free variable is x. For datay the free variable is y. I don’t want the problem to be due to this.

7

Hello Lorenzo,

ff is nothing but an empirical function to give an appoximation for your data.

At least it shows that your complicated integral seems to have something to do with the error-function.

And I think you can somewhat simplify your integral.

If I didn’t make a mistake your integrand is (could you please check this?)

jj = (Sqrt[1 + m^2]
    Erfc[Sqrt[(x - \[Alpha])^2 + (d + z - m \[Alpha])^2 + (y - \[Beta])^2]/(2 Sqrt[Dif t])])/
  Sqrt[(x - \[Alpha])^2 + (d + z - m \[Alpha])^2 + (y - \[Beta])^2];

Now introduce some abbreviations:

rule = {k1 -> d^2 + x^2 + 2 d z + z^2 - (d m + x + m z)^2/(1 + m^2), 
   k2 -> 1/(4 Dif t), \[Alpha]1 -> \[Alpha] - (d m + x + m z)/(1 + m^2), \[Beta]1 -> \[Beta] - y};

then you can write your integrand as

jj3 = Sqrt[1 + m^2]
    Erfc[Sqrt[k1 k2 + k2 (1 + m^2) \[Alpha]1^2 + k2 \[Beta]1^2 ]]/
   Sqrt[k1 + (1 + m^2) \[Alpha]1^2 + \[Beta]1^2];
jj - (jj3 /. rule) // Simplify

( I did not succeed to bring the arguments of the Erfc into the same form, but they are equal)

I define a coordinate-transformation for the variables beta1 = fb and alpha1 = fa

fb = Sqrt[(u + v)^2 - k1] Sin[v];
fa = Sqrt[((u + v)^2 - k1)/(1 + m^2)] Cos[v];
rulec = {\[Beta]1 -> fb, \[Alpha]1 -> fa};

and get as the new integrand

jj4 = (jj3 /. rulec) Det[( { {D[fb, u], D[fb, v]}, {D[fa, u], D[fa, v]} } )] // PowerExpand // FullSimplify
jj5 = jj4 /. Sqrt[k2 (u + v)^2] :> (u + v) Sqrt[k2] /. 1/Sqrt[(u + v)^2] :> 1/(u + v)

namely

-Erfc[Sqrt[k2] (u + v)]

and this can be integrated explicitly

Integrate[Integrate[-Erfc[Sqrt[k2] (u + v)], u], v]

or

Integrate[Integrate[-Erfc[Sqrt[k2] (u + v)], {u, u1, u2}], {v, v1, v2}]

If you do it you will see that there are a couple of error - functions to be seen in the result.

What I didn’t do is to look for the boundaries for u and v to be compatible with the region spanned by alpha1 and beta1 (if this is possible at all).

But I think you can do the backtransformation (with several solutions)

sol = Solve[{fb == bb, fa == aa}, {u, v}];

and look for example at

jj5 /. sol[[6]] /. rule

or

jj5 /. sol[[4]] /. rule
8

You are fitting multivariate data (function dependent on x and y) but your data only has one coordinate logged (either x or y). so either add the respective x and y values to your data or just generate the full data range for x and y.

Both are equally valid in this case.

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

t = 1;
b = 0.001;
l = 0.001;
d = 0.0001;
m = -Tan[45*Pi/180];
Intens = 25000;

step = 0.00006;

datax = Table[{x, 0, 
    T1[x, 0, t, b, l, d, m, Intens] + 
     Random[NormalDistribution[0, 0.1]]}, {x, -0.002, 0.002, step}];
datay = Table[{0, y, 
    T1[0, y, t, b, l, d, m, Intens] + 
     Random[NormalDistribution[0, 0.1]]}, {y, -0.002, 0.002, step}];
alldata = Join[datax, datay];

dataxy = Flatten[
   Table[{x, y, 
     T1[x, y, t, b, l, d, m, Intens] + 
      Random[NormalDistribution[0, 0.1]]}, {x, -0.002, 0.002, 
     4 step}, {y, -0.002, 0.002, 4 step}], 1];

Show[
 Plot[{T1[x, 0, t, b, l, d, m, Intens], 
   T1[0, x, t, b, l, d, m, Intens]}, {x, -0.002, 0.002}, 
  PlotStyle -> Red],
 ListPlot[{datax[[All, {1, 3}]], datay[[All, {2, 3}]]}, 
  PlotStyle -> {Black, Gray}]
 ]

Show[
 Plot3D[
  T1[x, y, t, b, l, d, m, Intens], {x, -0.002, 0.002}, {y, -0.002, 
   0.002},
  PlotRange -> Full, Mesh -> False, PlotStyle -> Gray, 
  Lighting -> "Neutral"],
 ListPointPlot3D[dataxy, PlotRange -> Full, PlotStyle -> Black],
 Graphics3D[{Thick, Red, Line[datax], Line[datay]}]
 ]

(*only lines*)
fitL = FindFit[alldata, 
  T1[x, y, 1, bf, lf, d, m, Intens], {{bf, 0.005}, {lf, 0.005}(*,d,m,
   Intens*)}, {x, y}, Method -> "LevenbergMarquardt"]

(*full surface*)
fitS = FindFit[dataxy, 
  T1[x, y, 1, bf, lf, d, m, Intens], {{bf, 0.005}, {lf, 0.005}(*d,m,
   Intens*)}, {x, y}, Method -> "LevenbergMarquardt"]
9

Now I get it. The problem is that I must necessarily use my original function. Thanks anyway.

10

You found the error. As usual, thank you!

11

I have another problem with the aforementioned fit.
I would like to perform the fit considering also as an independent variable “t”.
The latter is described by two similar equations, the first (T1) one for 0 <= t <= tau and the second (T2) for t > tau.

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

In producing noisy data I already get the first problems:

Block[{x = 0, y = 0, b = 0.001, l = 0.001, d = 0.0001, 
   m = -Tan[45*Pi/180], Intens = 25000}, 
  datat1 = Table[{0, 0, t, 
     T1[x, y, t, b, l, d, m, Intens] + 
      Random[NormalDistribution[0, 0.1]]}, {t, 0, 1, 0.01}]];

NIntegrate::inumr: The integrand (Sqrt[2] Erfc[ComplexInfinity])/Sqrt[\[Alpha]^2+(0.0001 +\[Alpha])^2+\[Beta]^2] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,0.0005},{0.,0.0005}}.
Power::infy: Infinite expression 1/0. encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.

Then I thought about doing the same thing for T2 and then using Join for the 4 data sets. But I think the problem starts already in the creation of noisy data.

I also tried with:

If[0 <= t <= 1, T1[0, 0, t, b, l, d, m, Intens], 
 T2[0, 0, t, b, l, d, 1, m, Intens]]

but nothing happens.

I wish I could fit the 4 datasets (for x, y and t) together with the 2 equations, always if this is possible.