多項式近似メモ

# Function to estimate polynomial coefficients and coordinates estimated
# 多項式近似の係数と、推定座標を返す関数
glmPolynomial<-function(y,x,dg,X){# y:dependent values (従属値列),x:descriptive values (説明変数),dg:degrees of polynomials (多項式の次数)
 # Xs: Matrix consisted of x^1,x^2,...x^{dg}
 # Xs: xの1,2,...,dg 乗の値でできた行列
 Xs<-matrix(0,length(X),dg)
 xs <- matrix(0,length(x),dg)
 for(i in 1:dg){
  Xs[,i]<-X^i
  xs[,i] <- x^i
 }
 # glm() is generalized linear regression function in R
 # glm() はＲの一般化線形回帰関数
 fit<-glm(y~xs[,1:dg]) 
 # beta: coefficients estimated
 # beta: 推定係数ベクトル
 beta<-fit$coefficients
 # Xs2 has an additional column for x^0
 # Xs2 には、x^0のための列を追加
 Xs2<-cbind(rep(1,length(X)),Xs)
 # predY is a matrix of values for individual degrees of x
 # predY は xの各次数の項
 predY<-t(t(Xs2)*beta)
 # apply(): Sum up all columns corresponding to 0 to dg degrees
 # apply()関数で0:dg次数のすべてを足し合わせる
 sumupPredY<-apply(predY,1,sum)
 list(beta=beta,x=X,predY=sumupPredY)
}

# Generate data
N <- 7
x<-seq(from=0,to=1,length=N) + rnorm(N)
y<-runif(N)+x^(N+2)*2+sin(x*6*pi)


# GLM is to be applied to multiple degrees (dgs)
# GLMを複数の次数に適用
dgs<-1:(N-1)
X <- seq(from=min(x)-0.2,to=max(x)+0.2,length=1000)
Ys <- matrix(0,length(dgs),length(X))
for(i in 1:length(dgs)){
 outglm<-glmPolynomial(y,x,dgs[i],X)
 Ys[i,] <- outglm$predY
}

xlim<-c(min(X),max(X))
ylim<-c(min(y)-10,max(y)+10)
par(mfcol=c(2,3))
for(i in 1:length(dgs)){
	plot(X,Ys[i,],xlim=xlim,ylim=ylim,col=dgs,type="l",main=i)
	points(x,y,col=2,pch=20,cex=2)
}
par(mfcol=c(1,1))