Fitting fractional complex data with NonlinearModelfit

Question

I've difficulites with the NonlinearModelFit function.

In principle Mathematica should be able to deal with complex data.

E.g. if I define the following table

set = Table[{i, (3*i^2*I + 1*I)}, {i, 0, 10}]

and perform NonlinearModelFit, Mathematica has no problem finding the right values

NonlinearModelFit[set, (a + b*x^2*I), {a, b}, x]
(* {a -> 0. + 1. I, b -> 3. + 0. I} *)

But I need to fit some fractional functions, so I tried to fit just the inverse of the above set of data (I don't change any syntax, besides putting 1/ in front)

set = Table[{i, 1/(3*i^2*I + 1*I)}, {i, 0, 10}]
NonlinearModelFit[set, 1/(a + b*x^2*I), {a, b}, x]

Then Mathematica comes up with an error:

FindFit::nrlnum: The function value {1.-1. I,0.5-4.5 I,<<7>>,0.000152393-244.012 I,
<<1>>} is not a list of real numbers with dimensions {11} at {a,b} = {1.,1.}. >>

So now I'm confused- there is no syntax or misspelling problem – is Mathematica just not able to perform such a fractional fit? Or are there any mathematical restrictions, which I have not considered? I'm using Mathematica 7.

Answer 1

This is, at least in principle, a duplicate of jVincent's answer and the one I gave here. The general approach has been suggested by various people over the years, although I first encountered it here courtesy of Daniel Lichtblau the first time I needed to fit several datasets simultaneously.

I've been meaning to post this package for a while, ideally after having generalized it further, but given that generalization is somewhat complicated, and yet this remains a common question, I think on balance it's probably worthwhile to post the code as it stands. Despite certain limitations (listed below), it seems good enough for most applications that require fitting complex data.

BeginPackage["TransformedFit`"];

ClearAll[TransformedParameter];
SetAttributes[TransformedParameter, HoldRest];

Unprotect[TransformedFit]; ClearAll[TransformedFit];

Unprotect[ComplexFit]; ClearAll[ComplexFit];

Begin["`Private`"];

(* Transform numeric quantities rather than renaming them *)
TransformedParameter[t_, num_?NumericQ] := t[num];

(* Generate unique symbols for each transformed parameter -- this avoids difficulties
caused by overzealous common subexpression elimination when the models are compiled *)
TransformedParameter[t_, p_] := TransformedParameter[t, p] =
   With[{sym = Unique["TransformedParameter$", Temporary]},
    (* Unset memo when cleared to facilitate garbage collection *)
    sym /: clear : (Clear | ClearAll | Remove)[___, sym, ___] :=
     clear /; (TransformedParameter[t, p] =.; True);

    (* Display as the parameter by itself if the transformation is Identity *)
    sym /: MakeBoxes[sym, form_] :=
     With[{boxes = MakeBoxes[p, form]},
       InterpretationBox[boxes, sym]
      ] /; t === Identity;

    sym /: MakeBoxes[sym, form_] :=
     With[{boxes = MakeBoxes[t[p], form]},
      InterpretationBox[boxes, sym]
     ];

    sym
   ];

ClearAll[$FitFunctions];
$FitFunctions = If[$VersionNumber >= 7,
   {FindFit, NonlinearModelFit},
   {FindFit}
  ];

Options[TransformedFit] = {
   "FitFunction" -> First[$FitFunctions],
   "Transformation" -> Identity,
   "ParameterTransformation" -> Identity,
   "Hold" -> False
  };

TransformedFit::cons = 
  "The constraint(s), `1`, should be given in terms of the transformed parameters only.";

TransformedFit[
    data_, {model_, cons_}, pars_, vars_,
    opts___
   ] /; Internal`DependsOnQ[cons, pars] :=
  Message[TransformedFit::cons, cons];

(* Deal with data given as ordinate values only *)
TransformedFit[
   data_?VectorQ, model_, pars_, vars_,
   opts : OptionsPattern[{TransformedFit, Sequence @@ $FitFunctions}]
  ] :=
  TransformedFit[
   Transpose[{Range@Length[data], data}], model, pars, vars,
   opts
  ];

TransformedFit[
   data_?MatrixQ, {model_, cons_} | {model_} | model_, pars_, vars_,
   opts : OptionsPattern[{TransformedFit, Sequence @@ $FitFunctions}]
  ] :=
  With[{
    fitFunction = If[MemberQ[$FitFunctions, #], #, First[$FitFunctions]] & @ OptionValue["FitFunction"],
    transformations = List@OptionValue["Transformation"]~Flatten~1
   },
   Block[{
     transformedData,
     transformedParameters, unusedParameterMask, parameterRules,
     transformedModel,
     i
    },
    (* TRANSFORM DATA *)
    With[{
      abscissae = Take[data, All, {1, -2}],
      ordinates = Take[data, All, {-1}]
     },
     transformedData = {
         ConstantArray[Range@Length[transformations], {Length[abscissae], 1}]~Transpose~{2, 3, 1},
         ConstantArray[abscissae, Length[transformations]],
         Through@transformations[ordinates]
        }~Flatten~{{2, 3}, {1, 4}};
     ];

    (* TRANSFORM PARAMETERS *)
    transformedParameters = Outer[
       TransformedParameter,
       transformations, {pars}~Flatten~1
      ] // Transpose;
    With[{
      (* Original and transformed parameters without initial guesses *)
      originalParameterNames = Replace[pars, {p_, __?NumericQ} :> p, {1}],
      transformedParameterNames = Replace[transformedParameters, {p_, __?NumericQ} :> p, {2}]
     },
     With[{
        (* Representations of the original parameters in terms of their transformations *)
        parameterRepresentations = OptionValue["ParameterTransformation"] @@@ transformedParameterNames
       },
       unusedParameterMask = MapThread[
         Composition[Thread, Unevaluated, FreeQ],
         {parameterRepresentations, transformedParameterNames}
        ];
       Clear @@  Flatten@Pick[transformedParameterNames, unusedParameterMask];
       parameterRules = Thread[originalParameterNames -> parameterRepresentations];
      ];
    ];

    (* TRANSFORM MODEL *)
    With[{
      reparameterizedModel = model /. parameterRules,
      KroneckerDelta = If[Equal[##], 1, 0] & (* compilable *)
     },
     transformedModel = Inner[
        #1[reparameterizedModel] KroneckerDelta[i, #2] &,
        transformations, Range@Length[transformations]
       ];
    ];

    (* PERFORM FIT *)
    If[TrueQ@OptionValue["Hold"], Composition[Hold, fitFunction], fitFunction][
     transformedData,
     {transformedModel, cons},
     Pick[transformedParameters, unusedParameterMask, False]~Flatten~1,
     {i, vars}~Flatten~1,
     FilterRules[{opts}, Options[fitFunction]]
    ]
   ]
  ];

Protect[TransformedFit];

ClearAll[coordinateSystemRules];
coordinateSystemRules["Cartesian"] = Sequence[
   "Transformation" -> {Re, Im},
   "ParameterTransformation" -> (#1 + I #2 &)
  ];
coordinateSystemRules["Polar"] = Sequence[
   "Transformation" -> {Abs, Arg},
   "ParameterTransformation" -> (#1 Exp[I #2] &)
  ];
coordinateSystemRules["Real"] = Sequence[
   "Transformation" -> {Re, Im},
   "ParameterTransformation" -> (#1 &)
  ];
coordinateSystemRules["Imaginary"] = Sequence[
   "Transformation" -> {Re, Im},
   "ParameterTransformation" -> (I #2 &)
  ];
(* Default to Cartesian coordinates *)
coordinateSystemRules[_] = coordinateSystemRules["Cartesian"];

Options[ComplexFit] = {
   "CoordinateSystem" -> Automatic
  };

ComplexFit[
   data_, model_, pars_, vars_,
   opts : OptionsPattern[{ComplexFit, TransformedFit, Sequence @@ $FitFunctions}]
  ] :=
  TransformedFit[
   data, model, pars, vars,
   coordinateSystemRules@OptionValue["CoordinateSystem"],
   FilterRules[{opts}, Except["CoordinateSystem" | "Transformation" | "ParameterTransformation"]]
  ];

Protect[ComplexFit];

End[];

EndPackage[];

Package (.m) and notebook files are also available.

The primary limitations are:

• the Weights option is not (directly) supported, because it isn't clear to me how one should transform the weights in general when splitting a complex-valued function into a multivalued real mapping
• the returned FittedModel objects still contain a reference to an index, i, that labels the coordinates (e.g. real line/imaginary line, modulus/argument, or whatever applies to any other coordinate system one may choose), because the structure of these objects is undocumented and I didn't yet figure out how to split them up
• the transformation is done quite rigidly and is not currently versatile enough to cater for all foreseeable situations

Anyway, let's give it a try:

ComplexFit[
  Table[{i, I + 3*i^2 I}, {i, 0, 10}],
  a + b*x^2 I, {a, b}, x,
  "FitFunction" -> NonlinearModelFit
 ]["ParameterConfidenceIntervalTable"]

Or, in polar coordinates:

ComplexFit[
  Table[{i, I + 3*i^2 I}, {i, 0, 10}],
  a + b*x^2 I, {a, b}, x,
  "FitFunction" -> NonlinearModelFit,
  "CoordinateSystem" -> "Polar"
 ]["ParameterConfidenceIntervalTable"]

And, just for fun, here's an example using FindFit instead of NonlinearModelFit, and where the parameters and the initial guesses of their values are explicitly complex:

ComplexFit[
 Table[{x, (17.381 + 53.249 I) x^(1.897 + 0.632 I)}, {x, -10, 10}],
 (a x^b), {{a, 20 + 50 I}, {b, 2 + 0.5 I}}, x
]
(* -> {Re[a] -> 17.381, Im[a] -> 53.249, Re[b] -> 1.897, Im[b] -> 0.632} *)

This is also useful for fitting real-valued data where the model may become erroneously complex-valued for certain values of the parameters. For example, from the other question:

ComplexFit[
 {{0.0, 100.0}, {0.02, 81.87}, {0.04, 67.03},
  {0.06, 54.88}, {0.08, 44.93}, {0.10, 36.76}},
 a b^t, {a, b}, t,
 "CoordinateSystem" -> "Real"
]
(* -> {Re[a] -> 100.004, Re[b] -> 0.0000452493} *)

Answer 2

NonlinearModelFit seeks to reduce the sum of the squares of the residuals. In this case, because you have used complex numbers, these residuals are complex values. The residuals need to be real for NonlinearModelFit to work.

I haven´t encountered a situation with complex residuals before. To be quite honest, I am bit confused by the idea. The complex numbers are not an order-able field, but I guess this doesn't really matter...

At any rate, what you probably want to do is to use FindFit instead of NonlinearModelFit. It has an option that allows you to set how the residuals are compared to each other:

FindFit[set, 1/(a + b*x^2*I), {a, b}, x, 
 NormFunction -> (Norm[Abs[#], 2] &)]

If this isn't what you want, you can probably build a fitting function really easily by using NMinimize. You just need to explicitly define how you want to minimize the residuals.

Answer 3

I can't say I don't know why NonlinearModelFit doesn't just work smoothly. As Searke mentions, having complex residuals is quite exotic, however typically you would simply use the norm of the difference. This however does not work still. So what I typically do is to take my complex data {{x[i],yR[i]+I yI[i]}...} and change it into {{1,x[1],yR[1]},{1,x[1],yR[1]}...,{2,x[1],yI[1]},{2,x[1],yI[1]}}. This way the first coordinate indicates whether you have the real or imaginary component, so you can just change your model to depend on one more parameter that switches between these. The first part can be put into a function as:

domainExtend[v_,domains_]:={
 Join@@Table[ConstantArray[i,Length@v\[Transpose][[1]]],{i,Length@domains}],
 Join@@Table[v\[Transpose][[1]],{i,Length@domains}],
 Join@@Table[domains[[i]][v\[Transpose][[2]]],{i,Length@domains}]
}//Transpose

To select the domain the following function is used, since it avoids the problems with Part printing errors when it's not yet evaluated and also rounding the domain variable:

domainSwitch[p_?NumericQ,domain_] := domain[[Round@p]]

And then your fit is just:

set = Table[{i, 1/(3*i^2*I + 1*I)}, {i, 0, 10}];
res=NonlinearModelFit[domainExtend[set, {Re, Im}], 
     domainSwitch[domain, {Re, Im}]@(1/(a + b*x^2*I)), {a, b}, {domain, x}]

So this now allows you to get all the typical functions of the fittedModel from NonlinearModelFit like "ANOVATable". The downside is that the data is suddenly on a different format, and so is the returned function, so you need to take this into account when plotting the results:

 Show[
 Plot[{res["Function"][1, x], res["Function"][2, x]}, {x, 0, 10}],
 ListPlot[{
           Cases[res["Data"], {1, _, _}][[1 ;;, 2 ;; 3]],
           Cases[res["Data"], {2, _, _}][[1 ;;, 2 ;; 3]]
         }] 
 ]