#1234567890123456789012345678901234567890123456789012345678901234567890123456789
#Ruth Meili

# Rg2e; Aufgabe 1 von Serie 2

# Vordefinitionen:
a<-0
b<-7
n1<-50
n2<-200
sd1<-0.2
sd2<-0.5

set.seed(155)

# Lösung für a):
f.reg<-function(x){dnorm(x,3,0.5)+dnorm(x,5,0.2)}
x.n1<-seq(a,b,len=n1)
Y1.n1<-f.reg(x.n1)
Y2.n1<-f.reg(x.n1)+rnorm(n1,mean=0,sd=sd1)
Y3.n1<-f.reg(x.n1)+rnorm(n1,mean=0,sd=sd2)

x.n2<-seq(a,b,len=n2)
Y1.n2<-f.reg(x.n2)
Y2.n2<-f.reg(x.n2)+rnorm(n2,mean=0,sd=sd1)
Y3.n2<-f.reg(x.n2)+rnorm(n2,mean=0,sd=sd2)

par(mfrow=c(2,3))
plot(x.n1,Y1.n1,main="kein Rauschen")
plot(x.n1,Y2.n1,main="Rauschen mit sd=0.2")
plot(x.n1,Y3.n1,main="Rauschen mit sd=0.5")
plot(x.n2,Y1.n2,main="kein Rauschen")
plot(x.n2,Y2.n2,main="Rauschen mit sd=0.2")
plot(x.n2,Y3.n2,main="Rauschen mit sd=0.5")

# Lösung für b)
library(lokern)
h<-c(.2,.5,1)
grid<-seq(a,b,0.01)
scale<-(b-a)/(length(grid)-1)

fitnp.Y1.n1.h1<-glkerns(x.n1,Y1.n1,x.out=grid,bandwidth=h[1])$est
fitnp.Y2.n1.h1<-glkerns(x.n1,Y2.n1,x.out=grid,bandwidth=h[1])$est
fitnp.Y3.n1.h1<-glkerns(x.n1,Y3.n1,x.out=grid,bandwidth=h[1])$est
fitnp.Y1.n2.h1<-glkerns(x.n2,Y1.n2,x.out=grid,bandwidth=h[1])$est
fitnp.Y2.n2.h1<-glkerns(x.n2,Y2.n2,x.out=grid,bandwidth=h[1])$est
fitnp.Y3.n2.h1<-glkerns(x.n2,Y3.n2,x.out=grid,bandwidth=h[1])$est

fitnp.Y1.n1.h2<-glkerns(x.n1,Y1.n1,x.out=grid,bandwidth=h[2])$est
fitnp.Y2.n1.h2<-glkerns(x.n1,Y2.n1,x.out=grid,bandwidth=h[2])$est
fitnp.Y3.n1.h2<-glkerns(x.n1,Y3.n1,x.out=grid,bandwidth=h[2])$est
fitnp.Y1.n2.h2<-glkerns(x.n2,Y1.n2,x.out=grid,bandwidth=h[2])$est
fitnp.Y2.n2.h2<-glkerns(x.n2,Y2.n2,x.out=grid,bandwidth=h[2])$est
fitnp.Y3.n2.h2<-glkerns(x.n2,Y3.n2,x.out=grid,bandwidth=h[2])$est

fitnp.Y1.n1.h3<-glkerns(x.n1,Y1.n1,x.out=grid,bandwidth=h[3])$est
fitnp.Y2.n1.h3<-glkerns(x.n1,Y2.n1,x.out=grid,bandwidth=h[3])$est
fitnp.Y3.n1.h3<-glkerns(x.n1,Y3.n1,x.out=grid,bandwidth=h[3])$est
fitnp.Y1.n2.h3<-glkerns(x.n2,Y1.n2,x.out=grid,bandwidth=h[3])$est
fitnp.Y2.n2.h3<-glkerns(x.n2,Y2.n2,x.out=grid,bandwidth=h[3])$est
fitnp.Y3.n2.h3<-glkerns(x.n2,Y3.n2,x.out=grid,bandwidth=h[3])$est

isenp.Y1.n1.h1<-sum((fitnp.Y1.n1.h1-f.reg(grid))^2)*scale
isenp.Y2.n1.h1<-sum((fitnp.Y2.n1.h1-f.reg(grid))^2)*scale
isenp.Y3.n1.h1<-sum((fitnp.Y3.n1.h1-f.reg(grid))^2)*scale
isenp.Y1.n2.h1<-sum((fitnp.Y1.n2.h1-f.reg(grid))^2)*scale
isenp.Y2.n2.h1<-sum((fitnp.Y2.n2.h1-f.reg(grid))^2)*scale
isenp.Y3.n2.h1<-sum((fitnp.Y3.n2.h1-f.reg(grid))^2)*scale

