Rでスプライン補完の実装

すべての観測点を通る場合のスプライン補完。
http://www.akita-nct.ac.jp/yamamoto/lecture/2004/5E/interpolation/text/html/node3.html
を参考に実装

x<-seq(1,10,length.out=20)
x1<-x
y<-3*x^4+5*x^3+x^2

plot(x,y)

h<-diff(x)
len<-length(h)
v<-6*(diff(y)[2:len]/h[2:len]-diff(y)[1:(len-1)]/h[1:(len-1)])


m<-matrix(nrow=len-1,ncol=len-1)
for(i in 1:(len-1)){
	if(i==1){
		 m[i,]<-c(rep(0,i-1),2*(h[i]+h[i+1]),h[i+2],rep(0,len-3-i+2-1))
	}
  	if(i>1 & i<len-1){
	   m[i,]<-c(rep(0,i-2),h[i],2*(h[i]+h[i+1]),h[i+2],rep(0,len-3-1-(i-2)))
	}
	if(i==(len-1)){
          m[i,]<-c(rep(0,i-2),h[i],2*(h[i]+h[i+1]))

	}

}

u<-solve((m),v)
b<-c(0,u,0)/2
a<-diff(c(0,u,0))/6/diff(x)
d<-y
c<-diff(y)/diff(x)-1.0/6.0*diff(x)*(2*c(0,u,0)[1:(length(c(0,u,0))-1)]+c(0,u,0)[2:length(c(0,u,0))])


f<-function(x){
	sum<-0
	for(i in 1:(length(x1)-1)){
		
		sum<-sum+dunif(x,min=x1[i],max=x1[i+1])*(x1[i+1]-x1[i])*(a[i]*(x-x1[i])^3+b[i]*(x-x1[i])^2+c[i]*(x-x1[i])+d[i])
　　　　　　
                        if(length(which(z==x1[i]))>=1 & i!=1){
				sum[which(z==x1[i])]<-sum[which(z==x1[i])]-d[i]
			}
	}

	sum
}

curve(f(x),col=2,add=T)