# Linear fit statistics for a complex data set

I am having some technical and conceptual troubles in using NonlinearModelFit with a complex data set. I have a 2 dimensional array where the two positions are complex numbers, such that:

data={{1 + I, 0.985402 + 1.08528 I}, {1 + 2 I, 2.09444 + 1.00236 I}, {1 + 3 I, 3.29011 + 0.969815 I}, {1 + 4 I,3.98937 + 0.969446 I}, {1 + 5 I, 4.71475 + 0.916196 I}, {1 + 6 I, 5.97777 + 0.994892 I}, {1 + 7 I, 6.35507 + 1.02171 I}, {1 + 8 I, 7.41285 + 0.9147 I}, {1 + 9 I, 8.73952 + 1.04088 I}, {1 + 10 I, 9.36015 + 0.947539 I}}

ansatz= a0 x


The data (indeed mock data in this case) is linear in x so using Fit I get:

Fit[data, {0, x}, x]=  (0.270025 - 0.902848 I) x


The same result is obtained now when using NonlinearModelFit as,

fit=NonlinearModelFit[data, a x, a, x]


However the issue appears when trying to extract some of the fit properties as AIC values, p-values and others. This is the error message,

fit["AIC"]


FittedModel::varnum: The estimated variance -0.836618 is not a positive number. Properties requiring division by the variance or standard error will not be computed.

So in my understanding, the function does not work well with complex numbers where, for instance, the variance seems to not be generically written as the absolute value of the residuals but somehow different. Therefore, my question here is: Is there any theoretical issue that makes the fitting statistics in the complex plane ill defined or it is just a bug in the function? Does it makes any sense the replacement of the variance by its absolute value?

Thanks!

• Questions about statistical models should be asked at CrossValidated. Here's a good reference for the techniques needed: stats.stackexchange.com/questions/66088/…. Then the advice from that reference can be implemented in Mathematica. – JimB Sep 25 '19 at 14:47
• You'll need to state the linear model that you're trying to fit. A problem with your mock data (for a simple "complex number" regression) is that all of the real values of the predictor variable are 1. Is the model $y=a+b x+\epsilon$ where $y$, $a$, $b$, $x$, and $\epsilon$ are all complex variables? If so, the error structure will involve a bivariate normal distribution. – JimB Sep 26 '19 at 4:02
• Oh thanks that link was very useful. I agree is a bivariate distribution but my point is, do i want to my log distribution as |d-y|^2 thus having cross (correlated) terms or i want to decorrelate them as re[d-y]^2+ im[d-y]^2? I guess is my only my choice and only some previous knowldge of the model would give the right answer – Francisco Jimenez Sep 26 '19 at 14:23
• I think the CrossValidated reference addresses the question in your comment (although if you are not a statistician, then it might not appear completely clear). If you can wait a few days, I can write up an example. However, if you could edit your question to include a specific linear model, that would be very helpful. Even if the example just clarifies if the "slope" and "intercept" parameters are real or complex would be very helpful. – JimB Sep 26 '19 at 14:52
• Sorry for the confusion guys but here the vertical bars on a complex number means the quadratic sum of the real & imaginary parts. – Francisco Jimenez Sep 27 '19 at 15:05

The linear model needs to be made explicit as there can be models that might look similar but are really parameterized very differently and would imply possibly very different analysis procedures and interpretation of results.

Using the notation from @whuber 's answer from CrossValidated, if we make the response variable ($$z=x+i y$$), predictor variable ($$w=u+i v$$), slope (b=$$\gamma_1+i \delta_1$$), intercept (a=$$\gamma_0+i \delta_0$$), and error ($$\epsilon_r+i \epsilon_i$$) all complex variables, then a simple linear model might look like the following:

$$z = a + b w + \epsilon$$

$$x + i y = \gamma_0 + i \delta_0 + (\gamma_1 + i \delta_1)*(y + i v) + \epsilon_r + i \epsilon_i$$

where $$\epsilon_r$$ and $$\epsilon_i$$ have a bivariate normal distribution.

Here is a Mathematica function using maximum likelihood to estimate the parameters:

(* X is the design matrix *)
(* z in the response variable *)
(* equalVariances is True if you want the real and imaginary components to have the same variance and False otherwise. *)
(* corr0 is True if you want to assume that the correlation between the real and imaginary error components to be zero.
Otherwise, the correlation will be estimated from the data. *)

