######################################################
## Aufgabe 1
#####################################################
library(lokern)
## Funktionen
f.m <- function(x) 2*tanh(8*(x-1)) + 2*tanh(2*(x-5)) - tanh(4*x-10) - tanh(4*(x-6.6))
c.sample <- function(n, truef, xmin=0, xmax=7, sd) {
    x <- seq(xmin, xmax, len=n)
    y <- truef(x) + rnorm(n, 0, sd=sd)
    list(x=x,y=y)
}

## "Abkürzungen", um den Code unten übersichtlicher zu halten:
p.sd <- function(s) paste0("sd=",s)
p.n  <- function(n) paste0("n=", n)
p.h  <- function(h) paste0("h=", h)

## Globale Werte
du <- 0.01
uu <- seq(0, 7, by = du)
sd.s <- c(0, 0.6, 1.2)
n.s <- c(50, 100)

## Stichproben erstellen - Trick: Konstruiere und arbeite mit einer Matrix von Listen:
Samp <- vector("list", length(sd.s)*length(n.s))
dim(Samp) <- c(length(sd.s), length(n.s))
dimnames(Samp) <- list(p.sd(sd.s), p.n(n.s))
set.seed(144)
for(n in n.s)
  for(s in sd.s)
    Samp[[p.sd(s), p.n(n)]] <- c.sample(n=n, truef=f.m, sd = s)

## a) Bias an den Stichprobenwerten berechnen - Ohne Array
#########################################################   ==========
## Lineare Regression:
dn1 <- Samp[["sd=0", "n=50"]]; x.n1 <- dn1$x
dn2 <- Samp[["sd=0","n=100"]]; x.n2 <- dn2$x
E1 <- predict(lm(y ~ x, data=dn1))
E2 <- predict(lm(y ~ x, data=dn2))
Bias1.n1 <- E1 - f.m(x.n1)
Bias1.n2 <- E2 - f.m(x.n2)
## Regression 5.Grades:
E1 <- predict(lm(y ~ poly(x,5), data= ....))
E2 <- predict(lm(y ~ poly(x,5), data= ....))
Bias5.n1 <- E1 - f.m(x.n1)
Bias5.n2 <- E2 - f.m(x.n2)
## Nichtparametrische Regression -- Aufgepasst, die S - Matrix hängt von
## der Bandbreite h ab, die von glkerns() automatisch gewählt wird;
## aber dies geht nicht gut für Daten ohne Rauschen.  Daher:
(h1 <- with(Samp[["sd=0.2", "n=50"]], glkerns(x,y, x.out=x)) $ bandwidth)
(h2 <- with(Samp[["sd=0.2","n=100"]], glkerns(x,y, x.out=x)) $ bandwidth)
## Anwendung *dieser* S_h - Matrix :
E1 <- with(dn1, glkerns(x,y, x.out=x, bandwidth=h1))$est
E2 <- with(dn2, glkerns(x,y, x.out=x, bandwidth=h2))$est
Biasnp.n1 <- E1  -  f.m(x.n1)
Biasnp.n2 <- E2  -  f.m(x.n2)
## Graphik
par(mfrow=c(1,1)) # .. "nur im Fall .."
plot (x.n1, Bias1.n1, type="l")
lines(x.n2, Bias1.n2,  lty=2)
lines(x.n1, Bias5.n1,  lty=3)
lines(x.n2, Bias5.n2,  lty=4)
lines(x.n1, Biasnp.n1, lty=5)
lines(x.n2, Biasnp.n2, lty=6)
regNms <- c("lin.Reg", "5.gr.Reg", "NonP.Reg")
legend("bottomleft", ## outer(x,y) : "jeder x mit jedem y"
       legend= t(outer(regNms, p.n(c(50,100)), paste, sep=",")),
       lty=1:6, col=rep(1:2, 3), title= ..., bty="n")

## b) Bias an den Grid-Punkten berechnen - mit Array
##################################################
## Predictions
h.s <- c("n=50" = .263,  "n=100" = .219) # ~= (h1, h2)  oben
Pred <- array(dim = c(length(uu), 3, length(sd.s), length(n.s)),
              dimnames = c(list(NULL, regNms), dimnames(Samp)))