isenp.Y1.n1.h2<-sum((fitnp.Y1.n1.h2-f.reg(grid))^2)*scale
isenp.Y2.n1.h2<-sum((fitnp.Y2.n1.h2-f.reg(grid))^2)*scale
isenp.Y3.n1.h2<-sum((fitnp.Y3.n1.h2-f.reg(grid))^2)*scale
isenp.Y1.n2.h2<-sum((fitnp.Y1.n2.h2-f.reg(grid))^2)*scale
isenp.Y2.n2.h2<-sum((fitnp.Y2.n2.h2-f.reg(grid))^2)*scale
isenp.Y3.n2.h2<-sum((fitnp.Y3.n2.h2-f.reg(grid))^2)*scale

isenp.Y1.n1.h3<-sum((fitnp.Y1.n1.h3-f.reg(grid))^2)*scale
isenp.Y2.n1.h3<-sum((fitnp.Y2.n1.h3-f.reg(grid))^2)*scale
isenp.Y3.n1.h3<-sum((fitnp.Y3.n1.h3-f.reg(grid))^2)*scale
isenp.Y1.n2.h3<-sum((fitnp.Y1.n2.h3-f.reg(grid))^2)*scale
isenp.Y2.n2.h3<-sum((fitnp.Y2.n2.h3-f.reg(grid))^2)*scale
isenp.Y3.n2.h3<-sum((fitnp.Y3.n2.h3-f.reg(grid))^2)*scale

cat("-------------------------------------------------","\n")
cat("ISE for Nonparametric Fit with h=",h[1],"\n")
cat("-------------------------------------------------","\n")
cat("n=",n1,"sd=",0,":",round(isenp.Y1.n1.h1,4),"\n")
cat("n=",n1,"sd=",sd1,":",round(isenp.Y2.n1.h1,4),"\n")
cat("n=",n1,"sd=",sd2,":",round(isenp.Y3.n1.h1,4),"\n")
cat("n=",n2,"sd=",0,":",round(isenp.Y1.n2.h1,4),"\n")
cat("n=",n2,"sd=",sd1,":",round(isenp.Y2.n2.h1,4),"\n")
cat("n=",n2,"sd=",sd2,":",round(isenp.Y3.n2.h1,4),"\n")
cat("-------------------------------------------------","\n")
cat("ISE for Nonparametric Fit with h=",h[2],"\n")
cat("-------------------------------------------------","\n")
cat("n=",n1,"sd=",0,":",round(isenp.Y1.n1.h2,4),"\n")
cat("n=",n1,"sd=",sd1,":",round(isenp.Y2.n1.h2,4),"\n")
cat("n=",n1,"sd=",sd2,":",round(isenp.Y3.n1.h2,4),"\n")
cat("n=",n2,"sd=",0,":",round(isenp.Y1.n2.h2,4),"\n")
cat("n=",n2,"sd=",sd1,":",round(isenp.Y2.n2.h2,4),"\n")
cat("n=",n2,"sd=",sd2,":",round(isenp.Y3.n2.h2,4),"\n")
cat("-------------------------------------------------","\n")
cat("ISE for Nonparametric Fit with h=",h[3],"\n")
cat("-------------------------------------------------","\n")
cat("n=",n1,"sd=",0,":",round(isenp.Y1.n1.h3,4),"\n")
cat("n=",n1,"sd=",sd1,":",round(isenp.Y2.n1.h3,4),"\n")
cat("n=",n1,"sd=",sd2,":",round(isenp.Y3.n1.h3,4),"\n")
cat("n=",n2,"sd=",0,":",round(isenp.Y1.n2.h3,4),"\n")
cat("n=",n2,"sd=",sd1,":",round(isenp.Y2.n2.h3,4),"\n")
cat("n=",n2,"sd=",sd2,":",round(isenp.Y3.n2.h3,4),"\n")
cat("-------------------------------------------------","\n")

# Lösung für c)

