# Analysis with complex data, anything different?

Say for example you are doing a linear model, but the data $y$ is complex.

$y = x \beta + \epsilon$

My data set is complex, as in all the numbers in $y$ are of the form $(a + bi)$. Is there anything procedurally different when working with such data?

I ask because, you will end up getting complex covariance matrices, and test statistics which are complex valued..

Do you need to use a conjugate transposes instead of transposes when doing least squares? is a complex valued covariance meaningful?

• Consider complex number as two separate variables, and that way remove i from all your equations. Otherwise it will be a nightmare... – sashkello Jul 31 '13 at 7:09
• Any information on $x$ or $\beta$? – Stijn Jul 31 '13 at 10:53
• @Sashkello What "nightmare"? The dimensions are halved when you use complex numbers, so arguably that's a simplification. Moreover, you have turned a bivariate DV into a univariate DV, which is a huge advantage. PeterRabbit: yes, conjugate transposes are needed. The complex covariance matrix is Hermitean positive-definite. Like its real counterpart, it still has positive real eigenvalues, which addresses the question of meaningfulness. – whuber Jul 31 '13 at 15:18
• @whuber It doesn't make sense to me whatsoever to go into complex numbers if the problem is as shown. It is not simpler to deal with complex numbers - otherwise there wouldn't be a question here at all. Not everything will work fine with complex numbers and it's not a straightforward change if you don't know what you are doing. Transforming this problem in real space is equivalent, and you can apply all the variety of statistical techniques then without worrying if it works or not in complex space. – sashkello Jul 31 '13 at 23:30
• @whuber Good answer and nice explanation. I'd say as soon as you get over the transform from one to another it's really not hard... – sashkello Aug 1 '13 at 23:25

### Summary

The generalization of least-squares regression to complex-valued variables is straightforward, consisting primarily of replacing matrix transposes by conjugate transposes in the usual matrix formulas. A complex-valued regression, though, corresponds to a complicated multivariate multiple regression whose solution would be much more difficult to obtain using standard (real variable) methods. Thus, when the complex-valued model is meaningful, using complex arithmetic to obtain a solution is strongly recommended. This answer also includes some suggested ways to display the data and present diagnostic plots of the fit.

For simplicity, let's discuss the case of ordinary (univariate) regression, which can be written

$$z_j = \beta_0 + \beta_1 w_j + \varepsilon_j.$$

I have taken the liberty of naming the independent variable $W$ and the dependent variable $Z$, which is conventional (see, for instance, Lars Ahlfors, Complex Analysis). All that follows is straightforward to extend to the multiple regression setting.

### Interpretation

This model has an easily visualized geometric interpretation: multiplication by $\beta_1$ will rescale $w_j$ by the modulus of $\beta_1$ and rotate it around the origin by the argument of $\beta_1$. Subsequently, adding $\beta_0$ translates the result by this amount. The effect of $\varepsilon_j$ is to "jitter" that translation a little bit. Thus, regressing the $z_j$ on the $w_j$ in this manner is an effort to understand the collection of 2D points $(z_j)$ as arising from a constellation of 2D points $(w_j)$ via such a transformation, allowing for some error in the process. This is illustrated below with the figure titled "Fit as a Transformation."

Note that the rescaling and rotation are not just any linear transformation of the plane: they rule out skew transformations, for instance. Thus this model is not the same as a bivariate multiple regression with four parameters.

### Ordinary Least Squares

To connect the complex case with the real case, let's write

$z_j = x_j + i y_j$ for the values of the dependent variable and

$w_j = u_j + i v_j$ for the values of the independent variable.

Furthermore, for the parameters write

$\beta_0 = \gamma_0 + i \delta_0$ and $\beta_1 = \gamma_1 +i \delta_1$.

Every one of the new terms introduced is, of course, real, and $i^2 = -1$ is imaginary while $j=1, 2, \ldots, n$ indexes the data.

OLS finds $\hat\beta_0$ and $\hat\beta_1$ that minimize the sum of squares of deviations,

$$\sum_{j=1}^n ||z_j - \left(\hat\beta_0 + \hat\beta_1 w_j\right)||^2 = \sum_{j=1}^n \left(\bar z_j - \left(\bar{\hat\beta_0} + \bar{\hat\beta_1} \bar w_j\right)\right) \left(z_j - \left(\hat\beta_0 + \hat\beta_1 w_j\right)\right).$$

