セットアップ日記

_ 多項式フィット

実験に必要なデータの処理は、rubyとRを使ってやっている。単純な計算はrubyで、グラフなどを扱うときにはRと使い分けている。カーブfittingもRで行うのだが、これが意外に面倒である。Rは統計のためのソフトなので、それほど複雑なfittingは想定していないのであろう。 まだ、Rに不慣れなためということもあるのだが、多項式でのfittingは、次数を変えるためにいちいち書き換えなければならなくて、煩雑に感じる。また、今回は、二変数でfittingしたかったので、さらに複雑になった。結局うまくいったfittingはこんな感じ。

result <- lm( z ~  
1+ x +I(x^2) +I(x^3) +I(x^4) +I(x^5) +
y+ I(y*x) +I(y*x^2) +I(y*x^3) +I(y*x^4) +I(y*x^5) +
I(y^2)+ I(y^2*x) +I(y^2*x^2) +I(y^2*x^3) +I(y^2*x^4) +I(y^2*x^5)
)

ここで、I()はfitting特有の式の解釈をしないために必要になるようです。result$coefで係数を取り出して、もとのデータとの整合性をチェックする。 あとはrubyで計算させれば良いわけだが、fittingに使ったxに関して五次、yに関して二次の式を計算するのに、rubyでは以下のように記述した。

coef=[
[ -1.741370e+00,  7.970617e+00, -4.912943e-01,  3.099862e-02, -6.283200e-04,  5.164492e-06 ],
[ -2.927893e-01,  7.592103e-02, -7.253879e-03,  3.186369e-04, -6.577961e-06,  5.171831e-08 ], #y
[ -1.592628e-02,  1.488542e-03, -3.180673e-05, -1.340528e-06,  6.018192e-08, -6.055244e-10 ], #y**2
]
sum=0.0
coef.each_with_index{|a,i|
  sum+=y**i*a.reverse.inject{|s,c| s*x+c}
}

このとき、injectを二回入れ子にしたら、エラーがでた。当然なのだが、一回目のときに、arrayと数値を計算しようとするためだ。しかし、これを回避する方法も思いついた。

coef=[
[ -1.741370e+00,  7.970617e+00, -4.912943e-01,  3.099862e-02, -6.283200e-04,  5.164492e-06 ],
[ -2.927893e-01,  7.592103e-02, -7.253879e-03,  3.186369e-04, -6.577961e-06,  5.171831e-08 ], #y
[ -1.592628e-02,  1.488542e-03, -3.180673e-05, -1.340528e-06,  6.018192e-08, -6.055244e-10 ], #y**2
0.0
]
coef.reverse.inject{|ss,a|
  ss*y+a.reverse.inject{|s,c| s*x+c}
}

として、最高次の項を0にしておけば、常に数の計算にできるので、多項式部分を単純にできる。しかし、どちらが良いのか微妙なところではある。今回はもう書いてしまったので、初めの方法で良いことにしよう。