Z perspektywy statystycznej: transformata Fouriera vs regresja na podstawie Fouriera

13

Próbuję zrozumieć, czy dyskretna transformata Fouriera daje taką samą reprezentację krzywej jak regresja z wykorzystaniem zasady Fouriera. Na przykład,

library(fda)
Y=daily$tempav[,1] ## my data
length(Y) ## =365

## create Fourier basis and estimate the coefficients
mybasis=create.fourier.basis(c(0,365),365)  
basisMat=eval.basis(1:365,mybasis)
regcoef=coef(lm(Y~basisMat-1))

## using Fourier transform
fftcoef=fft(Y)

## compare
head(fftcoef)
head(regcoef)

FFT daje liczbę zespoloną, podczas gdy regresja daje liczbę rzeczywistą.

Czy przekazują te same informacje? Czy istnieje mapa jeden do jednego między dwoma zestawami liczb?

(Byłbym wdzięczny, gdyby odpowiedź została napisana z perspektywy statystyka zamiast z perspektywy inżyniera. Wiele materiałów online, które mogę znaleźć, ma żargon techniczny w całym miejscu, co czyni mnie mniej smakowitym.)

qoheleth
źródło
Nie znam Twojego fragmentu kodu, więc nie mogę powiedzieć, czy dotyczy on następującego problemu. Jednak zazwyczaj podstawa DFT jest definiowana w kategoriach częstotliwości całkowych („liczba całkowita”), podczas gdy ogólna „podstawa czterokierunkowa” dla regresji mogłaby wykorzystywać dowolne stosunki częstotliwości (np. Włączając nieracjonalne, przynajmniej w ciągłej arytmetyce). Może to również być interesujące.
GeoMatt22,
Myślę, że każdy skorzystałby, gdybyś napisał pytanie w kategoriach matematycznych (w przeciwieństwie do fragmentów kodu). Jaki problem regresji rozwiązujesz? Z jakich podstawowych funkcji Fouriera korzystasz? Będziesz zaskoczony, jak poprawią się odpowiedzi na twoje pytanie.
Yair Daon

Odpowiedzi:

15

Oni są tacy sami. Oto jak...

Robić regresję

Załóżmy, że do modelu gdzie i . Nie jest to jednak odpowiednie dla regresji liniowej, więc zamiast tego używasz trygonometrii ( ) i dopasuj model równoważny: Uruchamianie regresji liniowej na wszystkich częstotliwościach Fouriera daje ci kilka ( ) bet: , . Dla każdego

yt=j=1nAjcos(2πt[j/N]+ϕj)
t=1,,Nn=floor(N/2)cos(a+b)=cos(a)cos(b)sin(a)sin(b)
yt=j=1nβ1,jcos(2πt[j/N])+β2,jsin(2πt[j/N]).
{j/N:j=1,,n}2n{β^i,j}i=1,2j, jeśli chcesz ręcznie obliczyć parę, możesz użyć:

β^1,j=t=1Nytcos(2πt[j/N])t=1Ncos2(2πt[j/N])
i Są to standardowe formuły regresji.
β^2,j=t=1Nytsin(2πt[j/N])t=1Nsin2(2πt[j/N]).

Wykonywanie dyskretnej transformaty Fouriera

Po uruchomieniu transformacji Fouriera obliczasz dla :j=1,,n

d(j/N)=N1/2t=1Nytexp[2πit[j/N]]=N1/2(t=1Nytcos(2πt[j/N])it=1Nytsin(2πt[j/N])).

Jest to liczba zespolona (zwróć uwagę na ). Aby zobaczyć, dlaczego ta równość się utrzymuje, należy pamiętać, że , i .ieix=cos(x)+isin(x)cos(x)=cos(x)sin(x)=sin(x)

Dla każdego biorąc kwadrat złożonej koniugatu daje „ periodogram :”j

