Análisis con datos complejos, ¿algo diferente?

Digamos, por ejemplo, que está haciendo un modelo lineal, pero los datos son complejos.$y$

$y = x \beta + \epsilon$

Mi conjunto de datos es complejo, ya que en todos los números en son de la forma ( a + b i ) . ¿Hay algo procesalmente diferente cuando se trabaja con dichos datos?$y$$(a + bi)$

Pregunto porque terminarás obteniendo matrices de covarianza complejas y pruebas estadísticas que son de valor complejo.

¿Necesitas usar transposiciones conjugadas en lugar de transposiciones cuando haces mínimos cuadrados? ¿Tiene sentido una covarianza compleja valorada?

Resumen

La generalización de la regresión de mínimos cuadrados a variables de valores complejos es sencilla, y consiste principalmente en reemplazar las transposiciones matriciales por transposiciones conjugadas en las fórmulas matriciales habituales. Sin embargo, una regresión de valor complejo corresponde a una regresión múltiple multivariada complicada cuya solución sería mucho más difícil de obtener utilizando métodos estándar (variable real). Por lo tanto, cuando el modelo de valor complejo es significativo, se recomienda utilizar la aritmética compleja para obtener una solución. Esta respuesta también incluye algunas formas sugeridas para mostrar los datos y presentar gráficas de diagnóstico del ajuste.

---

Para simplificar, analicemos el caso de la regresión ordinaria (univariante), que puede escribirse

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

Me he tomado la libertad de nombrar la variable independiente y la variable dependiente Z , que es convencional (ver, por ejemplo, Lars Ahlfors, Complex Analysis ). Todo lo que sigue es sencillo de extender a la configuración de regresión múltiple.$W$$Z$

Interpretación

Este modelo tiene una interpretación geométrica fácilmente visualizado: multiplicación por se reescalar w j por el módulo de β 1 y girar alrededor del origen por el argumento de β 1 . Posteriormente, la adición de β 0 traduce el resultado en esta cantidad. El efecto de ε j es "fluctuar" esa traducción un poco. Por lo tanto, hacer retroceder la z j en la w j de esta manera es un esfuerzo por comprender la colección de puntos 2D ( z j$\beta_1$ $w_j$$\beta_1$$\beta_1$$\beta_0$$\varepsilon_j$$z_j$$w_j$$(z_j)$como resultado de una constelación de puntos 2D través de dicha transformación, lo que permite algún error en el proceso. Esto se ilustra a continuación con la figura titulada "Ajustar como transformación".$(w_j)$

Tenga en cuenta que el cambio de escala y la rotación no son simplemente una transformación lineal del plano: descartan transformaciones asimétricas, por ejemplo. Por lo tanto, este modelo no es lo mismo que una regresión múltiple bivariada con cuatro parámetros.

Mínimos cuadrados ordinarios

Para conectar el caso complejo con el caso real, escribamos

para los valores de la variable dependiente y$z_j = x_j + i y_j$

para los valores de la variable independiente.$w_j = u_j + i v_j$

Además, para los parámetros escriba

y β 1 = γ 1 + i δ 1 . $\beta_0 = \gamma_0 + i \delta_0$$\beta_1 = \gamma_1 +i \delta_1$

Cada uno de los nuevos términos introducidos es, por supuesto, real, e es imaginario, mientras que j = 1 , 2 , ... , n indexa los datos.$i^2 = -1$$j=1, 2, \ldots, n$

MCO hallazgos ß 0 y β 1 que minimizan la suma de los cuadrados de las desviaciones,$\hat\beta_0$$\hat\beta_1$

$$\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).$$

Formalmente, esto es idéntico a la formulación matricial habitual: compárela con La única diferencia que encontramos es que la transposición de la matriz de diseño X ' se reemplaza por la transposición conjugada X ∗ = ˉ X ′ . En consecuencia, la solución de matriz formal es$\left(z - X\beta\right)'\left(z - X\beta\right).$$X'$ $X^* = \bar X '$

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