fitnp.Y1.n1.hopt<-glkerns(x.n1,Y1.n1,x.out=grid)$est
fitnp.Y2.n1.hopt<-glkerns(x.n1,Y2.n1,x.out=grid)$est
fitnp.Y3.n1.hopt<-glkerns(x.n1,Y3.n1,x.out=grid)$est
fitnp.Y1.n2.hopt<-glkerns(x.n2,Y1.n2,x.out=grid)$est
fitnp.Y2.n2.hopt<-glkerns(x.n2,Y2.n2,x.out=grid)$est
fitnp.Y3.n2.hopt<-glkerns(x.n2,Y3.n2,x.out=grid)$est

isenp.Y1.n1.hopt<-sum((fitnp.Y1.n1.hopt-f.reg(grid))^2)*scale
isenp.Y2.n1.hopt<-sum((fitnp.Y2.n1.hopt-f.reg(grid))^2)*scale
isenp.Y3.n1.hopt<-sum((fitnp.Y3.n1.hopt-f.reg(grid))^2)*scale
isenp.Y1.n2.hopt<-sum((fitnp.Y1.n2.hopt-f.reg(grid))^2)*scale
isenp.Y2.n2.hopt<-sum((fitnp.Y2.n2.hopt-f.reg(grid))^2)*scale
isenp.Y3.n2.hopt<-sum((fitnp.Y3.n2.hopt-f.reg(grid))^2)*scale

cat("ISE for Nonparametric Fit with h=hopt","\n")
cat("-------------------------------------------------","\n")
cat("n=",n1,"sd=",0,":",round(isenp.Y1.n1.hopt,4),"\n")
cat("n=",n1,"sd=",sd1,":",round(isenp.Y2.n1.hopt,4),"\n")
cat("n=",n1,"sd=",sd2,":",round(isenp.Y3.n1.hopt,4),"\n")
cat("n=",n2,"sd=",0,":",round(isenp.Y1.n2.hopt,4),"\n")
cat("n=",n2,"sd=",sd1,":",round(isenp.Y2.n2.hopt,4),"\n")
cat("n=",n2,"sd=",sd2,":",round(isenp.Y3.n2.hopt,4),"\n")
cat("-------------------------------------------------","\n")

# Grafische Darstellung von a)-c):

ymin<-min(Y1.n1,Y2.n1,Y3.n1,Y1.n2,Y2.n2,Y3.n2)
ymax<-max(Y1.n1,Y2.n1,Y3.n1,Y1.n2,Y2.n2,Y3.n2)

par(mfrow=c(3,2))
plot(x.n1,Y1.n1,ylim=c(ymin,ymax))
lines(grid,f.reg(grid),col=1)
lines(grid,fitnp.Y1.n1.h1,col=2)
lines(grid,fitnp.Y1.n1.h2,col=3)
lines(grid,fitnp.Y1.n1.h3,col=4)
lines(grid,fitnp.Y1.n1.hopt,col=5)
 title("kleine Stichprobe, kein Rauschen")
 legend(0,2.5,legend=c("wahre Kurve","Bandbreite 0.2","Bandbreite 0.5","Bandbreite1", "optimierte Bandbreite"),col = c(1,2,3,4,5))

plot(x.n2,Y1.n2,ylim=c(ymin,ymax))
lines(grid,f.reg(grid),col=1)
lines(grid,fitnp.Y1.n2.h1,col=2)
lines(grid,fitnp.Y1.n2.h2,col=3)
lines(grid,fitnp.Y1.n2.h3,col=4)
lines(grid,fitnp.Y1.n2.hopt,col=5)

plot(x.n1,Y2.n1,ylim=c(ymin,ymax))
lines(grid,f.reg(grid),col=1)
lines(grid,fitnp.Y2.n1.h1,col=2)
lines(grid,fitnp.Y2.n1.h2,col=3)
lines(grid,fitnp.Y2.n1.h3,col=4)
lines(grid,fitnp.Y2.n1.hopt,col=5)

plot(x.n2,Y2.n2,ylim=c(ymin,ymax))
lines(grid,f.reg(grid),col=1)
lines(grid,fitnp.Y2.n2.h1,col=2)
lines(grid,fitnp.Y2.n2.h2,col=3)
lines(grid,fitnp.Y2.n2.h3,col=4)
lines(grid,fitnp.Y2.n2.hopt,col=5)

plot(x.n1,Y3.n1,ylim=c(ymin,ymax))
lines(grid,f.reg(grid),col=1)
lines(grid,fitnp.Y3.n1.h1,col=2)
lines(grid,fitnp.Y3.n1.h2,col=3)
lines(grid,fitnp.Y3.n1.h3,col=4)
lines(grid,fitnp.Y3.n1.hopt,col=5)