for(n in n.s) {
  n. <- p.n(n)
  for(s in sd.s) {
    s. <- p.sd(s)
    Pred[, "lin.Reg", s., n.] <-
        predict(lm(y ~ x, data= Samp[[s., n.]]),
                newdata = data.frame(x=uu))
    Pred[, "kub.Reg", s., n.] <-
        predict(lm(y ~ poly(x,3), data= Samp[[s., n.]]),
                newdata = data.frame(x=uu))
    Pred[, "NonP.Reg", s., n.] <-
        with(Samp[[s., n.]],
             glkerns(x,y, x.out=uu, bandwidth = h.s[[n.]]))$est
  }
}
## Bias
Bias <- array(dim = c(length(uu), 3, l.n),
              dimnames = c(list(NULL, regNms), dimnames(Samp)[2]))
for(n in n.s) {
  n. <- p.n(n)
  for (i in 1:3)
     Bias[, i, n.] <- Pred[, i, p.sd(0), n.]  -  f.m(uu)
}
yl <- range(Bias)
## Die Graphiken dazu sind:

matplot (uu,Bias[,, "n=50"], type="l",lty=1:3,ylim=yl,col=1:3, ylab="Bias")
matlines(uu,Bias[,,"n=100"],          lty=4:6,lwd=2,  col=1:3)
legend("bottomleft", lty=1:6, col=c(1:3, 1:3), lwd=c(1,1,1, 2,2,2), bty="n",
       legend=c(paste(regNms,"n= 50", sep=","),
                paste(regNms,"n=100", sep=",")))

## c) ISB berechnen - mit Arrays
##############################
ISB <- array(dim = c(3, l.n),
             dimnames = list(regNms, dimnames(Samp)[[2]]))
for(n in n.s) {
  n. <- p.n(n)
  for (i in 1:3)
      ISB[i, n.] <- sum(Bias[, i, n.]^2)*du
}
ISB

#################################################################################


########################################################
## Aufgabe 2
########################################################
library(lokern)
## Funktionen, globale Werte, Abkürzungen [p.sd(),..] :
## ----> siehe oben, Aufgabe 1
##       =====================

l.h <- length(hh <- seq(0.1, 1, by=0.1))
## a) Stichproben erstellen - mit array
#####################################
set.seed(155)
l.s <- length(sd.s <- c(0, 0.2, 0.5))
l.n <- length(n.s <- c(40, 200))  ## andere n's als Aufgabe 1 !
Samp <- vector("list", l.s*l.n)
dim(Samp) <- c(l.s, l.n)
dimnames(Samp) <- list(p.sd(sd.s), p.n(n.s))
for(n in n.s)
  for(s in sd.s)
    Samp[[....]] <- c.sample(....)

## b) --------
## Gefittete Werte der nichtparametrischen Regression - variable f.n
## ISE - variable ise.n
## Bestimmung der Bandbreite des minimalen ISE - variable h.iseopt.n
f.n <- array(dim = c(length(uu), l.h, l.s, l.n),
             dimnames = c(list(NULL),list(p.h(hh)), dimnames(Samp)))
ise.n <- array(dim = c(l.h, l.s, l.n),
               dimnames = c(list(p.h(hh)),dimnames(Samp)))
h.ise.n <- h.iseopt.n <-
    array(dim = c(l.s, l.n), dimnames = dimnames(Samp))
for(n in n.s) {
  n. <- p.n(n)
  for(s in sd.s) {
    for(h in hh) {
      f.n[, p.h(h),p.sd(s), n.] <-
        with(Samp[[....]], glkerns(x,y, x.out = .... ,bandwidth = h))$est
      ise.n[p.h(h),p.sd(s), n.] <-
        sum((f.n[, .... , .... , n.]-f.m(....))^2)* du
    }
    boolean <- (ise.n[,p.sd(s), n.] == min(ise.n[, p.sd(s), n.]))
    h.iseopt.n[p.sd(s), n.] <- hin <- hh[boolean]
    h.ise.n   [p.sd(s), n.] <- ise.n[p.h(hin), p.sd(s), n.]
  }
}
yl <- range(ise.n)

matplot (hh, ise.n[,, "n=40"],type="l",lty=1:3, col=1:3, ylim=yl, lwd=2, ylab="ISE")
matlines(hh, ise.n[,,"n=200"],         lty=4:6, col=4:6,          lwd=2)
f.sds <- format(p.sd(sd.s))
legend("topright", lwd=2, lty=1:6, col=1:6,
       legend=c(paste(f.sds, "n= 40", sep=","),
                paste(f.sds, "n=200", sep=",")))
points( h....., h...., col=2)


