Dari perspektif statistik: Transformasi Fourier vs regresi dengan basis Fourier

13

Saya mencoba memahami apakah transformasi Fourier diskrit memberikan representasi kurva yang sama dengan regresi menggunakan basis Fourier. Sebagai contoh,

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 memberikan bilangan kompleks, sedangkan regresi memberikan bilangan real.

Apakah mereka menyampaikan informasi yang sama? Apakah ada peta satu ke satu di antara dua set angka?

(Saya akan sangat menghargai jika jawabannya ditulis dari perspektif ahli statistik alih-alih dari perspektif insinyur. Banyak bahan daring yang dapat saya temukan memiliki jargon teknik di mana-mana, yang membuat mereka kurang cocok untuk saya.)

qoheleth
sumber
Saya tidak terbiasa dengan cuplikan kode Anda, jadi tidak dapat mengatakan apakah masalah berikut ini berlaku di sana. Namun, biasanya basis DFT didefinisikan dalam hal frekuensi integral ("bilangan bulat"), sedangkan "basis fourier" umum untuk regresi dapat menggunakan rasio frekuensi sewenang-wenang (misalnya termasuk irasional, setidaknya dalam aritmatika kontinu). Ini mungkin juga menarik.
GeoMatt22
Saya pikir semua orang akan mendapat manfaat jika Anda menulis pertanyaan Anda dalam istilah matematika (berbeda dengan potongan kode). Apa masalah regresi yang Anda selesaikan? Apa fungsi dasar Fourier yang Anda gunakan? Anda akan terkejut melihat bagaimana jawaban untuk pertanyaan Anda akan meningkat.
Yair Daon

Jawaban:

15

Mereka sama. Begini caranya ...

Melakukan Regresi

Katakanlah Anda sesuai dengan model dimana t = 1 , ... , N dan n = lantai ( N / 2 ) . Ini tidak cocok untuk regresi linier, jadi alih-alih Anda menggunakan beberapa trigonometri ( cos ( a + b ) = cos

yt=j=1nAjcos(2πt[j/N]+ϕj)
t=1,,Nn=floor(N/2) ) dan sesuai dengan model yang setara: y t = 2 , j sin ( 2 π t [ j / N ] ) .cos(a+b)=cos(a)cos(b)sin(a)sin(b) Menjalankan regresi linier pada semua frekuensi Fourier{j/N:j=1,,n}memberi Anda banyak (2n) dari beta: { β i,j
yt=j=1nβ1,jcos(2πt[j/N])+β2,jsin(2πt[j/N]).
{j/N:j=1,,n}2n , i = 1 , 2 . Untuk setiap j , jika Anda ingin menghitung pasangan dengan tangan, Anda dapat menggunakan:{β^i,j}i=1,2j

dan β 2,j=Σ N t = 1 yt

β^1,j=t=1Nytcos(2πt[j/N])t=1Ncos2(2πt[j/N])
Ini adalah formula regresi standar.
β^2,j=t=1Nytsin(2πt[j/N])t=1Nsin2(2πt[j/N]).

Melakukan Transformasi Fourier Diskrit

Ketika Anda menjalankan transformasi Fourier, Anda menghitung, untuk :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])).

Ini adalah angka kompleks (perhatikan ). Untuk melihat mengapa persamaan itu berlaku, perlu diingat bahwa e i x = cos ( x ) + i sin ( x ) , cos ( - x ) = cos ( x ) dan sin ( - x ) = - sin ( x ) .ieix=cos(x)+isin(x)cos(x)=cos(x)sin(x)=sin(x)

Untuk setiap , mengambil kuadrat dari konjugat kompleks memberi Anda " periodogram :"j

|d(j/N)|2=N1(t=1Nytcos(2πt[j/N]))2+N1(t=1Nytsin(2πt[j/N]))2.
Dalam R, menghitung vektor ini adalah I <- abs(fft(Y))^2/length(Y), yang agak aneh, karena Anda harus skala itu.

Anda juga dapat menentukan " skala periodogram "

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

Koneksi Antara Keduanya

Ternyata hubungan antara regresi dan dua periodogram adalah:

P(j/N)=β^1,j2+β^2,j2.
jt=1Ncos2(2πt[j/N])=t=1Nsin2(2πt[j/N])=N/2

Sumber: https://www.amazon.com/Time-Analysis-Its-Applications-Statistics/dp/144197864X

Taylor
sumber
1
+1 untuk jawaban dan sumbernya. Akan lebih baik jika Anda dapat menunjukkan hasilnya dengan Robjek yang saya posting.
qoheleth
@ qoheleth aku akan menyerahkan itu padamu. Hanya menjadi bosan bagaimana fft()tidak skala cara saya menulis (saya sudah menyebutkan ini), bahwa saya belum membuktikan apa pun dengan penyadapan, dan yang create.fourier.basis()skala dasar berfungsi dengan aneh.
Taylor
6