Formally this is identical to the usual matrix formulation: compare it to $\left(z - X\beta\right)'\left(z - X\beta\right).$ The only difference we find is that the transpose of the design matrix $X'$ is replaced by the conjugate transpose $X^* = \bar X '$. Consequently the formal matrix solution is

$$\hat\beta = \left(X^*X\right)^{-1}X^* z.$$

At the same time, to see what might be accomplished by casting this into a purely real-variable problem, we may write the OLS objective out in terms of the real components:

$$\sum_{j=1}^n \left(x_j-\gamma_0-\gamma_1u_j+\delta_1v_j\right)^2 + \sum_{j=1}^n\left(y_j-\delta_0-\delta_1u_j-\gamma_1v_j\right)^2.$$

Evidently this represents two linked real regressions: one of them regresses $x$ on $u$ and $v$, the other regresses $y$ on $u$ and $v$; and we require that the $v$ coefficient for $x$ be the negative of the $u$ coefficient for $y$ and the $u$ coefficient for $x$ equal the $v$ coefficient for $y$. Moreover, because the total squares of residuals from the two regressions are to be minimized, it will usually not be the case that either set of coefficients gives the best estimate for $x$ or $y$ alone. This is confirmed in the example below, which carries out the two real regressions separately and compares their solutions to the complex regression.

This analysis makes it apparent that rewriting the complex regression in terms of the real parts (1) complicates the formulas, (2) obscures the simple geometric interpretation, and (3) would require a generalized multivariate multiple regression (with nontrivial correlations among the variables) to solve. We can do better.

### Example

As an example, I take a grid of $w$ values at integral points near the origin in the complex plane. To the transformed values $w\beta$ are added iid errors having a bivariate Gaussian distribution: in particular, the real and imaginary parts of the errors are not independent.

It is difficult to draw the usual scatterplot of $(w_j, z_j)$ for complex variables, because it would consist of points in four dimensions. Instead we can view the scatterplot matrix of their real and imaginary parts.

Ignore the fit for now and look at the top four rows and four left columns: these display the data. The circular grid of $w$ is evident in the upper left; it has $81$ points. The scatterplots of the components of $w$ against the components of $z$ show clear correlations. Three of them have negative correlations; only the $y$ (the imaginary part of $z$) and $u$ (the real part of $w$) are positively correlated.

For these data, the true value of $\beta$ is $(-20 + 5i, -3/4 + 3/4\sqrt{3}i)$. It represents an expansion by $3/2$ and a counterclockwise rotation of 120 degrees followed by translation of $20$ units to the left and $5$ units up. I compute three fits: the complex least squares solution and two OLS solutions for $(x_j)$ and $(y_j)$ separately, for comparison.

Fit            Intercept          Slope(s)
True           -20    + 5 i       -0.75 + 1.30 i
Complex        -20.02 + 5.01 i    -0.83 + 1.38 i
Real only      -20.02             -0.75, -1.46
Imaginary only          5.01       1.30, -0.92


It will always be the case that the real-only intercept agrees with the real part of the complex intercept and the imaginary-only intercept agrees with the imaginary part fo the complex intercept. It is apparent, though, that the real-only and imaginary-only slopes neither agree with the complex slope coefficients nor with each other, exactly as predicted.

Let's take a closer look at the results of the complex fit. First, a plot of the residuals gives us an indication of their bivariate Gaussian distribution. (The underlying distribution has marginal standard deviations of $2$ and a correlation of $0.8$.) Then, we can plot the magnitudes of the residuals (represented by sizes of the circular symbols) and their arguments (represented by colors exactly as in the first plot) against the fitted values: this plot should look like a random distribution of sizes and colors, which it does.

Finally, we can depict the fit in several ways. The fit appeared in the last rows and columns of the scatterplot matrix (q.v.) and may be worth a closer look at this point. Below on the left the fits are plotted as open blue circles and arrows (representing the residuals) connect them to the data, shown as solid red circles. On the right the $(w_j)$ are shown as open black circles filled in with colors corresponding to their arguments; these are connected by arrows to the corresponding values of $(z_j)$. Recall that each arrow represents an expansion by $3/2$ around the origin, rotation by $120$ degrees, and translation by $(-20, 5)$, plus that bivariate Guassian error.

These results, the plots, and the diagnostic plots all suggest that the complex regression formula works correctly and achieves something different than separate linear regressions of the real and imaginary parts of the variables.

### Code