Al mismo tiempo, para ver qué se puede lograr al convertir esto en un problema puramente variable real, podemos escribir el objetivo OLS en términos de los componentes reales:

$$\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.$$

Evidentemente, esto representa dos regresiones reales vinculadas : una de ellas regresa en u y v , la otra regresa y en u y v ; y requerimos que el coeficiente v para x sea ​​el negativo del coeficiente u para y y el coeficiente u para x sea igual al coeficiente v para y . Además, porque el total$x$$u$$v$$y$$u$$v$$v$$x$$u$$y$$u$$x$$v$$y$los cuadrados de los residuos de las dos regresiones se deben minimizar, por lo general no será el caso de que cualquiera de los conjuntos de coeficientes proporcione la mejor estimación para o y solo. Esto se confirma en el siguiente ejemplo, que realiza las dos regresiones reales por separado y compara sus soluciones con la regresión compleja.$x$$y$

Este análisis hace evidente que reescribir la regresión compleja en términos de las partes reales (1) complica las fórmulas, (2) oscurece la interpretación geométrica simple y (3) requeriría una regresión múltiple multivariada generalizada (con correlaciones no triviales entre las variables ) resolver. Podemos hacerlo mejor.

Ejemplo

Como ejemplo, tomo una cuadrícula de valores de en puntos integrales cerca del origen en el plano complejo. A los valores transformados w \ beta se agregan los errores iid que tienen una distribución gaussiana bivariada: en particular, las partes real e imaginaria de los errores no son independientes.$w$$w\beta$

Es difícil dibujar el diagrama de dispersión habitual de para variables complejas, ya que consistiría en puntos en cuatro dimensiones. En cambio, podemos ver la matriz de diagrama de dispersión de sus partes reales e imaginarias.$(w_j, z_j)$

Ignore el ajuste por ahora y mire las cuatro filas superiores y las cuatro columnas izquierdas: estas muestran los datos. La cuadrícula circular de es evidente en la esquina superior izquierda; Tiene 81 puntos. Los diagramas de dispersión de los componentes de w contra los componentes de z muestran correlaciones claras. Tres de ellos tienen correlaciones negativas; solo la y (la parte imaginaria de z ) yu (la parte real de w ) están positivamente correlacionadas.$w$$81$$w$$z$$y$$z$$u$$w$

