# Fitting fractional complex data with NonlinearModelfit

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.

-
Here's why NonlinearModelFit on complex data sometimes works and sometimes doesn't: mathematica.stackexchange.com/a/54876/484 –  Rahul Jul 15 at 16:27

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___ ] /; InternalDependsOnQ[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} *)
-
You could post this here or maybe here for more visibility. –  Szabolcs Mar 2 at 19:15
@Szabolcs thanks. I will post it in István's thread, which somehow I had missed until now. –  Oleksandr R. Mar 2 at 19:18
(+1) One additional limitation: in the most practical cases the conservative assumption is that the number of observations is equal to the number of complex data points, not twice the number. In such cases current implementation of ComplexFit significantly underestimates the confidence intervals and standard errors. –  Alexey Popkov Jul 15 at 7:09
@AlexeyPopkov very good point, and one that I had overlooked; thank you. I would obviously like to address this, but (a) I'm not sure how to without somehow hijacking functions that work on FittedModel objects and (b) the situation is arguably problem-dependent. I do agree with you though that a conservative estimate would be better. –  Oleksandr R. Jul 15 at 8:37
@AlexeyPopkov yes, that's possible at present (although these two sublists would simply be joined into one) but one cannot guess the weights automatically, so I simply didn't do anything with them. There is nothing stopping someone from using the Weights and VarianceEstimatorFunction options with ComplexFit provided they know how the transformation is done. For this you can specify "Hold" -> True and the transformed FindFit/NonlinearModelFit call is returned held rather than being evaluated. –  Oleksandr R. Jul 15 at 9:00

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.

-
I am just confused that NonlinearModelFit works in the first case with (a + b*x^2*I) , but not if I want to fit the inverse- which should be the same in principle. I know that there are different ways to solve this fitting problem - but i would like to use nonlinearmodelfit since it seems to be more effective for more difficult problems (which I have to deal with) and gives you some errors for the parameters –  Justus Mar 10 '12 at 18:44

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