library(fda) rm(list=ls()) str(CanadianWeather) str(daily) tag<-(1:365)-0.5 ort<-daily$place #b base<-create.fourier.basis(c(0,365),nbasis=101) base apply(daily$tempav,1,mean)->mittel Data2fd(tag,daily$tempav,base)->tempfdmm plot(tempfdmm) Data2fd(tag,mittel,base)->mitttfd lines(mitttfd,col=1,lwd=2) dim(daily$tempav) Mit<-matrix(rep(mittel,35),ncol=35) mtav<-daily$tempav-Mit Data2fd(tag,mtav,base)->tempfd plot(tempfd) #c tempfda<-pca.fd(tempfd,nharm=20) summary(tempfda) cumsum(tempfda$varprop) par(mfrow=c(1,3)) plot(tempfda$harmonics[1]->h1,ylim=c(-.1,.1)) plot(tempfda$harmonics[2]->h2,ylim=c(-.1,.1)) plot(tempfda$harmonics[3]->h3,ylim=c(-.1,.1)) par(mfrow=c(1,3)) plot(c(0,365),range(c(-40,20)),type="n") for (i in 1:35) lines(mitttfd+tempfda$scores[i,1]*h1,col=rainbow(35)[i]) plot(c(0,365),range(c(-40,20)),type="n") for (i in 1:35) lines(mitttfd+tempfda$scores[i,1]*h1+tempfda$scores[i,2]*h2,col=rainbow(35)[i]) plot(tempfda$scores[,1:2],pch=16,col=rainbow(35),cex=1.5) # d lok<-rep(0,35) for (i in 1:35){ if(CanadianWeather$region[i]=="Atlantic") lok[i]=1; if(CanadianWeather$region[i]=="Pacific") lok[i]=2; if(CanadianWeather$region[i]=="Continental") lok[i]=3; if(CanadianWeather$region[i]=="Arctic") lok[i]=4; } lok par(mfrow=c(1,3)) plot(c(0,365),range(c(-40,20)),type="n") abline(v=c(80,80+90,80+90+92,365-10),h=0,col=5) for (i in 1:35) lines(mitttfd+tempfda$scores[i,1]*h1,col=lok[i]) lines(mitttfd,col=6,lwd=2) plot(c(0,365),range(c(-40,20)),type="n") abline(v=c(80,80+90,80+90+92,365-10),h=0,col=5) for (i in 1:35) lines(mitttfd+tempfda$scores[i,1]*h1+tempfda$scores[i,2]*h2,col=lok[i]) lines(mitttfd,col=6,lwd=2) plot(tempfda$scores[,1:2],pch=16,col=lok,cex=1.5) library(rgl) plot3d(tempfda$scores[,1:3],col=lok,type="s")