Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fitting Cubic Splines
#FITTING SPLINES IN R
# D Menezes 20170531
#0. PRELIMINARIES
#install.packages("pacman") #only needed if you don't have it installed
pacman::p_load("splines")
#1. INPUT VALUES
input.x<-c(seq(from=0,to=120,by=12)) #assuming we have annual spline points
input.y<-c(0,0.48,1.08,1.03,1.01,1,1,1,1,1,1) # % Developed say
reqd.x<-c(seq(from=0,to=120,by=3)) #suppose we want quarterly spline output
n=length(input.x)
#2 KEY IDEA
#build basic spline
option.1<-spline(input.x, 100* input.y) #simplest - fits one cubic between each input x value points
option.1b<-spline(input.x, 100* input.y,xout=reqd.x) #fits one cubic between each required x value
#Compare
plot(input.x, 100*input.y,
main = paste("LDF Example: Spline Fit vs Actuals"),
xlab="Months Developed",
ylab="% Dev",
pch=16,
cex=1,
xaxt="n")
axis(1, at=seq(0,120,by=12), las=1)
points(option.1b)
lines(option.1b, col="red")
legend('bottomright',
legend=c("Annual-input",
"Quarterly - output"),
pch=c(16,1)
)
#3 ALTERNATIVES: Other spline options
option.2<-spline(input.x, 100* input.y,n=200) #creates a cubic spline with 200 points
option.3<-spline(input.x, 100* input.y,method="fmm") #different
option.4<-spline(input.x, 100* input.y,method="natural")
option.5<-spline(input.x, 100* input.y,method="periodic")
option.6<-smooth.spline(input.x, 100* input.y) #too smooth
#compare:
plot(option.1$x, option.1$y,
main = paste("Alternative Examples: Spline Fit vs Actuals"),
xlab="Months Developed",
ylab="% Dev",
pch=16,
cex=1,
col=1)
lines(option.1,col=1, lwd=1)
lines(option.2,col=2, lwd=1)
lines(option.3,col=3, lwd=1)
lines(option.4,col=4, lwd=1)
lines(option.5,col=5, lwd=1)
lines(option.6,col=6, lwd=1)
legend('bottom',
legend=c("Default",
"200 point",
"FMM",
"Natural",
"Periodic",
"Smoothed"),
lty=1,
lwd=2.5,
col=1:6
)
#4 RETURNING THE VALUE OF ARBITRARY X VALUE - USE interSpline
option.7<-interpSpline(input.x, 100* input.y)
predict(option.7,x=57)
#5 Blog specific chart output
plot(input.x, 100*input.y,
main = paste("Raw data"),
xlab="Months Developed",
ylab="% Dev",
pch=16,
cex=1,
xaxt="n")
axis(1, at=seq(0,120,by=12), las=1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.