Linear fit statistics for a complex data set

Question

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!

Answer 1

Edit: I added what I hope is a real answer (bad pun intended).

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.

Answer 2

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.