Rで多項式回帰をするとき、モデル式を作るためにpolyという関数を使いました。
poly関数は変数と次数を与えると、その次数の多項式モデルを返してくれる関数なのですが、我々が普通にイメージする多項式回帰を行うためにはオプションraw=TRUE
を指定しなければなりません。
ではraw
は何を指定するオプションなのか、とヘルプを引くと
raw: if true, use raw and not orthogonal polynomials.
とorthogonal polynomials、すなわち直交多項式などというパワーワードが出てきます。
ヘルプの中には直交多項式が何なのかに関するヒントは全く記されていないのでググって調べるとチェビシェフ多項式、ルジャンドル多項式とかいう新しい用語がヒットします。
それらは微積分が混じった恐ろしげな数式で定義されており、その定義式を解読しするだけで中々しんどい。
そしてようやくそれらを解読し計算できるようになったところでふと気づくのは、結局それが多項式回帰とどんな関係があるのかは分からない。
ここらで心が折れてきます。
そこをなんとか私なりに理解をしたつもりになったので、その内容をここにまとめておきます。
最初に結論を述べますと、ここでいうn次の直交多項式モデルとは、観測ベクトルについて各要素を累乗することで得られるn個のベクトルを、平均と分散とで正規化し、Gram-Schmidt直交化によってベクトルとして直交化したn個のベクトルを説明変数とするモデルのことです。
【統計モデル】
n次元ベクトルx
をpoly
関数に与え、オプションをdegree=3
, raw=FALSE
として目的変数y
の線形回帰をするとします。
すなわち
#R lm(y ~ poly(x, degree = 3, raw = FALSE))
このとき、実際に行われる回帰のモデル式は
になります。
ただしx'は正規化されたxであり
です。
またx''はx^2をGram-Schmidt直交化によってx'と直交にしたベクトル
を正規化したベクトル
です。
同様にしてx'''もx^3をGram-Schmidt直交化によってx',x''と直交にしたベクトルを正規化したベクトルとなります。
degreeが3以上のときにもGram-Schmidt直交可は一次独立なベクトルがあるかぎり手続きを継続できるので、観測ベクトルのベキ乗の系列が一次独立であるかぎり目的の次数まで観測ベクトルの変形ができます。
逆に観測ベクトルのベキ乗系列が一次従属になれば関数poly
はエラーを返します。
【計算】
具体的にpoly
関数の働き、すなわちどういう計算をしているかを示します。
ベクトルx
をpoly
関数に与え、degree=3
, raw=FALSE
としたとき
#R poly(x, degree = 2, raw = FALSE)
の戻り値は行列になり、その行列は一列目が
#R b1 <- (x - mean(x)) / sqrt(sum((x - mean(x))^2))
で計算されます。なおこれは平均と分散とで正規化されたxです。
また二列目は
#R b2a <- x^2 - b1 * sum(x^2 * b1)/sum(b1^2) b2 <- (b2 - mean(b2a)) / sqrt(sum((b2a - mean(b2a))^2))
と上節で述べたようにx^2
をGram-Schmidt直交化によってx
と直交なベクトルを正規化したベクトルになります。
degree
が2よりも大きくなったとしても同じ要領で計算されます。
【余談】
上記の結論はpoly
関数の内部定義を読み解くことで見出しました。
内部定義を見ると分かるのですが、実際にはqr
という関数を使って多項式を直交化しています。
関数qr
は行列をQR分解するR関数で、QR分解とは行列を高速に上三角化する処理です。
ベクトルの直交化がqr
でできるというテクニックは使い道があるのかもしれません。