「パターン認識と機械学習」1章 多項式曲線フィッティングをRで書いてみた

事の発端

先日「PRML読書会第一回 - Seeking for my unique color.」のソースコードが一部界隈で局地的に話題にあがり，やれ

  • takousikiってなんだよロシア人かよ

だとか

  • evalが多すぎて引いた

だのといった罵詈雑言が（主に僕の口から）飛び交っていたわけだけれども，さすがに言うだけではsyou6162先生に申し訳ないと感じ，改めて彼のコードをRで書きなおしてみることにした．

多項式曲線フィッティングのコードを書き直す

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

f <- function(x){sin(2*pi*x)}

takousiki <- function(n){
  x <- numeric(0)
  z <- numeric(0)
  for(i in 1:n){
    x <- cbind(x, seq(0,1,by=0.1)^i)
    z <- cbind(z, seq(0,1,by=0.01)^i)
  }
  data <- data.frame(f(seq(0,1,by=0.1))+rnorm(11, mean=0, sd=0.5), x)

  colnames(data) <- c("y",paste("x",as.character(1:n),sep=""))
  colnames(z) <- paste("x",as.character(1:n),sep="")

  # predict
  fit <- lm(y ~ ., data=data)

  plot(data[,2], data[,1], ylab="", ylim=c(-1.5,1.5), xlab="x", main=paste("M=",n))
  curve(f(x), add=T, lwd=2)
  lines(seq(0,1,by=0.01), predict(fit,data.frame(z)), col="red", type="l")
}

使い方としては，

1
2
3

takousiki(1)
takousiki(3)
takousiki(8)

とすると，多項式で曲線フィッティングした結果がプロットされる．ただし切片だけの0次は書けない．

サンプルデータの生成と推定した曲線の描写用のデータの生成を同時に行なっているため，コードがちょっと見にくいかもしれない．基本的にx^2 やx^3 などの多項式の各項の値をあらかじめデータフレームとして作っておいて，lm()関数で回帰するときのモデルを簡略化している．あと，関数名はオリジナルを踏襲しているので，あしからず．

他の方法としては

lm()関数などに渡すモデルは，as.formula()関数を使うことによって文字列から生成することができる．この方法もかなり簡単に使うことができるので，データフレームの列名を規則的に付けている場合などは使い勝手は良いと思われる．

1
2
3

> n <- 5
> as.formula(paste("y~",paste("x",1:n,sep="",collapse="+")))
y ~ x1 + x2 + x3 + x4 + x5