The R code to create the data, fits, and plots appears below. Note that the actual solution of $\hat\beta$ is obtained in a single line of code. Additional work--but not too much of it--would be needed to obtain the usual least squares output: the variance-covariance matrix of the fit, standard errors, p-values, etc.

#
# Synthesize data.
# (1) the independent variable w.
#
w.max <- 5 # Max extent of the independent values
w <- expand.grid(seq(-w.max,w.max), seq(-w.max,w.max))
w <- complex(real=w[[1]], imaginary=w[[2]])
w <- w[Mod(w) <= w.max]
n <- length(w)
#
# (2) the dependent variable z.
#
beta <- c(-20+5i, complex(argument=2*pi/3, modulus=3/2))
sigma <- 2; rho <- 0.8 # Parameters of the error distribution
library(MASS) #mvrnorm
set.seed(17)
e <- mvrnorm(n, c(0,0), matrix(c(1,rho,rho,1)*sigma^2, 2))
e <- complex(real=e[,1], imaginary=e[,2])
z <- as.vector((X <- cbind(rep(1,n), w)) %*% beta + e)
#
# Fit the models.
#
print(beta, digits=3)
print(beta.hat <- solve(Conj(t(X)) %*% X, Conj(t(X)) %*% z), digits=3)
print(beta.r <- coef(lm(Re(z) ~ Re(w) + Im(w))), digits=3)
print(beta.i <- coef(lm(Im(z) ~ Re(w) + Im(w))), digits=3)
#
# Show some diagnostics.
#
par(mfrow=c(1,2))
res <- as.vector(z - X %*% beta.hat)
fit <- z - res
s <- sqrt(Re(mean(Conj(res)*res)))
col <- hsv((Arg(res)/pi + 1)/2, .8, .9)
size <- Mod(res) / s
plot(res, pch=16, cex=size, col=col, main="Residuals")
plot(Re(fit), Im(fit), pch=16, cex = size, col=col,
main="Residuals vs. Fitted")

plot(Re(c(z, fit)), Im(c(z, fit)), type="n",
main="Residuals as Fit --> Data", xlab="Real", ylab="Imaginary")
points(Re(fit), Im(fit), col="Blue")
points(Re(z), Im(z), pch=16, col="Red")
arrows(Re(fit), Im(fit), Re(z), Im(z), col="Gray", length=0.1)

col.w <-  hsv((Arg(w)/pi + 1)/2, .8, .9)
plot(Re(c(w, z)), Im(c(w, z)), type="n",
main="Fit as a Transformation", xlab="Real", ylab="Imaginary")
points(Re(w), Im(w), pch=16, col=col.w)
points(Re(w), Im(w))
points(Re(z), Im(z), pch=16, col=col.w)
arrows(Re(w), Im(w), Re(z), Im(z), col="#00000030", length=0.1)
#
# Display the data.
#
par(mfrow=c(1,1))
pairs(cbind(w.Re=Re(w), w.Im=Im(w), z.Re=Re(z), z.Im=Im(z),
fit.Re=Re(fit), fit.Im=Im(fit)), cex=1/2)

• I had a follow up question, about the estimator $\hat{\beta}$ and its covariance. When I solve my problem with complex $y$, the covariance matrix (which I estimate using the fit residual) of my estimator has real and imaginary parts. I'm not sure how this works. Is the imaginary part of the covariance just concerning the imaginary part of the estimator (same for the real part)? If I want to plot CI's, I'm not sure how to go about this... Do the imaginary and real parts of the estimator have the same CI? Would it be possible to include a bit of info on this in your explanation? Thank you! – bill_e Sep 17 '13 at 4:49
• If all was computed correctly, the covariance will still be positive-definite. In particular, this implies that when you use it to compute the covariance of either the real part or the imaginary part of a variable, you will get a positive number, so all CI's will be well-defined. – whuber Sep 17 '13 at 15:01
• Cov matrix is positive semi-definite, but I guess what I am confused about is where you say: "when you use it to compute the covariance of either the real part or the imaginary part of a variable". I assumed that when I would compute a CI, it would have a real and imag part, which would correspond to the real and imag part of an element of $\hat{\beta}$. This doesn't seem to be the case though. Do you know why this is? – bill_e Sep 18 '13 at 5:59
• Also, if when I compute values for test statistics, I get numbers like say, 3 + .1*i. For this I was expecting the number to have no imaginary part. Is this normal? Or a sign I am doing something wrong? – bill_e Sep 18 '13 at 6:09
• When you compute test statistics with complex numbers, you should expect to get complex results! If you have a mathematical reason why the statistic should be real, then the calculation must be erroneous. When the imaginary part is really tiny compared to the real part, that's likely accumulated floating point error and it's usually safe to kill it off (zapsmall in R). Otherwise it is a sign something is fundamentally wrong. – whuber Sep 18 '13 at 13:35

