Difference between revisions of "Bsplines.R"
From Organic Design wiki
(R code for plotting B-splines) |
(fixed bs bug) |
||
| (4 intermediate revisions by the same user not shown) | |||
| Line 1: | Line 1: | ||
# {{R}} | # {{R}} | ||
| − | + | library(splines) | |
| + | |||
# p38 | # p38 | ||
pos <- function(x) ifelse(x>0, x, 0) | pos <- function(x) ifelse(x>0, x, 0) | ||
| Line 7: | Line 8: | ||
knot <- 4 | knot <- 4 | ||
plot(x,y) | plot(x,y) | ||
| − | X <- cbind(1,x,x^2, x^3, pos(x-knot)^3) | + | #X <- cbind("int"=1,"x"=x,"x2"=x^2,"x3"=x^3, "xknot"=pos(x-knot)^3) |
| − | fit <- lm(y~X-1) | + | #fit <- lm(y~X-1) |
| + | |||
| + | fit <- lm(y~1+x+x^2+x^3+pos(x-knot)^3) | ||
xx <- seq(1,7, length=200) | xx <- seq(1,7, length=200) | ||
| − | XX <- cbind( | + | XX <- cbind("x"=xx) |
lines(xx, predict(fit, as.data.frame(XX))) # Bug here | lines(xx, predict(fit, as.data.frame(XX))) # Bug here | ||
abline(v=knot, lty=2) | abline(v=knot, lty=2) | ||
X | X | ||
| + | |||
# p41 | # p41 | ||
| − | knots <- c( | + | par(mfrow=c(3,1)) |
| − | atx <- seq(0 | + | knots <- c(1:3,5,7,8,10) |
| − | for(ord in | + | atx <- seq(0, 11, by=0.1) |
| − | B <- | + | for(ord in 2:4) { # ord cannot be 1 |
| + | B <- bs(x=atx, degree=ord-1, knots=knots, intercept=TRUE, Boundary.knots=c(0,11)) | ||
| + | # B <- predict(B, seq(0,11, length=101) # works | ||
matplot(atx, B[,1], type="l", ylim=0:1, lty=2) | matplot(atx, B[,1], type="l", ylim=0:1, lty=2) | ||
matlines(atx, B[,-1], col=1, lty=1) | matlines(atx, B[,-1], col=1, lty=1) | ||
Latest revision as of 03:10, 17 October 2007
Code snipits and programs written in R, S or S-PLUS library(splines)
- p38
pos <- function(x) ifelse(x>0, x, 0) x <- 1:7 y <- c(8,3,8,5,9,14,11) knot <- 4 plot(x,y)
- X <- cbind("int"=1,"x"=x,"x2"=x^2,"x3"=x^3, "xknot"=pos(x-knot)^3)
- fit <- lm(y~X-1)
fit <- lm(y~1+x+x^2+x^3+pos(x-knot)^3) xx <- seq(1,7, length=200) XX <- cbind("x"=xx) lines(xx, predict(fit, as.data.frame(XX))) # Bug here abline(v=knot, lty=2) X
- p41
par(mfrow=c(3,1)) knots <- c(1:3,5,7,8,10) atx <- seq(0, 11, by=0.1) for(ord in 2:4) { # ord cannot be 1 B <- bs(x=atx, degree=ord-1, knots=knots, intercept=TRUE, Boundary.knots=c(0,11))
# B <- predict(B, seq(0,11, length=101) # works
matplot(atx, B[,1], type="l", ylim=0:1, lty=2) matlines(atx, B[,-1], col=1, lty=1) title(paste("B splines of order", ord)) abline(v=knots, lty=2) }
