2 replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/
source | link

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

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 questionother question:

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.

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:

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.

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:

1
source | link

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

Fit to fractional complex data in Cartesian coordinates

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

Fit to fractional complex data in polar coordinates

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} *)