#----------------------------------------------------------
#
#  Polynomial fit with matrix formalism
#
#  GdA 18/4/2016
#----------------------------------------------------------

polynom <- function(x, coef) {
  y = rep(0, length(x))
  for (i in 1:length(coef)) {
    if (!is.na(coef[i])) y = y + coef[i]*x^(i-1)
  }
  return(y)
}


poly.fit <- function(x, y, np) {
  # Build up design matrix
  X <- numeric()
  for(i in 0:np) X <- cbind(X, x^i)
  # E[beta]
  E.beta <- as.vector( solve( t(X) %*% X ) %*% t(X) %*% y  )
  # Covariance matrix for unitary sigma
  Cov.beta <- solve(t(X) %*% X )
  std.beta <- sqrt(diag(Cov.beta))
  Cor.beta <- Cov.beta/outer(std.beta, std.beta)
  return(list(E=E.beta, Cor=Cor.beta))
}

# example
beta.0 <- c(1, -2, 0.5, -0.01)
np <- length(beta.0) - 1       # degree of polynomial
sigma = 0.2

x <- seq(1, 5, len=11)         # predictors
n <-length(x)

#X <- numeric()                 # use design matrix 
#for(i in 0:np) X <- cbind(X, x^i)
#y <-  as.vector( X %*% beta.0 ) + rnorm(n,0,sigma) 

Dx <- max(x) - min(x)
# many points to plot the curve 'pl': plot
x.pl <- seq(min(x)-Dx*0.2, max(x)+Dx*0.2, len=101)
y.pl <- polynom(x.pl, beta.0)
plot( x.pl, y.pl, ty='l', col='blue', xlab='x', ylab='y',
      ylim=c(min(y.pl)-2*sigma, max(y.pl)+3*sigma) )

y <-  polynom(x,beta.0) + rnorm(n,0,sigma) 
points(x,y, pch=19, col='blue')

beta <- poly.fit (x, y, np)
points(x.pl, polynom(x.pl, beta$E), ty='l', lty=2, col='red')

print(beta.0)
print(beta$E)
print(beta$Cor)