Análise com dados complexos, algo diferente?

31

Digamos, por exemplo, que você esteja criando um modelo linear, mas os dados são complexos.y

y=xβ+ϵ

Meu conjunto de dados é complexo, pois todos os números em têm a forma ( a + b i ) . Existe algo processualmente diferente ao trabalhar com esses dados?y(a+bi)

Eu pergunto porque, você acabará recebendo matrizes de covariância complexas e estatísticas de teste com valores complexos.

Você precisa usar um transposto conjugado em vez de transposto ao fazer menos quadrados? uma covariância complexa com valor é significativa?

bill_e
fonte
3
Considere um número complexo como duas variáveis ​​separadas e, dessa forma, remova i de todas as suas equações. Caso contrário, será um pesadelo ...
sashkello
Alguma informação sobre ou β ? xβ
Stijn
3
@Sashkello Que "pesadelo"? As dimensões são reduzidas pela metade quando você usa números complexos; portanto, isso é uma simplificação. Além disso, você transformou um DV bivariado em um DV univariado , o que é uma grande vantagem. PeterRabbit: sim, são necessárias transposições conjugadas. A matriz de covariância complexa é definida positivamente por eremita. Como sua contraparte real, ele ainda possui autovalores reais positivos, que abordam a questão da significância.
whuber
2
@whuber Não faz sentido para mim entrar em números complexos se o problema for como mostrado. Não é mais simples lidar com números complexos - caso contrário, não haveria uma pergunta aqui. Nem tudo funcionará bem com números complexos e não será uma mudança direta se você não souber o que está fazendo. Transformar esse problema no espaço real é equivalente , e você pode aplicar toda a variedade de técnicas estatísticas sem se preocupar se funciona ou não em um espaço complexo.
Sashkello 31/07
1
Boa resposta e boa explicação. Eu diria que, logo que você começa sobre a transformação de um para outro não é realmente difícil ...
sashkello

Respostas:

40

Sumário

A generalização da regressão dos mínimos quadrados para variáveis ​​de valor complexo é direta, consistindo principalmente na substituição de transposições de matriz por transpostas de conjugado nas fórmulas usuais de matriz. Uma regressão de valor complexo, no entanto, corresponde a uma regressão múltipla multivariada complicada cuja solução seria muito mais difícil de obter usando métodos padrão (variáveis ​​reais). Assim, quando o modelo de valor complexo é significativo, é altamente recomendável usar aritmética complexa para obter uma solução. Esta resposta também inclui algumas maneiras sugeridas de exibir os dados e apresentar gráficos de diagnóstico do ajuste.


Para simplificar, vamos discutir o caso da regressão ordinária (univariada), que pode ser escrita

zj=β0+β1wj+εj.

Tomei a liberdade de nomear a variável independente e a variável dependente Z , que é convencional (veja, por exemplo, Lars Ahlfors, Complex Analysis ). Tudo o que se segue é simples de se estender à configuração de regressão múltipla.WZ

Interpretação

Este modelo tem uma interpretação geométrica facilmente visualizada: a multiplicação por irá redimensionar w j pelo módulo de β 1 e girá- lo em torno da origem pelo argumento de β 1 . Posteriormente, adicionar β 0 traduz o resultado por esse valor. O efeito de ε j é "tremer" essa tradução um pouco. Assim, a regredir o z j na w j desta maneira é um esforço para compreender a recolha de pontos em 2D ( z j )β1 wjβ1β1β0εjzjwj(zj)como resultante de uma constelação de pontos 2D por meio de uma tal transformação, permitindo algum erro no processo. Isso é ilustrado abaixo com a figura intitulada "Ajustar como transformação".(wj)

Observe que o redimensionamento e a rotação não são apenas transformações lineares do plano: elas excluem transformações inclinadas, por exemplo. Portanto, esse modelo não é o mesmo que uma regressão múltipla bivariada com quatro parâmetros.

Mínimos Quadrados Ordinários

Para conectar o caso complexo ao caso real, vamos escrever

para os valores da variável dependente ezj=xj+iyj

para os valores da variável independente.wj=uj+ivj

Além disso, para os parâmetros escreva

e β 1 = γ 1 + i δ 1 . β0=γ0+iδ0β1=γ1+iδ1

Todos os novos termos introduzidos são, obviamente, reais, e é imaginário enquanto j = 1 , 2 , , n indexa os dados.i2=1j=1,2,,n

MQO achados p 0 e β 1 que minimizam a soma dos quadrados dos desvios,β^0β^1

j=1n||zj(β^0+β^1wj)||2=j=1n(z¯j(β^0¯+β^1¯w¯j))(zj(β^0+β^1wj)).

Formalmente, isso é idêntico à formulação usual da matriz: compare com A única diferença que encontramos é que a transposição da matriz de projeto X é substituída pela transposta conjugada X = ˉ X . Consequentemente, a solução formal de matriz é(zXβ)(zXβ).X X=X¯

β^=(XX)1Xz.

Ao mesmo tempo, para ver o que pode ser conseguido ao converter isso em um problema de variável puramente real, podemos escrever o objetivo do OLS em termos dos componentes reais:

j=1n(xjγ0γ1uj+δ1vj)2+j=1n(yjδ0δ1ujγ1vj)2.

xuvyuvvxuyuxvyxy

Essa análise torna aparente que reescrever a regressão complexa em termos das partes reais (1) complica as fórmulas, (2) obscurece a interpretação geométrica simples e (3) exigiria uma regressão múltipla multivariada generalizada (com correlações não triviais entre as variáveis ) resolver. Nós podemos fazer melhor.