## c) ISE Berechnung - mit arrays
###############################
f.opt.n <- array(dim = c(length(uu),l.s, l.n),
                 dimnames = c(list(NULL),dimnames(Samp)))
ise.opt.n <-
  bw.n <- array(dim = c(l.s, l.n), dimnames = dimnames(Samp))
for(n in n.s) {
   n. <- p.n(n)
   for(s in sd.s) {
      f.opt.n[, .... , .... ] <-
        with(.... ,glkerns(x,y, x.out = uu))$est
      bw.n[.... , ....] <-
        with(Samp[[...., ....]],glkerns(x,y, x.out = ....))$bandwidth
      ise.opt.n[p.sd(s), n.] <- sum((f.opt.n[,p.sd(s), n.]-f.m(....))^2)* du
   }
}

matplot(hh,  ise.n[,, "n=40"], type = "l", lty = 1:3, col = 1:3,
        ylab = "ISE", xlim = c(0,1.1), ylim = yl)
matlines(hh, ise.n[,,"n=200"], type = "l", lty = 4:6, col = 4:6)
.....  ## legend; optimality points

## d)
n.sim <- 100 # = N
l.s <- length(sd.s <- c(0, 0.6))
l.n <- length(n.s <- c(40))
du <- 0.01
l.h <- length(hh <- seq(0.1, 1, by=0.1))
Samp <- vector("list", l.s*l.n)
dim(Samp) <- c(l.s, l.n)
dimnames(Samp) <- list(p.sd(sd.s), p.n(n.s))
## ISB berechnen -
                 #Datensatz genieren
for(n in n.s)
  for(s in sd.s)
    Samp[[....]] <- c.sample(....)

f.n <- array(dim = c(....),
             dimnames = c(list(NULL, p.h(hh)), dimnames(Samp)))
for(n in n.s) {
   n. <- p.n(n)
   for(s in sd.s) {
      for(h in hh) {
         f.n[, .... , .... , ....] <- ....
       }
   }
}
## ISB berechnen
isb.n <- array(dim = c(l.h, l.s, l.n),dimnames = c(.... , ....))
for(n in n.s) {
   n. <- p.n(n)
   for(s in sd.s) {
      for(h in hh)
         isb.n[....] <- sum((f.n[....]-f.m(....))^2)* du
   }
 }
## 100-Stichproben - IV und MISE berechnen
f.n <- array(dim = c(length(uu), l.h, l.s, l.n),
             dimnames = c(list(NULL, p.h(hh)) ,dimnames(Samp)))
iv.n <- array(0,dim = c(....), dimnames = c(list(p.h(hh)),dimnames(Samp)))
mise.n <- array(0,dim = c(l.h, l.s, l.n), dimnames = c(....))
for (k in 1:n.sim) {
   ## Stichprobe erstellen
   for(n in n.s)
      for(s in sd.s)
        Samp[....] <- c.sample(....)
  ## Schaetzer, IV und MISE berechnen:
   for(n in n.s) {
      n. <- p.n(n)
      for(s in sd.s) {
         for(h in hh) {
            f.n[, p.h(h),p.sd(s), n.] <- with(.... , bandwidth = h)$est
         }
      }
   }
   for(n in n.s) {
      n. <- p.n(n)
      for(s in sd.s) {
         s. <- p.sd(s)
         for(h in hh) {
            iv.n[p.h(h), s., n.] <- iv.n[....] +
              sum((f.n[.... , p.sd(0), ....]-f.n[.... , s., ....])^2)* du
            mise.n[....] <- mise.n[....] + sum((f.n[....] - f.m(uu))^2)* du
         }
      }
   }
}
iv.n   <- iv.n  /n.sim
mise.n <- mise.n/n.sim

## Graphik
library(sfsmisc)
mult.fig(2,main = "MISE=ISB+IV")
yl <- range(isb.n[,"sd=0","n=40"], iv.n[,"sd=0","n=40"], mise.n[,"sd=0","n=40"])
matplot(isb.n[,"sd=0","n=40"], type = "l", col = 1, ylab = "y",
        xlim = c(0,11), ylim = yl)
....

yl <- range(isb.n[,"sd=0.6","n=40"],iv.n[,"sd=0.6","n=40"],mise.n[,"sd=0.6","n=40"])
matplot(isb.n[,"sd=0.6","n=40"],type = "l",col = 1,ylab = "y",
        xlim = c(0,11),ylim = yl)
....