After a nice long google sesh, I found some relevant information on understanding the problem in an alternative manner. It turns out similar problems are somewhat common in statistical signal processing. Instead of starting with a Gaussian likelihood which corresponds to a linear least squares for real data, one starts with a:

http://en.wikipedia.org/wiki/Complex_normal_distribution

Specifically, if you can assume that the distribution of your estimator $$\hat{\beta}$$ is multivariate normal, then in the case of complex data one would use the complex normal. The computation of the covariance of this estimator is a bit different, and given on the wiki page.

The textbook by Giri, Multivariate Statistical Analysis also covers this.

This issue has come up again on the Mathematica StackExchange and my answer/extended comment there is that @whuber 's excellent answer should be followed.

My answer here is an attempt to extend @whuber 's answer just a little bit by making the error structure a little more explicit. The proposed least squares estimator is what one would use if the bivariate error distribution has a zero correlation between the real and imaginary components. (But the data generated has a error correlation of 0.8.)

If one has access to a symbolic algebra program, then some of the messiness of constructing maximum likelihood estimators of the parameters (both the "fixed" effects and the covariance structure) can be eliminated. Below I use the same data as in @whuber 's answer and construct the maximum likelihood estimates by assuming $$\rho=0$$ and then by assuming $$\rho\neq0$$. I've used Mathematica but I suspect any other symbolic algebra program can do something similar. (And I've first posted a picture of the code and output followed by the actual code in an appendix as I can't get the Mathematica code to look as it should with just using text.)

Now for the maximum likelihood estimates assuming $$\rho=0$$...

We see that the maximum likelihood estimates which assume that $$\rho=0$$ match perfectly with the total least squares estimates.

Now let the data determine an estimate for $$\rho$$:

