Analiza ze złożonymi danymi, coś innego?

31

Powiedzmy na przykład, że robisz model liniowy, ale dane są złożone.y

y=xβ+ϵ

Mój zestaw danych jest złożony, ponieważ we wszystkich liczbach mają postać . Czy jest coś proceduralnie odmiennego podczas pracy z takimi danymi?y(a+bi)

Pytam, bo skończysz na otrzymywaniu złożonych macierzy kowariancji i statystykach testowych, które są złożone.

Czy podczas wykonywania najmniejszych kwadratów musisz użyć transpozycji sprzężonej zamiast transpozycji? czy złożona kowariancja ma znaczenie?

bill_e
źródło
3
Rozważ liczbę zespoloną jako dwie oddzielne zmienne, w ten sposób usuń i ze wszystkich równań. W przeciwnym razie będzie to koszmar ...
sashkello
Wszelkie informacje na temat lub ? xβ
Stijn
3
@Sashkello Jaki „koszmar”? Wymiary są zmniejszone o połowę, gdy używasz liczb zespolonych, więc prawdopodobnie jest to uproszczenie. Co więcej, zamieniłeś dwuwariantowy DV w jednoczynnikowy DV, co jest ogromną zaletą. PeterRabbit: tak, potrzebne są transpozycje sprzężone. Złożoną macierzą kowariancji jest Hermitean-pozytywnie określona. Podobnie jak jego prawdziwy odpowiednik, nadal ma pozytywne rzeczywiste wartości własne, które odnoszą się do pytania o sens.
whuber
2
@ whuber Nie ma dla mnie sensu wchodzenie w liczby zespolone, jeśli problem jest taki, jak pokazano. Nie jest łatwiej radzić sobie z liczbami zespolonymi - inaczej nie byłoby tutaj żadnego pytania. Nie wszystko będzie dobrze działać z liczbami zespolonymi i nie jest to prosta zmiana, jeśli nie wiesz, co robisz. Przekształcenie tego problemu w rzeczywistą przestrzeń jest równoważne i możesz zastosować całą różnorodność technik statystycznych, nie martwiąc się, czy to działa, czy nie w złożonej przestrzeni.
sashkello
1
@whuber Dobra odpowiedź i ładne wyjaśnienie. Powiedziałbym, że jak tylko przejdziesz przez transformację z jednej na drugą, to naprawdę nie jest trudne ...
sashkello

Odpowiedzi:

40

Podsumowanie

Uogólnienie regresji metodą najmniejszych kwadratów do zmiennych o wartościach zespolonych jest proste, polegające przede wszystkim na zastąpieniu transpozycji macierzy transpozycjami sprzężonymi w zwykłych formułach macierzy. Jednak regresja o złożonej wartości odpowiada skomplikowanej regresji wielowymiarowej wieloczynnikowej, której rozwiązanie byłoby znacznie trudniejsze do uzyskania przy użyciu standardowych metod (zmiennej rzeczywistej). Dlatego, gdy model o wartościach zespolonych ma sens, zdecydowanie zaleca się stosowanie złożonej arytmetyki w celu uzyskania rozwiązania. Ta odpowiedź zawiera także kilka sugerowanych sposobów wyświetlania danych i prezentacji wykresów diagnostycznych dopasowania.


Dla uproszczenia omówmy przypadek regresji zwykłej (jednoczynnikowej), którą można zapisać

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

Pozwoliłem sobie nazwać zmienną niezależną i zmienną zależną , która jest umowna (patrz na przykład Lars Ahlfors, Analiza złożona ). Wszystko, co następuje, można łatwo rozszerzyć na ustawienie regresji wielokrotnej.ZWZ

Interpretacja

Model ten łatwo uwidocznić interpretacji geometryczne: mnożenie przez będzie przeskalowanie przez moduł i obrócić ją wokół pochodzeniu argumentu . Następnie dodanie tłumaczy wynik o tę kwotę. Efektem jest „drżenie” tego tłumaczenia. Zatem regresowanie na w ten sposób jest próbą zrozumienia zbioru punktów 2D wynikającego z konstelacji punktów 2Dw j β 1 β 1 β 0 ε j z j w j ( z j ) ( w j )β1 wjβ1β1β0εjzjwj(zj)(wj)poprzez taką transformację, dopuszczając pewien błąd w procesie. Ilustruje to rysunek zatytułowany „Dopasuj jako transformację”.

