rm(list=ls())
sample<-read.table("sample.dat")

ev<-eigen(var(sample)->V)
ev
V

evectors <- NULL
for (i in 1:3)
evectors <- rbind(evectors, rbind(c(0,0,0),ev$vectors[1:3,i]))
evectors


par(mfrow=c(2,2))
plot(sample[,1:2],asp=1,cex=0.5,pch=16)
abline(0,ev$vectors[2,1]/ev$vectors[1,1],col=2)
abline(0,ev$vectors[2,2]/ev$vectors[1,2],col=3)
abline(0,ev$vectors[2,3]/ev$vectors[1,3],col=4)
lines(evectors[1:2,1],evectors[1:2,2],col=5,lwd=2)
lines(evectors[3:4,1],evectors[3:4,2],col=6,lwd=2)
lines(evectors[5:6,1],evectors[5:6,2],col=7,lwd=2)
legend("bottomright",legend=c("EV 1", "EV 2", "EV 3"),col=5:7,lwd=2)

plot(sample[,2:3],asp=1,pch=16,cex=0.5)
abline(0,ev$vectors[3,1]/ev$vectors[2,1],col=2)
abline(0,ev$vectors[3,2]/ev$vectors[2,2],col=3)
abline(0,ev$vectors[3,3]/ev$vectors[2,3],col=4)
lines(evectors[1:2,2],evectors[1:2,3],col=5,lwd=2)
lines(evectors[3:4,2],evectors[3:4,3],col=6,lwd=2)
lines(evectors[5:6,2],evectors[5:6,3],col=7,lwd=2)

plot(sample[,c(1,3)],asp=1,pch=16,cex=0.5)
abline(0,ev$vectors[3,1]/ev$vectors[1,1],col=2)
abline(0,ev$vectors[3,2]/ev$vectors[1,2],col=3)
abline(0,ev$vectors[3,3]/ev$vectors[1,3],col=4)
lines(evectors[1:2,1],evectors[1:2,3],col=5,lwd=2)
lines(evectors[3:4,1],evectors[3:4,3],col=6,lwd=2)
lines(evectors[5:6,1],evectors[5:6,3],col=7,lwd=2)


library(rgl)
plot3d(sample,type="s",size=.75,col="gray")
#plot3d(sample)
lines3d(evectors, col=1, lwd=2)
for (i in 1:3)
abclines3d(0,0,0,evectors[2*i,],col=i+1)

plot3d(ellipse3d(x=V,centre= apply(sample,2,mean)),level=.95,col="orchid2",alpha=0.4,add=TRUE)

# nicht gefragt: 
plot3d(sample,type="s",size=.75,col=c(2:4,rep("gray",197)))
for (i in 1:3)
abclines3d(0,0,0,evectors[2*i,],col=i+1)
par(mfrow=c(2,2))
plot(sample[,1:2],asp=1,pch=16,col=c(2:4,rep("grey",197)))
plot(sample[,2:3],asp=1,pch=16,col=c(2:4,rep("grey",197)))
plot(sample[,c(1,3)],asp=1,pch=16,col=c(2:4,rep("grey",197)))

# 2 Dimensionen:  die erste 2 Hauptkomponenten beschreiben 99% der Änderung im Datensatz
cumsum(ev$values)/sum(ev$values)

pc<-princomp(sample)
summary(pc)

str(pc)
pc$loadings # ist die "Rotationsmatrix"
# also genau die Matrix, die die Eigenvektoren von V enthält

# pc$scores sind Elemente aus "sample" im neuen Koordinatensystem
plot(pc$scores[,1:2],asp=1,pch=16,col=c(2:4,rep("grey",197)))

par(mfrow=c(2,2))
r<-c(-6,6)
plot(pc$scores[,1:2],xlim=r,ylim=r,pch=16,col=c(2:4,rep("grey",197)))
plot(pc$scores[,c(3,2)],xlim=r,ylim=r,pch=16,col=c(2:4,rep("grey",197)))
plot(pc$scores[,c(1,3)],xlim=r,ylim=r,pch=16,col=c(2:4,rep("grey",197)))


##### A2
set.seed(369) #funktioniert gut
# set.seed(444) #weniger, aber gleiche Folgerungen
# set.seed(446) Panik!
t<-seq(-1,1,by=0.01)
length(t)
#n<-length(t)-1
n<-length(t)
m<-(1-t^2)^2
y<-m+rnorm(length(t),0,0.2)
plot(t,y)
lines(t,m)
# b,c und d

K=seq(1,10,by=1)
length(K)->nk
cv<-rep(0,nk)
appr<-rep(0,nk)

#par(mfrow=c(1,3))
#plot(t,y)
#lines(t,m)

mdax<-matrix(0,201,nk)

for (i in 1:nk){
	k<-K[i]
	phi=matrix(0,length(t),k+1)
	for (l in 1:(k+1))
		phi[,l]=cos(pi*(l-1)*t)
	inve<-solve(t(phi)%*%phi)
	S<-(phi)%*%inve%*%t(phi)
	mdach<-S%*%y
	mdax[,i]<-mdach

#	lines(t,mdach,col=rainbow(nk)[i],lwd=1)
	appr[i]<-sum(((y-mdach)/(1-diag(S)))^2)/n

	summanden<-rep(0,n)
	for (j in 1:n){
		Smj<-S[,-j]
		Ymj<-y[-j]
		mdachj<-Smj[j,]%*%Ymj
		summanden[j]<-(y[j]-mdachj)^2
		mm<-mdachj
	}
	cv[i]<-sum(summanden)/n
}
 
#legend("bottom",title="Anzahl cos", legend=1:10,col=rainbow(10),lwd=2,ncol=2)
mi<-min(appr,cv)
M<-max(appr,cv)
par(mfrow=c(1,3))
plot(appr,pch=16,col=rainbow(nk),cex=2,ylim=c(mi,M))
legend("bottomright",title="Anzahl cos", legend=1:10,col=rainbow(10),pch=16,cex==2,ncol=2)
lines(appr)
plot(cv,pch=16,col=rainbow(nk),cex=2,ylim=c(mi,M))
lines(cv)
plot(appr,cv,pch=16,col=rainbow(nk),cex=2)
abline(0,1)