plot(x.n2,Y3.n2,ylim=c(ymin,ymax))
lines(grid,f.reg(grid),col=1)
lines(grid,fitnp.Y3.n2.h1,col=2)
lines(grid,fitnp.Y3.n2.h2,col=3)
lines(grid,fitnp.Y3.n2.h3,col=4)
lines(grid,fitnp.Y3.n1.hopt,col=5)

# Lösung fuer d)

N<-100
ise.h1<-rep(0,N)
ise.h2<-rep(0,N)
ise.h3<-rep(0,N)
fit.h1<-matrix(0,nrow=N,ncol=length(grid))
fit.h2<-matrix(0,nrow=N,ncol=length(grid))
fit.h3<-matrix(0,nrow=N,ncol=length(grid))

for(i in 1:N){
Y2.n1<-f.reg(x.n1)+rnorm(n1,mean=0,sd=sd1)
fit.h1[i,]<-glkerns(x.n1,Y2.n1,x.out=grid,bandwidth=h[1])$est
fit.h2[i,]<-glkerns(x.n1,Y2.n1,x.out=grid,bandwidth=h[2])$est
fit.h3[i,]<-glkerns(x.n1,Y2.n1,x.out=grid,bandwidth=h[3])$est
ise.h1[i]<-sum((fit.h1[i,]-f.reg(grid))^2)*scale
ise.h2[i]<-sum((fit.h2[i,]-f.reg(grid))^2)*scale
ise.h3[i]<-sum((fit.h3[i,]-f.reg(grid))^2)*scale
}

meanfit.h1<-apply(fit.h1,2,"mean")
meanfit.h2<-apply(fit.h2,2,"mean")
meanfit.h3<-apply(fit.h3,2,"mean")

dev2fit.h1<-(fit.h1-matrix(rep(meanfit.h1,N),byrow=T,nrow=N))^2
dev2fit.h2<-(fit.h2-matrix(rep(meanfit.h2,N),byrow=T,nrow=N))^2
dev2fit.h3<-(fit.h3-matrix(rep(meanfit.h3,N),byrow=T,nrow=N))^2

bias.h1<-sum((meanfit.h1-f.reg(grid))^2)*scale
bias.h2<-sum((meanfit.h2-f.reg(grid))^2)*scale
bias.h3<-sum((meanfit.h3-f.reg(grid))^2)*scale

var.h1<-sum(dev2fit.h1)*scale/N
var.h2<-sum(dev2fit.h2)*scale/N
var.h3<-sum(dev2fit.h3)*scale/N

mise.h1<-mean(ise.h1)
mise.h2<-mean(ise.h2)
mise.h3<-mean(ise.h3)

cat("Integrated Squared Bias","\n")
cat("-------------------------------------------------","\n")
cat("h=",h[1],":",round(bias.h1,4),"\n")
cat("h=",h[2],":",round(bias.h2,4),"\n")
cat("h=",h[3],":",round(bias.h3,4),"\n")
cat("-------------------------------------------------","\n")
cat("Integrated Variance","\n")
cat("-------------------------------------------------","\n")
cat("h=",h[1],":",round(var.h1,4),"\n")
cat("h=",h[2],":",round(var.h2,4),"\n")
cat("h=",h[3],":",round(var.h3,4),"\n")
cat("-------------------------------------------------","\n")
cat("MISE","\n")
cat("-------------------------------------------------","\n")
cat("h=",h[1],":",round(mise.h1,4),"\n")
cat("h=",h[2],":",round(mise.h2,4),"\n")
cat("h=",h[3],":",round(mise.h3,4),"\n")
cat("-------------------------------------------------","\n")

# Lösung fuer e)

# Wiederholung von a)-d) mit folgenden Parametern für die Regression-Funktion:
f.reg<-function(x){exp(-x)}
a<-0
b<-7

# Oder mit
f.reg<-function(x){exp(x)/(1+exp(x))}
a<--7
b<-7


# Teilaufgabe f)
source("ftp://stat.ethz.ch/NDK/Source-NDK-7/rg2e.R")

f.compare.bdw()


#1234567890123456789012345678901234567890123456789012345678901234567890123456789
#Ruth Meili

# Rg2e; Aufgabe 2 von Serie 2

library(lokern)

source("ftp://stat.ethz.ch/NDK/Source-NDK-7/rg2e.R")
f.compare.veloc(nboys=10,ngirls=10)
# This program plots distance curves and velocity curves of nboys boys and
# ngirls girls (randomly selected from a population of 120 boys and 112 girls)

# The program also allows to extract features characterizing the
# puberty from the velocity curves