|d(j/N)|2=N1(t=1Nytcos(2πt[j/N]))2+N1(t=1Nytsin(2πt[j/N]))2.
W R obliczenie tego wektora byłoby I <- abs(fft(Y))^2/length(Y)dziwne, ponieważ trzeba go skalować.

Możesz także zdefiniować „ skalowany periodogram ” Wyraźnie . W R byłoby to .

P(j/N)=(2Nt=1Nytcos(2πt[j/N]))2+(2Nt=1Nytsin(2πt[j/N]))2.
P(j/N)=4N|d(j/N)|2P <- (4/length(Y))*I[(1:floor(length(Y)/2))]

Połączenie między dwoma

Okazuje się, że związek między regresją a dwoma periodogramami jest:

j N t =

P(j/N)=β^1,j2+β^2,j2.
Dlaczego? Ponieważ wybrana podstawa jest ortogonalna / ortonormalna. Możesz pokazać dla każdego że . Podłącz to do mianowników formuł dla współczynników regresji i voila.jt=1Ncos2(2πt[j/N])=t=1Nsin2(2πt[j/N])=N/2

Źródło: https://www.amazon.com/Time-Analysis-Its-Applications-Statistics/dp/144197864X

Taylor
źródło
1
+1 za odpowiedź i źródło. Byłoby również dobrze, jeśli możesz zademonstrować wynik za pomocą Robiektów, które opublikowałem.
qoheleth
@ qoheleth Zostawię to tobie. Po prostu bądźcie zmęczeni tym, jak fft()nie skaluje się tak, jak napisałem (już o tym wspomniałem), że niczego nie udowodniłem za pomocą przechwytywania i że create.fourier.basis()dziwnie skaluje funkcje podstawy.
Taylor
6

Są silnie powiązane. Twój przykład nie jest powtarzalny, ponieważ nie podałeś swoich danych, dlatego zrobię nowy. Po pierwsze, stwórzmy funkcję okresową:

T <- 10
omega <- 2*pi/T
N <- 21
x <- seq(0, T, len = N)
sum_sines_cosines <- function(x, omega){
    sin(omega*x)+2*cos(2*omega*x)+3*sin(4*omega*x)+4*cos(4*omega*x)
}
Yper <- sum_sines_cosines(x, omega)
Yper[N]-Yper[1] # numerically 0

x2 <- seq(0, T, len = 1000)
Yper2 <- sum_sines_cosines(x2, omega)
plot(x2, Yper2, col = "red", type = "l", xlab = "x", ylab = "Y")
points(x, Yper)

wprowadź opis zdjęcia tutaj

Teraz utwórzmy podstawę Fouriera dla regresji. Zauważ, że przy tak naprawdę nie ma sensu tworzyć więcej niż funkcji bazowych, tj. niesynchroniczne sinus i cosinus, ponieważ składniki o wyższej częstotliwości są aliasowane na takiej siatce. Na przykład sinus o częstotliwości jest nie do odróżnienia od costant (sinus): rozważ przypadek , tj. . W każdym razie, jeśli chcesz dwukrotnie sprawdzić, po prostu przejdź do fragmentu kodu poniżej i spójrz na dwie ostatnie kolumny: zobaczysz, że są one faktycznie bezużyteczne (i powodują problemy z dopasowaniem, ponieważ macierz projektowa jest teraz pojedyncza ).N=2k+1N2N3=2(k1)kωN=3k=1N-2N

# Fourier Regression with fda
library(fda)
mybasis <- create.fourier.basis(c(0,T),N-2)
basisMat <- eval.basis(x, mybasis)
FDA_regression <- lm(Yper ~ basisMat-1)
FDA_coef <-coef(FDA_regression)
barplot(FDA_coef)

wprowadź opis zdjęcia tutaj