We see that $$\gamma_0$$ and $$\delta_0$$ are essentially identical whether or not we allow for the estimation of $$\rho$$. But $$\gamma_1$$ is much closer to the value that generated the data (although inferences with a sample size of 1 shouldn't be considered definitive to say the least) and the log of the likelihood is much higher.

My point in all of this is that the model being fit needs to be made completely explicit and that symbolic algebra programs can help alleviate the messiness. (And, of course, the maximum likelihood estimators assume a bivariate normal distribution which the least squares estimators do not assume.)

Appendix: The full Mathematica code

(* 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};

(* whuber 's least squares estimates *)
{a, b} = Inverse[ConjugateTranspose[w1].w1].ConjugateTranspose[w1].z
(* {-20.0172+5.00968 \[ImaginaryI],-0.830797+1.37827 \[ImaginaryI]} *)

(* Break up into the real and imaginary components *)
x = Re[z];
y = Im[z];
u = Re[w];
v = Im[w];
n = Length[z]; (* Sample size *)

(* Construct the real and imaginary components of the model *)
(* This is the messy part you probably don't want to do too often with paper and pencil *)
model = \[Gamma]0 + I \[Delta]0 + (\[Gamma]1 + I \[Delta]1) (u + I v);
modelR = Table[
Re[ComplexExpand[model[[j]]]] /. Im[h_] -> 0 /. Re[h_] -> h, {j, n}];
(* \[Gamma]0+u \[Gamma]1-v \[Delta]1 *)
modelI = Table[
Im[ComplexExpand[model[[j]]]] /. Im[h_] -> 0 /. Re[h_] -> h, {j, n}];
(* v \[Gamma]1+\[Delta]0+u \[Delta]1 *)

(* Construct the log of the likelihood as we are estimating the parameters associated with a bivariate normal distribution *)
logL = LogLikelihood[
BinormalDistribution[{0, 0}, {\[Sigma]1, \[Sigma]2}, \[Rho]],
Transpose[{x - modelR, y - modelI}]];

mle0 = FindMaximum[{logL /. {\[Rho] ->
0, \[Sigma]1 -> \[Sigma], \[Sigma]2 -> \[Sigma]}, \[Sigma] >
0}, {\[Gamma]0, \[Delta]0, \[Gamma]1, \[Delta]1, \[Sigma]}]
(* {-357.626,{\[Gamma]0\[Rule]-20.0172,\[Delta]0\[Rule]5.00968,\[Gamma]1\[Rule]-0.830797,\[Delta]1\[Rule]1.37827,\[Sigma]\[Rule]2.20038}} *)

(* Now suppose we don't want to restrict \[Rho]=0 *)
mle1 = FindMaximum[{logL /. {\[Sigma]1 -> \[Sigma], \[Sigma]2 -> \[Sigma]}, \[Sigma] > 0 && -1 < \[Rho] <
1}, {\[Gamma]0, \[Delta]0, \[Gamma]1, \[Delta]1, \[Sigma], \[Rho]}]
(* {-315.313,{\[Gamma]0\[Rule]-20.0172,\[Delta]0\[Rule]5.00968,\[Gamma]1\[Rule]-0.763237,\[Delta]1\[Rule]1.30859,\[Sigma]\[Rule]2.21424,\[Rho]\[Rule]0.810525}} *)


While @whuber has a beautifully-illustrated and well-explained answer, I think it's a simplified model that misses some of the power of the complex space.

Linear least-squares regression on reals is equivalent to the following model with inputs $w$, parameters $\beta$, and target $x$:

$$z = \beta_0 + \beta_1 w + \epsilon$$

where $\epsilon$ is normally-distributed with zero mean and some (typically constant) variance.

I suggest that complex linear regression be defined as follows:

$$z = \beta_0 + \beta_1 w + \beta_2 \overline w + \epsilon$$

There are two major differences.

First, there is an additional degree of freedom $\beta_2$ that allows phase sensitivity. You might not want that, but you can easily have that.

Second, $\epsilon$ is a complex normal distribution with zero mean and some variance and “pseudo-variance”.

Going back to the real model, the ordinary least squares solution comes out minimizing the loss, which is the negative log-likelihood. For a normal distribution, this is the parabola:

$$y = ax^2 + cx + d.$$

where $x = z - (\beta_0 + \beta_1 w)$, $a$ is fixed (typically), $c$ is zero as per the model, and $d$ doesn't matter since loss functions are invariant under constant addition.

Back to the complex model, the negative log-likelihood is \begin{align} y = a{|x|}^2 + \Re({bx^2 + cx}) + d. \end{align}

$c$ and $d$ are zero as before. $a$ is the curvature and $b$ is the “pseudo-curvature”. $b$ captures anisotropic components. If the $\Re$ function bothers you, then an equivalent way of writing this is \begin{align} {\begin{bmatrix}x-\mu \\ \overline{x-\mu}\end{bmatrix}}^H \begin{bmatrix}s & u \\ \overline{u} & \overline{s}\end{bmatrix}^{-1}\! \begin{bmatrix}x-\mu \\ \overline{x-\mu}\end{bmatrix} + d \end{align} for another set of parameters $s, u, \mu, d$. Here $s$ is the variance and $u$ is the pseudo-variance. $\mu$ is zero as per our model.

Here's an image of a complex normal distribution's density:

Notice how it's asymmetric. Without the $b$ parameter, it can't be asymmetric.

This complicates the regression although I'm pretty sure the solution is still analytical. I solved it for the case of one input, and I'm happy to transcribe my solution here, but I have a feeling that whuber might solve the general case.

• Thank you for this contribution. I don't follow it, though, because I'm not sure (a) why you introduce a quadratic polynomial, (b) what you actually mean by "corresponding" polynomial, or (c) what statistical model you are fitting. Would you be able to elaborate on those? – whuber Apr 10 '17 at 13:45
• @whuber I've rewritten it as a statistical model. Please let me know if makes sense to you. – Neil G Apr 10 '17 at 14:34
• Thank you: That clears it up (+1). Your model is no longer an analytic function of the variables. But because it is an analytic function of the parameters, it can be conceived of as a multiple regression of $z$ against the two complex variables $w$ and $\bar w$. In addition, you allow $\epsilon$ to have a more flexible distribution: that's not comprehended within my solution. As far as I can tell, your solution is equivalent to converting everything into its real and imaginary parts and conducting a multivariate multiple real regression. – whuber Apr 10 '17 at 14:49
• @whuber Right, with the two changes I suggested, I think it is as you said multivariate real regression. $\Beta_2$ can be removed to constrain the transformation as you describe in your solution. However, the pseudo-curvature term has some realistic practical applications such as trying to do regression to predict an AC voltage with a nonzero ground state? – Neil G Apr 10 '17 at 15:01
• It's true, Neil, that the loss is not analytic, but the solution is analytic. – whuber Apr 10 '17 at 15:03