1
$\begingroup$
FindFit[data, {a/(1 + i*b*x^2), Element[{a, b}, Reals]}, {a, b}, x]

FindFit::nrgnum: The gradient is not a vector of real numbers at {a,b} = {1.,1.}. FindFit::grad: Evaluation of the gradient of function Experimental`NumericalFunction[{Hold[<<1>>],Block},{0,{{{},1,0,Hold[a],0,0},{{},1,1,Hold[b],0,0}}},{{{1,2,817},{{{Hold[Block[{<<4>>},CompoundExpression[<<3>>]]],Block},Automatic,None,1,Automatic},{}}}},{8,3,{},0},{904,MachinePrecision,{None,{Hold[Block[{x={<<1000>>},Optimization`FindFit`y$45202={<<1000>>}},1/2 > Subtract[<<2>>].Subtract[<<2>>]]],Block}},True,{{Automatic,CleanUpRegisters->False,WarningMessages->False,EvaluateSymbolically->False,RuntimeErrorHandler->($Failed&)},{},Automatic,WVM},FindFit,Automatic,None},{None,None,None}] failed at {1.,1.}.

$\endgroup$
  • 1
    $\begingroup$ The issue seems to be the i. Do you mean the complex unit I? $\endgroup$ – Henrik Schumacher Nov 26 '17 at 8:58
  • $\begingroup$ yes , I am trying to fit this complex function with real parameters . $\endgroup$ – Subhamita Sengupta Nov 26 '17 at 18:29
3
$\begingroup$

To my suprise there seems to be no method for nonlinear model fits to maps into $\mathbb{R}^n$ or to $\mathbb{C}$. Here is a simple implementation of the Gauss-Newton method that works at least in the given example:

MyNonlinearModelFit[pts0_, vals_, model_, pars_, vars_, pguess_] := 
 Module[{p, modelfun, cmodel, cDmodel, error, Derror, gradfit, 
   residual, iter, pts},
  pts = If[VectorQ[pts0],
    Transpose[{pts0}],
    pts0
    ];

  Quiet[Block[{p, pp, x, xx},
    pp = Table[p[[i]], {i, 1, Length[pars]}];
    xx = Table[x[[i]], {i, 1, Length[vars]}];
    modelfun = {x, p} \[Function] 
      Evaluate[model /. Thread[pars -> pp] /. Thread[vars -> xx]];
    With[{code = N[modelfun[xx, pp]] /. Part -> Compile`GetElement},
     cmodel = Compile[{{x, _Real, 1}, {p, _Real, 1}},
       code,
       CompilationTarget -> "C",
       Parallelization -> True,
       RuntimeAttributes -> {Listable}
       ]
     ];

    With[{code = 
       N[D[modelfun[xx, pp], {pp, 1}]] /. 
        Part -> Compile`GetElement},
     cDmodel = Compile[{{x, _Real, 1}, {p, _Real, 1}},
       code,
       CompilationTarget -> "C",
       Parallelization -> True,
       RuntimeAttributes -> {Listable}
       ]
     ];
    ]];

  p = pguess;
  error = Flatten[cmodel[pts, p] - vals];
  Derror = Flatten[cDmodel[pts, p], 1];
  iter = 0;
  gradfit = error.Derror;
  residual = Max[Abs[gradfit]]/Length[error];
  While[residual > 1. 10^-12,
   ++iter;
   p = p - 
     LinearSolve[Derror\[Transpose].Derror, gradfit, 
      Method -> "Cholesky"];
   error = Flatten[cmodel[pts, p] - vals];
   Derror = Flatten[cDmodel[pts, p], 1];
   gradfit = error.Derror;
   residual = Max[Abs[gradfit]]/Length[error]
   ];
  Thread[pars -> p]
  ]

Here is an example for the model given by the OP:

n = 100;
ptrue = {3., 2.};
SeedRandom[1234];
model = ComplexExpand[ReIm[a/(1 + I b x^2)]];
ftrue = x \[Function] Evaluate[model /. Thread[{a, b} -> ptrue]];
pts = RandomReal[{-1, 1}, {n}];
vals = ftrue /@ pts + RandomVariate[
    MultinormalDistribution[{0, 0}, DiagonalMatrix[{1, 1}]], {n}];
pars = {a, b};
vars = {x};
pguess = {1., 0.};
fit = MyNonlinearModelFit[pts, vals, model, pars, vars, pguess]; //AbsoluteTiming
fit
ffitted = x \[Function] Evaluate[model /. fit];
Show[
 ParametricPlot3D[
  Evaluate[{Join[{x}, ftrue[x]], Join[{x}, ffitted[x]]}],
  {x, Min[pts], Max[pts]}, 
  PlotLegends -> {"fitted curve", "true curve"}],
 Graphics3D[Point[Join[Transpose[{pts}], vals, 2]]],
 PlotRange -> All
 ]

{0.150965, Null}

{a -> 3.04508, b -> 1.93957}

enter image description here

Since there are only two parameters involved, this works fine even for quite high amounts of data. (I tested it for up to 10 million data points.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.