Para estos datos, el verdadero valor de es ( - 20 + 5 i , - 3 / 4 + 3 / 4 $\beta$. Representa una expansión por3/2y una rotación en sentido antihorario de 120 grados, seguidos por la traducción de20unidades a la izquierda y5unidades hacia arriba. Calculo tres ajustes: la solución de mínimos cuadrados complejos y dos soluciones OLS para(xj)e(yj) porseparado, para comparar.$(-20 + 5i, -3/4 + 3/4\sqrt{3}i)$$3/2$$20$$5$$(x_j)$$(y_j)$

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

Siempre será el caso que la intersección de solo real esté de acuerdo con la parte real de la intersección compleja y la intersección de solo imaginario esté de acuerdo con la parte imaginaria de la intersección compleja. Sin embargo, es evidente que las pendientes solo real e imaginaria no están de acuerdo con los coeficientes de pendiente complejos ni entre sí, exactamente como se predijo.

Echemos un vistazo más de cerca a los resultados del ajuste complejo. Primero, una gráfica de los residuos nos da una indicación de su distribución gaussiana bivariada. (La distribución subyacente tiene desviaciones estándar marginales de y una correlación de 0.8 .) Luego, podemos trazar las magnitudes de los residuos (representados por tamaños de los símbolos circulares) y sus argumentos (representados por colores exactamente como en el primer diagrama) en contra de los valores ajustados: este gráfico debe verse como una distribución aleatoria de tamaños y colores, lo que hace.$2$$0.8$

Finalmente, podemos representar el ajuste de varias maneras. El ajuste apareció en las últimas filas y columnas de la matriz de diagrama de dispersión ( qv ) y puede merecer una mirada más cercana en este punto. A continuación, a la izquierda, los ajustes se trazan como círculos azules abiertos y las flechas (que representan los residuos) los conectan a los datos, que se muestran como círculos rojos sólidos. A la derecha, los se muestran como círculos negros abiertos rellenos con los colores correspondientes a sus argumentos; estos están conectados por flechas a los valores correspondientes de ( z j ) . Recordemos que cada flecha representa una expansión por 3 / 2 alrededor del origen, la rotación por $(w_j)$$(z_j)$$3/2$$120$grados, y traducción por , más ese error bivariado de Guassian.$(-20, 5)$

Estos resultados, las gráficas y las gráficas de diagnóstico sugieren que la fórmula de regresión compleja funciona correctamente y logra algo diferente que las regresiones lineales separadas de las partes real e imaginaria de las variables.

Código

El Rcódigo para crear los datos, ajustes y trazados aparece a continuación. Tenga en cuenta que la solución real de β se obtiene en una sola línea de código. Se necesitaría trabajo adicional, pero no demasiado, para obtener la salida de mínimos cuadrados habitual: la matriz de varianza-covarianza del ajuste, errores estándar, valores p, etc.$\hat\beta$

#
# 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)

Después de un largo y largo google sesh, encontré información relevante para entender el problema de una manera alternativa. Resulta que problemas similares son algo comunes en el procesamiento de señales estadísticas. En lugar de comenzar con una probabilidad gaussiana que corresponde a mínimos cuadrados lineales para datos reales, uno comienza con un:

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

Esta página de wikipedia ofrece un resumen satisfactorio de este objeto.

$\hat{\beta}$

Otra fuente que encontré que llega a la misma conclusión que Whuber, pero explora otros estimadores como la máxima probabilidad es: "Estimaciones de los modelos de regresión lineal inadecuados", de Yan et al.

Si bien @whuber tiene una respuesta bellamente ilustrada y bien explicada, creo que es un modelo simplificado que pierde parte del poder del espacio complejo.

$w$$\beta$$x$

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

$\epsilon$

Sugiero que la regresión lineal compleja se defina de la siguiente manera:

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

Hay dos diferencias principales.

$\beta_2$

$\epsilon$

Volviendo al modelo real, la solución de mínimos cuadrados ordinarios sale minimizando la pérdida, que es la probabilidad de registro negativa. Para una distribución normal, esta es la parábola:

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

$x = z - (\beta_0 + \beta_1 w)$$a$$c$$d$

$$\begin{align} y = a{|x|}^2 + \Re({bx^2 + cx}) + d. \end{align}$$

$c$$d$$a$$b$$b$$\Re$

$$\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}$$ $s, u, \mu, d$$s$$u$$\mu$

Aquí hay una imagen de la densidad de una distribución normal compleja:

$b$

Esto complica la regresión, aunque estoy bastante seguro de que la solución sigue siendo analítica. Lo resolví para el caso de una entrada, y estoy feliz de transcribir mi solución aquí, pero tengo la sensación de que Whuber podría resolver el caso general.

Este problema ha vuelto a aparecer en el Mathematica StackExchange y mi respuesta / comentario extendido es que se debe seguir la excelente respuesta de @whuber.

Mi respuesta aquí es un intento de extender la respuesta de @whuber solo un poco haciendo que la estructura del error sea un poco más explícita. El estimador de mínimos cuadrados propuesto es el que se usaría si la distribución de error bivariada tuviera una correlación cero entre los componentes real e imaginario. (Pero los datos generados tienen una correlación de errores de 0.8).

$\rho=0$$\rho\neq0$

$\rho=0$

$\rho=0$

$\rho$

$\gamma_0$$\delta_0$$\rho$$\gamma_1$

Mi punto en todo esto es que el modelo que se ajusta debe hacerse completamente explícito y que los programas de álgebra simbólica pueden ayudar a aliviar el desorden. (Y, por supuesto, los estimadores de máxima verosimilitud asumen una distribución normal bivariada que los estimadores de mínimos cuadrados no asumen).

Apéndice: el código completo de Mathematica

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