Initial pull

Author: roberto maestre

DESCRIPTION

@@ -0,0 +1,16 @@
+Package: additiveRegressionSpline
+Type: Package
+Title: C++ implementation of regression splines with academic purpose
+Version: 0.1.0
+Author: Roberto Maestre
+Maintainer: Roberto Maestre <roberto.maestre@bbvadata.com>
+Description: C++ implementation of examples from "Generalized Additive Models: An Introduction with R, Simon N. Wood".
+License: GPL-2
+LazyData: TRUE
+Imports: Rcpp (>= 0.11.0)
+LinkingTo: Rcpp,RcppArmadillo
+Depends:
+    R (>= 3.0.3),
+    data.table,
+    ggplot2,
+    mgcv

NAMESPACE

@@ -0,0 +1,3 @@
+useDynLib(additiveRegressionSpline)
+exportPattern("^[[:alpha:]]+")
+importFrom(Rcpp, evalCpp)

R/RcppExports.R

@@ -0,0 +1,27 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+hatMatrix <- function(m, length) {
+    .Call('additiveRegressionSpline_hatMatrix', PACKAGE = 'additiveRegressionSpline', m, length)
+}
+
+cubicSpline <- function(x, knots) {
+    .Call('additiveRegressionSpline_cubicSpline', PACKAGE = 'additiveRegressionSpline', x, knots)
+}
+
+splineModel <- function(x, knots, y) {
+    .Call('additiveRegressionSpline_splineModel', PACKAGE = 'additiveRegressionSpline', x, knots, y)
+}
+
+splineModelPenalized <- function(y, x, knots, lambda) {
+    .Call('additiveRegressionSpline_splineModelPenalized', PACKAGE = 'additiveRegressionSpline', y, x, knots, lambda)
+}
+
+splineModelPMatrix <- function(x, z, knots_x, knots_z) {
+    .Call('additiveRegressionSpline_splineModelPMatrix', PACKAGE = 'additiveRegressionSpline', x, z, knots_x, knots_z)
+}
+
+splineModelPMatrixFit <- function(y, x, modelMatrix, sX, sZ, lambda_x, lambda_z) {
+    .Call('additiveRegressionSpline_splineModelPMatrixFit', PACKAGE = 'additiveRegressionSpline', y, x, modelMatrix, sX, sZ, lambda_x, lambda_z)
+}
+

R/gam.R

@@ -0,0 +1,76 @@
+rk <- function(x,z) {
+  ((z-0.5)^2 - 1/12)*((x-0.5)^2 - 1/12)/4 -
+    ((abs(x-z)-0.5)^4-(abs(x-z)-0.5)^2/2 + 7/240) / 24
+}
+
+spl.X <- function(x,knots){
+  q <- length(knots) + 2                # number of parameters
+  n <- length(x)                        # number of observations
+  X <- matrix(1, n, q)                  # initialized model matrix
+  X[,2] <- x                            # set second column to x
+  X[,3:q] <- outer(x, knots, FUN = rk)  # remaining to cubic spline
+  X
+}
+prs.fit <- function(y, x, knots, lambda) {
+  q = length(knots) + 2       # dimension of basis
+  n = length(x)               # number of observations
+  Xa = rbind(spl.X(x, knots), mat.sqrt(spl.S(knots))*sqrt(lambda))  # augmented model matrix
+  y[(n+1):(n+q)] = 0          # augment the data vector
+  lm(y ~ Xa - 1)              # fit and return penalized regression spline
+}
+spl.S <- function(knots) {
+  q = length(knots) + 2
+  S = matrix(0, q, q)                           # initialize matrix
+  S[3:q, 3:q] = outer(knots, knots, FUN=rk)     # fill in non-zero part
+  S
+}
+mat.sqrt <- function(S){
+  d = eigen(S, symmetric=T)
+  rS = d$vectors %*% diag(d$values^.5) %*% t(d$vectors)
+  rS
+}
+GVC <- function(y,x,knots) {
+  n <- length(x)
+  # GCV (R)
+  mod <- prs.fit(y, x, knots, lambda = 0.0001)
+  trA <- sum(hatvalues(mod)[1:n])
+  trA
+  rss <- sum((y- fitted(mod)[1:n])^2)
+  rss
+  
+  (length(x)*rss)/(length(x)-trA)^2
+}
+trunc <- function(x, ..., prec = 0) base::trunc(x * 10^prec, ...) / 10^prec
+
+am.setup<-function(x,z,q=10)
+  # Get X, S_1 and S_2 for a simple 2 term AM
+{ # choose knots ...
+  xk <- quantile(unique(x),1:(q-2)/(q-1))
+  zk <- quantile(unique(z),1:(q-2)/(q-1))
+  # get penalty matrices ...
+  S <- list()
+  S[[1]] <- S[[2]] <- matrix(0,2*q-1,2*q-1)
+  S[[1]][2:q,2:q] <- spl.S(xk)[-1,-1]
+  S[[2]][(q+1):(2*q-1),(q+1):(2*q-1)] <- spl.S(zk)[-1,-1]
+  # get model matrix ...
+  n<-length(x)
+  X<-matrix(1,n,2*q-1)
+  X[,2:q]<-spl.X(x,xk)[,-1]           # 1st smooth
+  X[,(q+1):(2*q-1)]<-spl.X(z,zk)[,-1] # 2nd smooth
+  list(X=X,S=S)
+}
+fit.am<-function(y,X,S,sp)
+  # function to fit simple 2 term additive model
+{ # get sqrt of total penalty matrix ...
+  rS <- mat.sqrt(sp[1]*S[[1]]+sp[2]*S[[2]])
+  q.tot <- ncol(X)                # number of params
+  n <- nrow(X)                    # number of data
+  X1 <- rbind(X,rS)               # augmented X
+  y1 <- c(y,rep(0,q.tot))         # augment data
+  b<-lm(y1~X1-1)                  # fit model
+  trA<-sum(influence(b)$hat[1:n]) # tr(A)
+  norm<-sum((y-fitted(b)[1:n])^2) # RSS
+  list(model=b,gcv=norm*n/(n-trA)^2,sp=sp)
+}
+
+

additiveRegressionSpline.Rproj

@@ -0,0 +1,20 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source