Należy zauważyć, że przeskalowanie i obrót nie są po prostu żadną liniową transformacją płaszczyzny: wykluczają na przykład transformacje skośne. Zatem ten model nie jest tym samym co dwuwymiarowa regresja wielokrotna z czterema parametrami.

Zwykłe najmniejsze kwadraty

Aby połączyć złożoną sprawę ze sprawą prawdziwą, napiszmy

zj=xj+iyj dla wartości zmiennej zależnej i

wj=uj+ivj dla wartości zmiennej niezależnej.

Ponadto, dla parametrów napisz

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

Każdy z wprowadzonych nowych terminów jest oczywiście prawdziwy, a jest wyimaginowane, zaś indeksuje dane.j = 1 , 2 , , ni2=1j=1,2,,n

OLS znajduje i które minimalizują sumę kwadratów odchyleń, β 1β^0β^1

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

Formalnie jest to identyczne ze zwykłym sformułowaniem macierzowym: porównaj to z Jedyną różnicą, którą widzimy, jest to, że transpozycja macierzy projektowej jest zastąpiona sprzężoną transpozycją . W związku z tym formalnym rozwiązaniem macierzy jestX X = ˉ X(zXβ)(zXβ).X X=X¯

β^=(XX)1Xz.

Jednocześnie, aby zobaczyć, co można osiągnąć, umieszczając to w problemie o wyłącznie rzeczywistej zmiennej, możemy napisać cel OLS pod względem rzeczywistych składników:

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

Widocznie ta obejmuje dwie połączone rzeczywiste regresji: jeden z nich ulega zmniejszeniu o i , do pozostałych cofa o i ; i wymagamy, aby współczynnik dla był ujemny współczynnika dla a współczynnik dla równy współczynnik dla . Ponadto, ponieważ ogółemu v y u v v x u y u x v y x yxuvyuvvxuyuxvykwadraty reszt z dwóch regresji należy zminimalizować, zwykle nie będzie tak, że którykolwiek zestaw współczynników da najlepsze oszacowanie dla samego lub . Potwierdza to poniższy przykład, w którym oddzielnie przeprowadza się dwie prawdziwe regresje i porównuje ich rozwiązania z regresją złożoną.xy

Ta analiza pokazuje, że przepisanie złożonej regresji w kategoriach części rzeczywistych (1) komplikuje formuły, (2) przesłania prostą interpretację geometryczną i (3) wymagałoby uogólnionej regresji wieloczynnikowej wielorakiej (z nietrywialnymi korelacjami między zmiennymi ) rozwiązać. Możemy zrobić lepiej.

Przykład

Jako przykład biorę siatkę wartościach w integralnych punktów pobliżu pochodzenia w płaszczyźnie zespolonej. Do przekształconych wartości dodaje się błędy id mające dwuwymiarowy rozkład Gaussa: w szczególności rzeczywiste i urojone części błędów nie są niezależne.w βwwβ

Trudno jest narysować zwykły wykres rozproszenia dla zmiennych złożonych, ponieważ składałby się on z punktów w czterech wymiarach. Zamiast tego możemy zobaczyć matrycę wykresu rozrzutu ich rzeczywistych i urojonych części.(wj,zj)

Matryca punktowa

Na razie zignoruj ​​dopasowanie i spójrz na cztery górne wiersze i cztery lewe kolumny: wyświetlają one dane. Okrągła siatka widoczna jest w lewym górnym rogu; ma punktów. Wykresy rozrzutu składników względem składników wykazują wyraźne korelacje. Trzy z nich mają ujemne korelacje; tylko (urojona część ) (rzeczywista część ) są dodatnio skorelowane.81 w z y z u ww81wzyzuw

