#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Tutorium zu Schaetzen und Testen 2 im SoSe 2016
# zu Tutorium 09 am 30.06.2016
# Autor: Almond Stoecker
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#### Truncated polynom basis ####

tp <- function(x, degree = 3, knots = 3, boundary.knots = range(x)  ) {
  
  breaks <- seq(boundary.knots[1], boundary.knots[2], length.out = knots + 2)
  
  B0 <- sapply( 0:max(degree-1, 0), function(d) x^d )
    
  B1 <- sapply( breaks[-length(breaks)], function(knot) (x-knot)^degree *(x>=knot) )
   
  B <- cbind(B0, B1)
  
  return(B)
  }


#### plot TP ####

x <- seq(0,1, l = 60)

TP <- tp(x, degree = 3)
dim(TP)

matplot(x, TP, type = "l")

set.seed(3023)
(gamma <- rnorm(ncol(TP)))

matplot(x, scale(TP, scale = gamma, center = FALSE), type = "l")

matplot(x, TP%*%gamma, type = "l")



#### B-spline basis ####

library(splines)

# Funktion bs aus splines muss leicht angepasst werden
bs2 <- function(x, degree = 3, knots = 3, boundary.knots = range(x)) {
  h <- diff(boundary.knots)/(knots+1)
  Boundary.knots <- c(boundary.knots[1] - degree*h, boundary.knots[2] + degree*h)
  breaks <- seq(Boundary.knots[1], Boundary.knots[2], len = knots+2+2*degree)
  
  B <- bs(x, knots = breaks[-c(1,length(breaks))], degree = degree, intercept = TRUE, Boundary.knots = Boundary.knots)
  B <- B[,(1+degree):(knots+2*degree+1), drop = FALSE]
  return(B)
}


#### plot BS ####

BS <- bs2(x, degree = 3)
dim(BS)

matplot(x, BS, type = "l")

matplot(x, scale(BS, scale = gamma, center = FALSE), type = "l")

matplot(x, BS%*%gamma, type = "l")


#### Transform BS into TP ####

# Gedanke: TP = BS %*% Omega -> 'solve(BS)%*%TP = Omega' -> ABER: BS, TP nicht symmetrisch

x0 <- seq(min(x), max(x), l = ncol(TP))

TPsym <- tp(x0, degree = 3)
dim(TPsym)

BSsym <- bs2(x0, degree = 3)
dim(BSsym)

Omega <- solve(BSsym) %*% TPsym

## und damit:

opar <- par(mfrow = c(1,2))
matplot(x, BS %*% Omega, type = "l")
matplot(x, TP, type = "l")
par(opar)

# oder andersherum

opar <- par(mfrow = c(1,2))
matplot(x, BS, type = "l")
matplot(x, TP %*% solve(Omega), type = "l")
par(opar)



#### B-splines fitten ####

## Funktion aus letztem Tutorium:
  
  # Wahrer funktionaler Zusammenhang
  f <- function(x) {
    
    sin(2*(4*x - 2)) + 2*exp(-(16^2)*(x - 0.5)^2) + 2
    
  }
  
  n <- 200 
  
  x <- seq(0,1,length=n)
  
  # Zufallszahlen
  eps <- rnorm(n,mean=0, sd=0.3)
  
  y <- f(x) + eps

  plot(x,f(x), type = "l")
  points(x,y)
  
## mit gam fitten ##

library(mgcv)

gam1 <- gam(y ~ s(x, bs = "ps", k = 40))

plot(gam1, all.terms = TRUE)  
points(x,y-coef(gam1)[1], col = "cornflowerblue")
lines(x,f(x)-coef(gam1)[1], col = "cornflowerblue", lty = "dashed")  

## Penalisierung ##

# Differenzenmatrix
D <- diff(diag(ncol(BS)), diff = 2)
# Penalty Matrix
K <- t(D)%*%D

## Penalty Matrix für zur TP-Basis transformierte B-Spline-Basis 
round( t(Omega)%*%K%*%Omega, 5 )