Mereka sangat terkait. Contoh Anda tidak dapat direproduksi karena Anda tidak memasukkan data Anda, jadi saya akan membuat yang baru. Pertama-tama, mari kita buat fungsi periodik:

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)

enter image description here

Sekarang, mari kita buat basis Fourier untuk regresi. Perhatikan bahwa, denganN=2k+1, tidak masuk akal untuk membuat lebih dari N-2 fungsi dasar, yaitu, N-3=2(k-1)sinus dan cosinus non-konstan, karena komponen frekuensi yang lebih tinggi alias pada kisi-kisi seperti itu. Misalnya, frekuensi sinuskω tidak dapat dibedakan dari costant (sinus): pertimbangkan kasus N=3yaitu k=1. Lagi pula, jika Anda ingin memeriksa ulang, cukup ganti N-2ke Ndalam cuplikan di bawah ini dan lihat dua kolom terakhir: Anda akan melihat bahwa mereka sebenarnya tidak berguna (dan mereka membuat masalah untuk kecocokan, karena matriks desain sekarang tunggal ).

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

enter image description here

Perhatikan bahwa frekuensinya tepat, tetapi amplitudo komponen bukan nol tidak (1,2,3,4). Alasannya adalah bahwa fdafungsi-fungsi basis Fourier diskalakan dengan cara yang aneh: nilai maksimumnya bukan 1, seperti untuk basis Fourier biasa1,dosaωx,cosωx,.... Ini bukan1π baik, karena akan menjadi dasar Fourier ortonormal, 12π,dosaωxπ,cosωxπ,....

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

enter image description here

Anda jelas melihat bahwa:

  1. nilai maksimum kurang dari 1π
  2. dasar Fourier (terpotong ke yang pertama N-2 istilah) berisi fungsi konstan (garis hitam), sinus frekuensi yang meningkat (kurva yang sama dengan 0 pada batas domain) dan cosinus dari peningkatan frekuensi (kurva yang sama dengan 1 pada batas domain), karena seharusnya

Penskalaan basis Fourier yang diberikan oleh fda, sehingga basis Fourier yang biasa diperoleh, mengarah ke koefisien regresi yang memiliki nilai yang diharapkan:

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

enter image description here

Mari kita coba fftsekarang: perhatikan bahwa karena Yperadalah urutan periodik, titik terakhir tidak benar-benar menambahkan informasi (DFT urutan selalu periodik). Dengan demikian kita dapat membuang poin terakhir saat menghitung FFT. Juga, FFT hanyalah algoritma numerik cepat untuk menghitung DFT, dan DFT dari urutan bilangan real atau kompleks adalah kompleks . Jadi, kami benar-benar ingin modulus koefisien FFT:

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

Kami berkembang biak dengan 2N-1 untuk memiliki skala yang sama dengan basis Fourier 1,dosaωx,cosωx,.... Jika kami tidak skala, kami masih akan memulihkan frekuensi yang benar, tetapi amplitudo semua akan diskalakan oleh faktor yang sama sehubungan dengan apa yang kami temukan sebelumnya. Sekarang mari kita plot koefisien fft:

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

enter image description here

Ok: frekuensinya benar, tetapi perhatikan bahwa sekarang fungsi dasar bukan sinus dan cosinus lagi (mereka eksponensial kompleks expnsayaωx, dimana dengan sayaSaya menunjukkan unit imajiner). Perhatikan juga bahwa alih-alih satu set frekuensi bukan nol (1,2,3,4) seperti sebelumnya, kami mendapat satu set (1,2,5). Alasannya adalah istilah ituxnexpnsayaωx dalam ekspansi koefisien kompleks ini (dengan demikian xn kompleks) sesuai dengan dua istilah nyata Sebuahnssayan(nωx)+bncHais(nωx) dalam ekspansi dasar trigonometri, karena rumus Euler expsayax=cosx+sayadosax. Modulus koefisien kompleks sama dengan jumlah dalam kuadrature dari dua koefisien nyata, yaitu,|xn|=Sebuahn2+bn2. Faktanya,5=33+42.

DeltaIV
sumber
1
terima kasih DeltaIV, datanya dailydilengkapi dengan fdapaket.
qoheleth
@ qoheleth aku tidak tahu. Malam ini saya akan mengubah jawaban saya menggunakan dataset Anda, dan saya akan mengklarifikasi beberapa poin.
DeltaIV