complexLinearModelFit[{X_, z_}, equalVariances_, corr0_] := Module[
{model, modelR, modelI, logL, vR, vI, v, lmR, lmI, estimates,
initialValues, ρ0, σR0, σI0, mle, cov, aic, γ, δ, xR, xI, conditions},

(* Regression parameters *)
γ = Table[ToExpression["γ" <> ToString[i]], {i, Dimensions[X][[2]]}];
δ = Table[ToExpression["δ" <> ToString[i]], {i, Dimensions[X][[2]]}];

model = Sum[(γ[[j]] + I δ[[j]]) X[[All, j]], {j, Dimensions[X][[2]]}];
modelR = Re[ComplexExpand[model]] /. Im[h_] -> 0 /. Re[h_] -> h;
modelI = Im[ComplexExpand[model]] /. Im[h_] -> 0 /. Re[h_] -> h;

(* Get lists of variables *)
vR = Variables[modelR];
vI = Variables[modelI];
v = Variables[model];

(* Get regression coefficients associated with each variable *)
xR = Table[Coefficient[modelR, vR[[i]]], {i, Length[vR]}];
xI = Table[Coefficient[modelI, vI[[i]]], {i, Length[vI]}];

(* Perform linear models on real and imaginary components to get initial estimates of coefficients *)
lmR = LinearModelFit[{Transpose[xR], Re[z]}];
lmI = LinearModelFit[{Transpose[xI], Im[z]}];
(* Get the mean of parameter estimates estimate *)
estimates = Transpose[{Join[vR, vI], Join[lmR["BestFitParameters"], lmI["BestFitParameters"]]}];
initialValues = Table[Mean[Select[estimates, #[[1]] == v[[i]] &][[All, 2]]], {i, Length[v]}];
(* Now get initial values for the covariance matrix *)
ρ0 = Correlation[lmR["FitResiduals"], lmI["FitResiduals"]];
σR0 = lmR["EstimatedVariance"]^0.5;
σI0 = lmI["EstimatedVariance"]^0.5;

(* Put together all of the initial estimates,
log likelihood function, and determine conditions on parameters *)
logL = LogLikelihood[
BinormalDistribution[{0, 0}, {σR, σI}, ρ],
Transpose[{Re[z] - modelR, Im[z] - modelI}]];
If[equalVariances,
logL = logL /. σR -> σ /. σI -> σ;
conditions = {σ > 0};
v = Join[v, {σ}];
initialValues = Join[initialValues, {(σR0 + σI0)/2}],

conditions = {σR >= 0, σI >= 0};
v = Join[v, {σR, σI}];
initialValues = Join[initialValues, {σR0, σI0}]
];

If[corr0,
logL = logL /. ρ -> 0,

conditions = Join[conditions, {-1 <= ρ <= 1}];
v = Join[v, {ρ}];
initialValues = Join[initialValues, {ρ0}]
];
initialValues = Transpose[{v, initialValues}];

(* Find maximum likelihood estimates *)
mle = FindMaximum[{logL, conditions}, initialValues];

(* Estimates of standard errors *)
(cov = -Inverse[(D[logL, {v, 2}]) /. mle[[2]]]);

(* AIC *)
aic = -2 mle[[1]] + 2*Length[v];

(* Return results from FindMaximum, covariance matric estimate, and AIC *)
{mle, cov, aic}
]


Now for some examples. First consider the data posted:

data = {{1 + I, 0.985402 + 1.08528 I}, {1 + 2 I, 2.09444 + 1.00236 I},
{1 + 3 I, 3.29011 + 0.969815 I}, {1 + 4 I, 3.98937 + 0.969446 I},
{1 + 5 I, 4.71475 + 0.916196 I}, {1 + 6 I, 5.97777 + 0.994892 I},
{1 + 7 I, 6.35507 + 1.02171 I}, {1 + 8 I, 7.41285 + 0.9147 I},
{1 + 9 I, 8.73952 + 1.04088 I}, {1 + 10 I, 9.36015 + 0.947539 I}};

response = data[[All, 2]];
X = Transpose[{ConstantArray[1, Length[data]], data[[All, 1]]}];

lm = complexLinearModelFit[{X, response}, False, False];
(* Parameter estimates *)
clm[[1, 2]]
(* {γ1 -> 0.269328, γ2 -> -0.00612192, δ1 -> 1.93427, δ2 -> -0.914316,
σR -> 0.193555, σI -> 0.0483416, ρ -> 0.0186871} *)
(* Covariance matrix *)
clm[[2]] // TableForm


(* AIC *)
clm[[3]]
(* -22.6791 *)


Now for the example given on CrossValidated:

(* Predictor variable *)
w = {0 - 5 I, -3 - 4 I, -2 - 4 I, -1 - 4 I, 0 - 4 I, 1 - 4 I, 2 - 4 I, 3 - 4 I, -4 - 3 I, -3 - 3 I, -2 - 3 I, -1 - 3 I, 0 - 3 I, 1 - 3 I, 2 - 3 I, 3 - 3 I, 4 - 3 I, -4 - 2 I, -3 - 2 I, -2 - 2 I, -1 - 2 I, 0 - 2 I, 1 - 2 I, 2 - 2 I, 3 - 2 I, 4 - 2 I, -4 - 1 I, -3 - 1 I, -2 - 1 I, -1 - 1 I, 0 - 1 I, 1 - 1 I, 2 - 1 I, 3 - 1 I, 4 - 1 I, -5 + 0 I, -4 + 0 I, -3 + 0 I, -2 + 0 I, -1 + 0 I, 0 + 0 I, 1 + 0 I, 2 + 0 I, 3 + 0 I, 4 + 0 I, 5 + 0 I, -4 + 1 I, -3 + 1 I, -2 + 1 I, -1 + 1 I, 0 + 1 I, 1 + 1 I, 2 + 1 I, 3 + 1 I, 4 + 1 I, -4 + 2 I, -3 + 2 I, -2 + 2 I, -1 + 2 I, 0 + 2 I, 1 + 2 I, 2 + 2 I, 3 + 2 I, 4 + 2 I, -4 + 3 I, -3 + 3 I, -2 + 3 I, -1 + 3 I, 0 + 3 I, 1 + 3 I, 2 + 3 I, 3 + 3 I, 4 + 3 I, -3 + 4 I, -2 + 4 I, -1 + 4 I, 0 + 4 I, 1 + 4 I, 2 + 4 I, 3 + 4 I, 0 + 5 I};
(* Add in a "1" for the intercept *)
w1 = Transpose[{ConstantArray[1 + 0 I, Length[w]], w}];

z = {-15.83651 + 7.23001 I, -13.45474 + 4.70158 I, -13.63353 + 4.84748 I, -14.79109 + 4.33689 I, -13.63202 + 9.75805 I, -16.42506 + 9.54179 I, -14.54613 + 12.53215 I, -13.55975 + 14.91680 I, -12.64551 + 2.56503 I, -13.55825 + 4.44933 I, -11.28259 + 5.81240 I, -14.14497 + 7.18378 I, -13.45621 + 9.51873 I, -16.21694 + 8.62619 I, -14.95755 + 13.24094 I, -17.74017 + 10.32501 I, -17.23451 + 13.75955 I, -14.31768 + 1.82437 I, -13.68003 + 3.50632 I, -14.72750 + 5.13178 I, -15.00054 + 6.13389 I, -19.85013 + 6.36008 I, -19.79806 + 6.70061 I, -14.87031 + 11.41705 I, -21.51244 + 9.99690 I, -18.78360 + 14.47913 I, -15.19441 + 0.49289 I, -17.26867 + 3.65427 I, -16.34927 + 3.75119 I, -18.58678 + 2.38690 I, -20.11586 + 2.69634 I, -22.05726 + 6.01176 I, -22.94071 + 7.75243 I, -28.01594 + 3.21750 I, -24.60006 + 8.46907 I, -16.78006 - 2.66809 I, -18.23789 - 1.90286 I, -20.28243 + 0.47875 I, -18.37027 + 2.46888 I, -21.29372 + 3.40504 I, -19.80125 + 5.76661 I, -21.28269 + 5.57369 I, -22.05546 + 7.37060 I, -18.92492 + 10.18391 I, -18.13950 + 12.51550 I, -22.34471 + 10.37145 I, -15.05198 + 2.45401 I, -19.34279 - 0.23179 I, -17.37708 + 1.29222 I, -21.34378 - 0.00729 I, -20.84346 + 4.99178 I, -18.01642 + 10.78440 I, -23.08955 + 9.22452 I, -23.21163 + 7.69873 I, -26.54236 + 8.53687 I, -16.19653 - 0.36781 I, -23.49027 - 2.47554 I, -21.39397 - 0.05865 I, -20.02732 + 4.10250 I, -18.14814 + 7.36346 I, -23.70820 + 5.27508 I, -25.31022 + 4.32939 I, -24.04835 + 7.83235 I, -26.43708 + 6.19259 I, -21.58159 - 0.96734 I, -21.15339 - 1.06770 I, -21.88608 - 1.66252 I, -22.26280 + 4.00421 I, -22.37417 + 4.71425 I, -27.54631 + 4.83841 I, -24.39734 + 6.47424 I, -30.37850 + 4.07676 I, -30.30331 + 5.41201 I, -28.99194 - 8.45105 I, -24.05801 + 0.35091 I, -24.43580 - 0.69305 I, -29.71399 - 2.71735 I, -26.30489 + 4.93457 I, -27.16450 + 2.63608 I, -23.40265 + 8.76427 I, -29.56214 - 2.69087 I};

(* Estimation assuming ρ=0 (which is what the CrossValidated example assumes *)
complexLinearModelFit[{w1, z}, True, True][[1, 2]]
(* {γ1 -> -20.0172, γ2 -> -0.830797, δ1 -> 5.00968, δ2 -> 1.37827, σ -> 2.20038} *)

(* Now allow the estimation of ρ *)
complexLinearModelFit[{w1, z}, True, False][[1, 2]]
(* {γ1 -> -20.0172, γ2 -> -0.763237, δ1 -> 5.00968, δ2 -> 1.30859, σ -> 2.21424, ρ -> 0.810525} *)


The true values are γ1 -> -20, γ2 -> -0.75, δ1 -> 5, δ2 -> 1.299038, σ -> 2, and ρ -> 0.8 which are very close to the estimates when $$\rho$$ is allowed to be estimated.

The code should work for any number of complex predictors. But, again, that assumes that the models being fit with the code are the models that you need.

The default assumption that NonlinearModelFit makes to compute these statistics, is that the residuals are normally distributed. You can access the residuals easily with fit["FitResiduals"] and obviously these are also complex values now, so what does it mean exactly to say that these residuals are normally distributed? Normal distributions generate real numbers, so obviously the default assumptions fail here.

I'm not saying nothing can be done, but you'll need to define your own model to analyze the complex residuals with a complex-valued probability distribution and estimate its parameters. Once you've fitted such a distribution, you can just calculate the log-likelihood of the errors using this distribution and calculated the AIC using the formula:

AIC = 2 k - 2 LogLikelihood[dist, residuals] (* with k the number of free fitting parameters *)


It's probably not difficult, but I've never fitted complex data myself, so I don't know exactly how people generally extend the normal distribution into the complex plane.

edit

As a simple example, you could estimate the normal distribution of the absolute values of the error like this:

fit = NonlinearModelFit[data, a x, a, x]
residuals = Abs[fit["FitResiduals"]];
logLike = LogLikelihood[
EstimatedDistribution[residuals, NormalDistribution[mu, sigma]],
residuals
]
aic = 2 * 2 - 2*logLike


edit2

As JimB remarked, estimating the absolute values of the errors with a NormalDistribution is a bit senseless, so let me try to give a few better possibilities.

For one, we could imagine squaring the absolute values of the residuals and assume that the residuals follow a MultiNormalDistribution in the complex plane. A 2D multi normal distribution is a BinormalDistribution, so we can make the following estimate of the log-likelihood:

residuals = ReIm[fit["FitResiduals"]];
errorDist = EstimatedDistribution[
residuals,
BinormalDistribution[{sigma1, sigma2}, rho]
]
logLike = LogLikelihood[errorDist, residuals]


However, we need to estimate 3 parameters to do this, which might be a bit overkill. Instead, we can simplify the error model and say that both the real and imaginary parts of the errors have to be generated by NormalDistribution[0, sigma] (i.e., the real and imaginary parts have the same standard deviation):

errorDist = EstimatedDistribution[
residuals,
ProductDistribution[{NormalDistribution[0, sigma], 2}]
]


As you can see, there are different choices can make for analysing the errors. If you want to use error models other than those from the Gaussian family, you probably need to go back to square one and also rethink how to fit the data in the first place, since it's quite likely that least square estimation wouldn't be appropriate for such error models.

• Hi Sjoerd, thanks for your answer. True on the normal distribution issue. Indeed , taking abs[residuals] you can have again a well defined problem however i was unsure about that due to the poor bibliography i found. – Francisco Jimenez Sep 26 '19 at 14:15
• @FranciscoJimenez I'll update the answer with a brief example of how to implement your suggestion to take the absolute values. – Sjoerd Smit Sep 26 '19 at 15:46
• Would you add in what linear model is being estimated? I don't see that in your current answer or in the OP. Also, I might be misinterpreting what the OP means by "absolute values" but is the "norm" what is being requested? – JimB Sep 26 '19 at 16:04
• Isn't "absolute value" a common synonym for the norm when it comes to complex numbers? What else could it be? – Sjoerd Smit Sep 26 '19 at 21:59
• Sorry, I didn't phrase that very well. I don't understand why one would want to look at the absolute value of the residuals and assume that those positive numbers would have a normal distribution. – JimB Sep 26 '19 at 23:33