##
## 1-D smoothing
##

X = runif(200, min=0, max=10)
X = sort(X)
YT = log(X) + 1/(0.7+sin(X)^2)
Y = YT + 0.2*rnorm(200)
plot(X, Y)
lines(X, YT, col='orange')

## Simple Kernel smoothing.
Yhat = array(0, c(length(X),3))
S = c(0.05, 0.5,1)
for (k in 1:3) {
  for (i in 1:length(X)) {
    W = (X-X[i])^2/S[k]^2
    W = exp(-W/2)
    W = W/sum(W)
    Yhat[i,k] = sum(W*Y)
  }
}

lines(X, Yhat[,1], col='blue')
lines(X, Yhat[,2], col='green')
lines(X, Yhat[,3], col='yellow')

YL = loess(Y ~ X)
lines(X, YL$fitted, col='grey')


##
## Generalized additive models
##

library(gam)

X = rnorm(200*4)
X = array(X, c(200,4))

XC = X
XC[,1] = 2*X[,1]
XC[,2] = X[,2]^2
XC[,3] = log(1+X[,3]^2)
XC[,4] = 1/(1+abs(X[,4]))

YT = 1 + XC[,1] + XC[,2] + XC[,3] + XC[,4]
Y = YT + rnorm(length(YT))

g = gam(Y ~ s(X[,1]) + s(X[,2]) + s(X[,3]) + s(X[,4]))
c = g$coefficients

for (k in 1:4) {
  XX = cbind(X[,k], X[,k]*c[k+1] + g$smooth[,k], XC[,k])
  pdf(sprintf("gam%d.pdf", k))
  plot(XX[,1], XX[,2], ylim=c(min(XX[,2:3]),max(XX[,2:3])))
  ii = order(XX[,1])
  lines(XX[ii,1], XX[ii,3], col='orange')
  dev.off()
}

stop()
  
##
## Dimension reduction
##

## A single index model.
n = 200
X = rnorm(n*4)
X = array(X, c(n,4))
LP = X %*% c(1, -1, 2, -2)
YT = LP / (1 + abs(LP)^1.5)
Y = YT + 0.2*rnorm(length(YT))
m = lm(Y ~ X)

sir = function(Y, X) {

  ii = order(Y)

  ns = floor(length(ii)/20)
  M = array(0, c(ns, dim(X)[2]))
  for (j in 1:ns) {
    jj = seq((j-1)*20+1, j*20)
    M[j,] = apply(X[ii[jj],], 2, mean)
  }

  CM = cov(M)
  CX = cov(X)
  C = chol(CX)
  CI = solve(C)
  Q = t(CI) %*% CM %*% CI
  eg = eigen(Q)

  F = list()
  F$vectors = CI %*% eg$vectors
  F$values = eg$values
  return(F)
}

F = sir(Y, X)

## A double index model.
n = 200
X = rnorm(n*4)
X = array(X, c(n,4))
LP1 = X %*% c(1, -1, 2, -2)
LP2 = X %*% c(1.5, 1.5, 1.5, 1.5)
YT = (1 + LP1)^3/10 + (1 - LP2)^2
Y = YT + 10*rnorm(length(YT))

F = sir(Y,X)

LP1_est = X %*% F$vectors[,1]
LP2_est = X %*% F$vectors[,2]

print(cor(LP1, LP1_est))
print(cor(LP2, LP2_est))