Dla tych danych prawdziwą wartością jest . Reprezentuje rozszerzenie o i obrót o 120 stopni w kierunku przeciwnym do ruchu wskazówek zegara, a następnie przesunięcie o jednostek w lewo i jednostek w górę. Obliczam trzy pasowania: złożone rozwiązanie najmniejszych kwadratów i dwa rozwiązania OLS dla i osobno, dla porównania.( - 20 + 5 i , - 3 / 4 + 3 / 4 β3/2205(xJ)(rj)(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

Zawsze będzie tak, że przechwytywanie tylko rzeczywiste zgadza się z rzeczywistą częścią przechwytywania złożonego, a przechwytywanie tylko wyobrażeniowe zgadza się z częścią urojoną przechwytywania złożonego. Oczywiste jest jednak, że zbocza tylko rzeczywiste i wyobrażone nie zgadzają się ze złożonymi współczynnikami nachylenia ani ze sobą, dokładnie tak, jak przewidywano.

Przyjrzyjmy się bliżej wynikom złożonego dopasowania. Po pierwsze, wykres reszt zawiera wskazanie ich dwuwymiarowego rozkładu Gaussa. (Podstawowy rozkład ma marginalne odchylenia standardowe i korelację .) Następnie możemy wykreślić wielkości reszt (reprezentowane przez rozmiary okrągłych symboli) i ich argumenty (reprezentowane przez kolory dokładnie tak, jak na pierwszym wykresie) w stosunku do dopasowanych wartości: ta fabuła powinna wyglądać jak losowy rozkład rozmiarów i kolorów, co robi.0,820.8

Działka resztkowa

Wreszcie możemy przedstawić dopasowanie na kilka sposobów. Dopasowanie pojawiło się w ostatnich wierszach i kolumnach macierzy wykresu rozrzutu ( qv ) i może być warte bliższego przyjrzenia się temu punktowi. Poniżej po lewej pasowania są wykreślone jako otwarte niebieskie kółka, a strzałki (reprezentujące resztki) łączą je z danymi, pokazanymi jako ciągłe czerwone kółka. Po prawej stronie są pokazane jako otwarte czarne kółka wypełnione kolorami odpowiadającymi ich argumentom; są one połączone strzałkami z odpowiednimi wartościami . Przypomnij sobie, że każda strzałka przedstawia rozszerzenie o wokół początku, obrót o stopni i tłumaczenie o , plus ten dwuwymiarowy błąd Guassiana.( z j ) 3 / 2 120 ( - 20 , 5 )(wj)(zj)3/2120(20,5)

Dopasuj jako transformację

Te wyniki, wykresy i wykresy diagnostyczne wszystkie sugerują, że formuła regresji złożonej działa poprawnie i osiąga coś innego niż oddzielne regresje liniowe rzeczywistych i urojonych części zmiennych.

Kod

RKod do tworzenia danych, drgawki, a działki znajduje się poniżej. Zauważ, że rzeczywiste rozwiązanie uzyskuje się w jednym wierszu kodu. Dodatkowa praca - ale nie za duża - byłaby potrzebna do uzyskania zwykłego wyniku najmniejszych kwadratów: macierzy wariancji-kowariancji dopasowania, błędów standardowych, wartości p itp.β^

#
# 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
źródło
Miałem kolejne pytanie dotyczące estymatora i jego kowariancji. Kiedy rozwiązuję mój problem ze złożonym , macierz kowariancji (którą szacuję za pomocą dopasowania resztkowego) mojego estymatora ma części rzeczywiste i urojone. Nie jestem pewien, jak to działa. Czy wyimaginowana część kowariancji dotyczy tylko wyimaginowanej części estymatora (tak samo w przypadku części rzeczywistej)? Jeśli chcę wykreślić CI, nie jestem pewien, jak sobie z tym poradzić ... Czy wyimaginowane i rzeczywiste części estymatora mają ten sam CI? Czy w wyjaśnieniu można by podać trochę informacji na ten temat? Dziękuję Ci! β^y
bill_e
Jeśli wszystko zostało poprawnie obliczone, kowariancja nadal będzie dodatnia. W szczególności oznacza to, że kiedy użyjesz go do obliczenia kowariancji części rzeczywistej lub części urojonej zmiennej, otrzymasz liczbę dodatnią, więc wszystkie CI będą dobrze zdefiniowane.
whuber
Macierz Cov jest dodatnia, półokreślona, ​​ale myślę, że jestem zdezorientowany w tym, co mówisz: „kiedy używasz jej do obliczania kowariancji albo części rzeczywistej, albo części urojonej zmiennej”. Zakładałem, że kiedy obliczę CI, będzie on miał część real i imag, która odpowiada rzeczywistej i imag części elementu . Wydaje się, że tak nie jest. Czy wiesz dlaczego tak jest? β^
bill_e
Ponadto, jeśli obliczę wartości dla statystyki testowej, otrzymam liczby takie jak powiedzmy 3 + .1 * i. W tym celu spodziewałem się, że liczba nie będzie miała części wyobrażonej. Czy to normalne? Czy znak, że robię coś źle?
bill_e
Kiedy obliczasz statystyki testowe z liczbami zespolonymi, powinieneś spodziewać się złożonych wyników! Jeśli masz matematyczny powód, dla którego statystyki powinny być prawdziwe, obliczenia muszą być błędne. Gdy część urojona jest naprawdę niewielka w porównaniu do części rzeczywistej, jest to prawdopodobnie nagromadzony błąd zmiennoprzecinkowy i zwykle bezpiecznie jest go zabić ( zapsmallw R). W przeciwnym razie jest to znak, że coś jest zasadniczo nie tak.
whuber
5

Po ładnym długim wyszukiwaniu w Google znalazłem kilka istotnych informacji na temat alternatywnego zrozumienia problemu. Okazuje się, że podobne problemy są dość powszechne w statystycznym przetwarzaniu sygnałów. Zamiast zaczynać od prawdopodobieństwa gaussowskiego, które odpowiada liniowym najmniejszym kwadratom dla rzeczywistych danych, zaczyna się od:

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

Ta strona Wikipedii daje zadowalające podsumowanie tego obiektu.

W szczególności, jeśli można założyć, że rozkład estymatora jest wieloramienny gaussowski, to w przypadku złożonych danych można użyć złożonej normalnej. Obliczenie kowariancji tego estymatora jest nieco inne i podane na stronie wiki. β^

Innym źródłem, które znalazłem, który dochodzi do tego samego wniosku co whuber, ale bada inne estymatory, takie jak maksymalne prawdopodobieństwo, to: „Oszacowania niewłaściwych modeli regresji liniowej”, autorstwa Yan i in.

bill_e
źródło
1

Podczas gdy @whuber ma pięknie zilustrowaną i dobrze wyjaśnioną odpowiedź, myślę, że jest to uproszczony model, w którym brakuje pewnej mocy złożonej przestrzeni.

Liniowa regresja najmniejszych kwadratów na liczbach rzeczywistych odpowiada poniższemu modelowi z danymi wejściowymi , parametrami i celem :wβx

z=β0+β1w+ϵ

gdzie jest normalnie rozłożony z zerową średnią i pewną (zwykle stałą) wariancją.ϵ

Proponuję zdefiniować złożoną regresję liniową w następujący sposób:

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

Istnieją dwie główne różnice.

Po pierwsze, istnieje dodatkowy stopień swobody który pozwala na czułość fazową. Możesz tego nie chcieć, ale możesz to łatwo mieć.β2

Po drugie, jest złożonym rozkładem normalnym o zerowej średniej i pewnej wariancji i „pseudo-wariancji”.ϵ

Wracając do prawdziwego modelu, wychodzi zwykłe rozwiązanie najmniejszych kwadratów minimalizujące straty, które są ujemnym prawdopodobieństwem logarytmicznym. Dla rozkładu normalnego jest to parabola:

y=ax2+cx+d.

gdzie , jest stałe (zwykle), jest zerowe jak w modelu, a nie ma znaczenia, ponieważ funkcje straty są niezmienne przy stałym dodawaniu.x=z(β0+β1w)acd

Wracając do modelu złożonego, ujemne prawdopodobieństwo logarytmu to

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

c i mają wartość zero, jak poprzednio. to krzywizna, a to „pseudo-krzywizna”. wychwytuje składniki anizotropowe. Jeśli przeszkadza Ci funkcja , równoważny sposób pisania to: dla innego zestawu parametrów . Oto wariancja, a pseudo-wariancja. wynosi zero zgodnie z naszym modelem.dabb

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

Oto obraz gęstości złożonego rozkładu normalnego:

Gęstość złożonego jednowymiarowego rozkładu normalnego

Zauważ, jak to jest asymetryczne. Bez parametru nie może być asymetryczny.b

To komplikuje regresję, chociaż jestem pewien, że rozwiązanie jest nadal analityczne. Rozwiązałem to dla przypadku jednego wejścia i cieszę się, że mogę tutaj przepisać moje rozwiązanie, ale mam wrażenie, że whuber może rozwiązać ogólny przypadek.

Neil G.
źródło
Dziękuję za ten wkład. Nie podążam za tym, ponieważ nie jestem pewien (a) dlaczego wprowadzasz kwadratowy wielomian, (b) co tak naprawdę rozumiesz przez „odpowiadający” wielomian, lub (c) jaki model statystyczny pasujesz. Czy byłbyś w stanie rozwinąć te kwestie?
whuber
@ whuber Przepisałem go jako model statystyczny. Daj mi znać, jeśli ma to dla Ciebie sens.
Neil G
Dziękuję: To wyjaśnia (+1). Twój model nie jest już funkcją analityczną zmiennych. Ponieważ jednak jest to funkcja analityczna parametrów, można ją traktować jako regresję wielokrotną względem dwóch zmiennych złożonych i . Ponadto zezwalasz na bardziej elastyczną dystrybucję: nie jest to zrozumiałe w moim rozwiązaniu. O ile mogę stwierdzić, twoje rozwiązanie jest równoważne przekształceniu wszystkiego w jego rzeczywiste i wymyślone części i przeprowadzeniu wielowymiarowej wielokrotnej regresji rzeczywistej . w ˉ w ϵzww¯ϵ
whuber
@ whuber Racja, z dwiema sugerowanymi przeze mnie zmianami, myślę, że jest tak, jak powiedziałeś regresję wielowymiarową. można usunąć, aby ograniczyć transformację zgodnie z opisem w rozwiązaniu. Jednak pseudo-krzywizna ma kilka realistycznych praktycznych zastosowań, takich jak próba regresji, aby przewidzieć napięcie prądu przemiennego o niezerowym stanie podstawowym? \Beta2
Neil G
Jeśli chodzi o funkcję analityczną, twoja nie jest ani analityczna, ponieważ twoja strata jest paraboloidem , co nie jest analityczne. Siodło jest analityczne, ale samo w sobie nie można go zminimalizować, ponieważ jest rozbieżne. x 2|x|2x2
Neil G
1

Ten problem pojawił się ponownie w Mathematica StackExchange, a moja odpowiedź / rozszerzony komentarz jest taki, że należy postępować zgodnie z doskonałą odpowiedzią @whuber.

Moja odpowiedź tutaj jest próbą rozszerzenia odpowiedzi @whuber tylko przez uszczegółowienie struktury błędu. Proponowany estymator najmniejszych kwadratów byłby wykorzystany, gdyby dwuwymiarowy rozkład błędów miał zerową korelację między składową rzeczywistą a urojoną. (Ale wygenerowane dane mają korelację błędów wynoszącą 0,8.)

Jeśli ktoś ma dostęp do programu algebry symbolicznej, to można wyeliminować część bałaganu związanego z konstruowaniem estymatorów maksymalnego prawdopodobieństwa parametrów (zarówno efektów „ustalonych”, jak i struktury kowariancji). Poniżej używam tych samych danych, co w odpowiedzi @whuber i konstruuję szacunki maksymalnego prawdopodobieństwa, przyjmując a następnie przyjmując . Użyłem Mathematiki, ale podejrzewam, że każdy inny program algebry symbolicznej może zrobić coś podobnego. (I najpierw opublikowałem zdjęcie kodu i wyniku, a następnie rzeczywisty kod w dodatku, ponieważ nie mogę sprawić, by kod Mathematica wyglądał tak, jak powinien, używając tylko tekstu.)ρ=0ρ0

Estymator danych i najmniejszych kwadratów

Teraz dla maksymalnych oszacowań prawdopodobieństwa przy założeniu ...ρ=0

szacunki maksymalnego prawdopodobieństwa przy założeniu, że rho wynosi zero

Widzimy, że szacunki maksymalnego prawdopodobieństwa, które zakładają, że idealnie pasują do oszacowań całkowitej liczby najmniejszych kwadratów.ρ=0

Teraz pozwól danym określić szacunkową wartość dla :ρ

Oszacowania maksymalnego prawdopodobieństwa, w tym rho

Widzimy, że i są zasadniczo identyczne bez względu na to, czy zezwalamy na oszacowanie . Ale jest znacznie bliższa wartości, która wygenerowała dane (chociaż wnioskowania o wielkości próby 1 nie powinny być uważane za ostateczne co najmniej), a log prawdopodobieństwa jest znacznie wyższy.γ0δ0ργ1

Chodzi mi o to, że dopasowanie modelu musi być całkowicie jednoznaczne, a symboliczne programy algebry mogą pomóc złagodzić bałagan. (I oczywiście estymatory maksymalnego prawdopodobieństwa zakładają dwuwymiarowy rozkład normalny, którego nie przyjmują estymatory najmniejszych kwadratów.)

Dodatek: Pełny kod 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
źródło