Exemplo

wwβ

(wj,zj)

Matriz de dispersão

w81wzyzuw

β(20+5i,3/4+3/43i)3/2205(xj)(yj)

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

Será sempre o caso de que a interceptação somente real concorda com a parte real da interceptação complexa e a interceptação somente imaginária concorda com a parte imaginária da interceptação complexa. É aparente, no entanto, que as inclinações reais e imaginárias não concordam com os coeficientes de inclinação complexos nem entre si, exatamente como previsto.

20.8

Lote residual

(wj)(zj)3/2120(20,5)

Ajustar como transformação

Esses resultados, os gráficos e os diagnósticos sugerem que a fórmula de regressão complexa funciona corretamente e alcança algo diferente de regressões lineares separadas das partes real e imaginária das variáveis.

Código

Rβ^

#
# 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)
whuber
fonte
β^y
Se tudo foi computado corretamente, a covariância ainda será positiva-definida. Em particular, isso implica que, quando você o usa para calcular a covariância da parte real ou da parte imaginária de uma variável, você obtém um número positivo, portanto todos os ICs serão bem definidos.
whuber
β^
Além disso, se quando eu calcular valores para estatísticas de teste, obter números como, por exemplo, 3 + .1 * i. Para isso, esperava que o número não tivesse parte imaginária. Isso é normal? Ou um sinal de que estou fazendo algo errado?
bill_e
Ao calcular estatísticas de teste com números complexos, você deverá obter resultados complexos! Se você tem um motivo matemático para que a estatística seja real, o cálculo deve ser incorreto. Quando a parte imaginária é realmente pequena em comparação à parte real, é provável que ocorra um erro de ponto flutuante acumulado e geralmente é seguro eliminá-la ( zapsmallin R). Caso contrário, é um sinal de que algo está fundamentalmente errado.
whuber
5

Após uma longa e agradável pesquisa no google, encontrei algumas informações relevantes sobre como entender o problema de uma maneira alternativa. Acontece que problemas semelhantes são comuns no processamento estatístico de sinais. Em vez de começar com uma probabilidade gaussiana que corresponde a mínimos quadrados lineares para dados reais, começa-se com a:

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

Esta página da Wikipedia fornece um resumo satisfatório sobre esse objeto.

β^

Outra fonte que encontrei que chega à mesma conclusão que o whuber, mas explora outros estimadores como a probabilidade máxima é: "Estimativas dos modelos de regressão linear inadequada", de Yan et al.

bill_e
fonte
1

Embora o @whuber tenha uma resposta bem ilustrada e bem explicada, acho que é um modelo simplificado que perde parte do poder do espaço complexo.

wβx

z=β0+β1w+ϵ

ϵ

Sugiro que a regressão linear complexa seja definida da seguinte forma:

z=β0+β1w+β2w¯+ϵ

Existem duas grandes diferenças.

β2

ϵ

Voltando ao modelo real, a solução ordinária de mínimos quadrados sai minimizando a perda, que é a probabilidade logarítmica negativa. Para uma distribuição normal, esta é a parábola:

y=ax2+cx+d.

x=z(β0+β1w)acd

y=a|x|2+(bx2+cx)+d.

cdabb

[xμxμ¯]H[suu¯s¯]1[xμxμ¯]+d
s,u,μ,dsuμ

Aqui está uma imagem da densidade de uma distribuição normal complexa:

A densidade de uma distribuição normal univariada complexa

b

Isso complica a regressão, embora eu tenha certeza de que a solução ainda é analítica. Eu o resolvi para o caso de uma entrada e estou feliz em transcrever minha solução aqui, mas tenho a sensação de que o whuber pode resolver o caso geral.

Neil G
fonte
Obrigado por esta contribuição. Mas não o sigo, porque não tenho certeza (a) por que você introduz um polinômio quadrático, (b) o que você realmente quer dizer com polinômio "correspondente" ou (c) qual modelo estatístico você está ajustando. Você seria capaz de elaborar sobre isso?
whuber
@whuber Eu reescrevi como um modelo estatístico. Por favor, deixe-me saber se faz sentido para você.
Neil G
zww¯ϵ
\Beta2
|x|2x2
1

Esse problema surgiu novamente no Mathematica StackExchange e minha resposta / comentário extenso é que a resposta excelente do @whuber deve ser seguida.

Minha resposta aqui é uma tentativa de estender a resposta do @whuber um pouco, tornando a estrutura do erro um pouco mais explícita. O estimador de mínimos quadrados proposto é o que se usaria se a distribuição de erro bivariada tivesse uma correlação zero entre os componentes reais e imaginários. (Mas os dados gerados têm uma correlação de erro de 0,8.)

ρ=0ρ0

Estimador de dados e mínimos quadrados

ρ=0

estimativas de probabilidade máxima assumindo que rho é zero

ρ=0

ρ

Estimativas de probabilidade máxima, incluindo rho

γ0δ0ργ1

O que quero dizer com tudo isso é que o modelo que está sendo ajustado precisa ser completamente explícito e que os programas de álgebra simbólica podem ajudar a aliviar a confusão. (E, é claro, os estimadores de probabilidade máxima assumem uma distribuição normal bivariada que os estimadores de mínimos quadrados não assumem.)

Apêndice: O código completo do 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}} *)
JimB
fonte