2009年
7月
1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31

setup diary

2007|12|
2008|01|02|03|04|05|06|07|08|09|10|11|12|
2009|01|02|03|04|05|06|07|08|09|10|11|12|
2010|01|02|03|04|05|06|07|08|09|10|11|12|
2011|01|02|03|04|05|06|07|08|09|10|11|12|
2012|01|02|03|04|05|06|07|08|10|11|12|
2013|01|02|03|04|05|06|07|08|09|10|11|12|
2014|01|02|03|04|06|08|11|
2015|01|02|03|04|05|06|07|08|10|11|12|
2016|01|02|03|04|05|06|07|08|09|10|11|12|
2017|01|02|03|04|05|06|07|08|09|10|11|12|
2018|01|02|03|04|05|06|07|08|09|10|11|12|
2019|01|02|03|04|05|06|07|08|09|10|11|12|
2020|01|02|03|04|05|06|07|08|09|10|11|12|
2021|01|02|03|04|05|06|07|08|09|10|11|12|
2022|01|02|03|04|05|06|07|08|09|10|11|12|
2023|01|02|03|04|05|06|07|08|09|10|11|

2009-07-09 Rでfitting

_ 多項式フィット

実験に必要なデータの処理は、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にしておけば、常に数の計算にできるので、多項式部分を単純にできる。しかし、どちらが良いのか微妙なところではある。今回はもう書いてしまったので、初めの方法で良いことにしよう。