--------------------------------- # Fourierreihen als # # lineares Regressionsproblem # --------------------------------- n = 1000 # devide [-pi,pi] into n or n-1 pieces p = 20 # number of basis functions (sines and cosines) x = seq(from=-pi,to=pi,length=n) # Delta_x = 2*pi/(n-1) # # set up regressors: # X = outer(x,1:p,FUN="*") # n times p matrix mit # Eintraegen x[i]*j Xsin = sin(X) Xcos = cos(X) # that's all to be done, now we can Fourier expand any y with # lm(y ~ Xsin + Xcos), this is super-easy!! # # example 1: y(x)=x, ungerade Funktion # y = x plot(x,y,type="l") res = lm(y ~ Xsin) summary(res) plot(x,res\$fit,type="l") lines(x,y,col="red") plot(res\$coeff) plot(2/res\$coeff[-1]) lines(1:20,col="red") lines((-1):(-20),col="red") res = lm(y ~ Xcos) summary(res) # alle betas 0 (oder e-16) # da ungerade Funktion res = lm(y ~ Xsin + Xcos) summary(res) # let's see how the terms add up: plot(x,y,type="l",col="red") for(pp in 1:p) { X = outer(1:pp,x,FUN="*") X = t(X) Xsin = sin(X) res = lm(y ~ Xsin) lines(x,res\$fit) readline("press enter to continue..") } # die Bilder werden vielleicht uebersichtlicher, # wenn man jeweils nur eine Fourier-Summe plottet: for(pp in 1:p) { X = outer(1:pp,x,FUN="*") X = t(X) Xsin = sin(X) res = lm(y ~ Xsin) info = paste("p =",pp) plot(x,y,type="l",col="red",main=info) lines(x,res\$fit) readline("press enter to continue..") } # # example 2: ungerade Rechteck-Funktion: # y = ifelse(x>=0,1,-1) plot(x,y,type="l") for(pp in 1:p) { X = outer(1:pp,x,FUN="*") X = t(X) Xsin = sin(X) res = lm(y ~ Xsin) info = paste("p =",pp) plot(x,y,type="l",col="red",ylim=c(-1.5,1.5),main=info) lines(x,res\$fit) readline("press enter to continue..") } # # example 3: y=|x|, gerade Funktion # y = abs(x) plot(x,y,type="l") res = lm(y ~ Xcos) summary(res) plot(x,res\$fit,type="l") lines(x,y,col="red") # nearly perfect fit plot(res\$coeff) # let's see how the terms add up: for(pp in 1:p) { X = outer(1:pp,x,FUN="*") X = t(X) Xcos = cos(X) res = lm(y ~ Xcos) info = paste("p =",pp) plot(x,y,type="l",col="red",main=info) lines(x,res\$fit) readline("press enter to continue..") } # # example 4: y=exp(x), weder gerade noch ungerade # -> we need sines and cosines # y = exp(x) plot(x,y,type="l",col="red") res = lm(y ~ Xsin + Xcos) summary(res) lines(x,res\$fit,type="l") plot(res\$coeff) # plottet erst die Konstante, # dann die sin-Koeffizienten, # dann die cos-Koeffizienten # let's see how the terms add up: for(pp in 1:p) { X = outer(1:pp,x,FUN="*") X = t(X) Xsin = sin(X) Xcos = cos(X) res = lm(y ~ Xsin + Xcos) info = paste("p =",pp) plot(x,y,type="l",col="red",ylim=c(-5,25),main=info) lines(x,res\$fit) readline("press enter to continue..") } # alles in einem Bild: plot(x,y,type="l",col="red",ylim=c(-5,25)) for(pp in 1:p) { X = outer(1:pp,x,FUN="*") X = t(X) Xsin = sin(X) Xcos = cos(X) res = lm(y ~ Xsin + Xcos) lines(x,res\$fit) } # die letzte Fourier-Summe in gruen und etwas # fetter zum bestehenden plot hinzufuegen: points(x,res\$fit, pch=".", cex=4, col="green") # pch = point character # cex = character expansion, scheint wie folgt # zu wirken: # geplotteteGroesse = cex * defaultGroesse --------------------------------------------- # Noch ein paar Infos zu den Parametern # # pch und cex # --------------------------------------------- plot(1:10, pch=".") plot(1:10, pch=".", cex=5) plot(1:10, pch=".", cex=10) # "." scheint offensichtlich ein Quadrat zu # sein, kein Kreis. Schauen wir uns mal einige # plot characters an: for(i in 1:30) { info = paste("pch =",i) plot(1:10,pch=i,main=info) readline("press enter to continue..") } # also: pch=20 sind kleine ausgefuellte Kreise, # pch=21 sind kleine unausgefuellte Kreise, # kann man die mit cex beliebig vergroessern # und verkleinern? # vergroessern: for(scale in 1:20) { info = paste("cex =",scale) plot(1:10,pch=20,cex=scale,main=info) readline("press enter to continue..") } # verkleinern: for(scale in 1/(1:10) ) { info = paste("cex =",scale) plot(1:10,pch=20,cex=scale,main=info) readline("press enter to continue..") } # fuer den default plot character: # vergroessern: for(scale in 1:20) { info = paste("cex =",scale) plot(1:10,cex=scale,main=info) readline("press enter to continue..") } # verkleinern: for(scale in 1/(1:10) ) { info = paste("cex =",scale) plot(1:10,cex=scale,main=info) readline("press enter to continue..") }