Zauważ, że częstotliwości są dokładnie właściwe, ale amplitudy składników niezerowych nie są (1,2,3,4). Powodem jest to, że fdafunkcje podstawowe Fouriera są skalowane w dziwny sposób: ich maksymalna wartość nie wynosi 1, tak jak w przypadku zwykłej podstawy Fouriera . Nie jest to również , jak byłoby w przypadku ortonormalnej podstawy Fouriera, .1,sinωx,cosωx,1π12π,sinωxπ,cosωxπ,

# FDA basis has a weird scaling
max(abs(basisMat))
plot(mybasis)

wprowadź opis zdjęcia tutaj

Wyraźnie widać, że:

  1. maksymalna wartość jest mniejsza niż1π
  2. podstawa Fouriera (skrócona do pierwszych składników ) zawiera stałą funkcję (czarna linia), sinusy o rosnącej częstotliwości (krzywe, które są równe 0 na granicach domen) i cosinus o rosnącej częstotliwości (krzywe, które są równa 1 na granicach domeny), tak jak powinno byćN2

Proste skalowanie podanej przez Fouriera podstawy Fouriera fda, aby uzyskać zwykłą podstawę Fouriera, prowadzi do współczynników regresji o oczekiwanych wartościach:

basisMat <- basisMat/max(abs(basisMat))
FDA_regression <- lm(Yper ~ basisMat-1)
FDA_coef <-coef(FDA_regression)
barplot(FDA_coef, names.arg = colnames(basisMat), main = "rescaled FDA coefficients")

wprowadź opis zdjęcia tutaj

Spróbujmy fftteraz: zwróć uwagę, że ponieważ Yperjest to sekwencja okresowa, ostatni punkt tak naprawdę nie dodaje żadnych informacji (DFT sekwencji jest zawsze okresowa). W ten sposób możemy odrzucić ostatni punkt przy obliczaniu FFT. Ponadto FFT to szybki algorytm numeryczny do obliczania DFT, a DFT ciągu liczb rzeczywistych lub zespolonych jest złożony . Tak więc naprawdę chcemy modułów współczynników FFT:

# FFT
fft_coef <- Mod(fft(Yper[1:(N-1)]))*2/(N-1)

Mnożymy przez , aby uzyskać takie samo skalowanie jak w przypadku podstawy Fouriera . Gdybyśmy nie skalowali, nadal odzyskalibyśmy prawidłowe częstotliwości, ale wszystkie amplitudy byłyby skalowane o ten sam współczynnik w stosunku do tego, co znaleźliśmy wcześniej. Teraz wykreślmy współczynniki fft:2N11,sinωx,cosωx,

fft_coef <- fft_coef[1:((N-1)/2)]
terms <- paste0("exp",seq(0,(N-1)/2-1))
barplot(fft_coef, names.arg = terms, main = "FFT coefficients")

wprowadź opis zdjęcia tutaj

Ok: częstotliwości są poprawne, ale zauważ, że teraz podstawowymi funkcjami nie są już sinus i cosinus (są złożonymi wykładnikami , gdzie oznaczam jednostkę urojoną). Zauważ też, że zamiast zestawu niezerowych częstotliwości (1,2,3,4), jak poprzednio, otrzymaliśmy zestaw (1,2,5). Powodem jest to, że termin w tym złożonym rozwinięciu współczynnika (a więc jest złożony) odpowiada dwóm rzeczywistym warunkom w trygonometryczne rozszerzenie podstawy, ze względu na formułę Eulera . Moduł współczynnika zespolonego jest równy sumie kwadratury dwóch rzeczywistych współczynników, tj.expniωxixnexpniωxxnansin(nωx)+bncos(nωx)expix=cosx+isinx 5=|xn|=an2+bn2 . W rzeczywistości .5=33+42

DeltaIV
źródło
1
dzięki DeltaIV dane są dailydostarczane wraz z fdapakietem.
qoheleth
@ qoheleth Nie wiedziałem. Tego wieczoru zmodyfikuję moją odpowiedź za pomocą waszego zestawu danych i wyjaśnię kilka punktów.
DeltaIV