2008年05月04日

ScilabとRの処理速度

やることは山積みだがヒマ感がすごいのでRとScilabで処理速度を比較してみた。

RもScilabも繰り返しは苦手っぽいけど繰り返しの処理速度は気になるし、 分かりやすいのでEular法で微分方程式を解かせてみた。

問題は 「ゼロから学ぶ物理数学 (ゼロから学ぶシリーズ)」 より拝借。食べると1秒あたり4パーセント身長が縮むキノコをアリス(身長125cm)が食べたとき、 身長の変化はどうなるのか?という問題。

0.005秒ごとに値を計算させ、ベクトルとして記録、プロットまで終わらせた時点の時間を計測。用いたマシンはOS: Windows Vista、CPU:Core 2 Duo 2.20 GHz、RAM:2GByte。

まずはScilabから。

timer()
x = 0;
y = 125;
i = 1;
dX = 0.005;
X = 0;
Y = 0;
while x < 100
X(i) = x;
Y(i) = y;
y = y-y*0.04*dX;
x = x+dX;
i = i+1;
end
plot(X,Y)
timer()

入れたばかりで使い方が怪しいけど、これで実行時間が分かるはず。

つぎにR。

f <- function(){
 x <- 0
 y <- 125
 i <- 1
 dX <- 0.005
 X <- numeric(0)
 Y <- numeric(0)
 while(x < 100){
X[i] <- x
Y[i] <- y
y <- y-y*0.04*dX
x <- x+dX
i <- i+1
}
 plot(X,Y)
 }
system.time(f())

system.time関数の使い方がよく分からないので一連の処理を関数にまとめてから実行時間を計った。

結果。

  • Scilab: 7.46
  • R:       2.81

もうちょっとくらい調べてからやればよかったけど、MATLAB系の言語ってのはもともとそんなに早くなくて、 特に繰り返し文を多用すると遅さが顕著になるらしい。ちなみに、Rでは

X <- numeric(0)
Y <- numeric(0)

この部分を

X <- numeric(20000)
Y <- numeric(20000)

といった具合に必要なベクトルの長さをちゃんと指定する形に直してやると、実行速度は1.15秒まで縮まる。

何にも指定していない状態だとグラフはScilabの方が綺麗かな。gnuplotぽい。

RにしろScilabにしろ演算をベクトル同士に帰着できればもっと早くできるはず。画像処理なんかはどうなんだろう。 Rでやったときはプロットがやたら重くなったけど。またヒマなときにやろう。

posted by Rion778 at 14:34 | Comment(6) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
1年も前の記事にコメント失礼します。
通りすがりのものです。

Scilabについてですが、微分方程式でしたら関数odeを使って計算させると劇的に高速化されますよ。

function res=alice(t,y)
res=-y*0.04;
endfunction

timer()
x0 = 0;
y0 = 125;
dX = 0.005;
X=0:dX:100;
t0=0;
Y=ode(y0,t0,X,alice);
plot(X,Y)
timer()

私の環境ですと、これで
 15秒→0.5秒
になりました。

おそらくRでも同等の関数があると思いますよ。
Posted by IC at 2009年06月14日 19:30
有益な情報ありがとうございます。
あまりScilabについては詳しくないので助かります。
Rも正直あまり詳しくないというか、微分方程式を解くようなことはあまりしないので分からないのですが、たしかに同様の関数あるかもしれませんね。
また調べてみようかと思います。
Posted by Rion778 at 2009年06月20日 17:42
いまさらの通りすがりですが、ちょっと気になったことがあります。
普通に考えると、1/200秒刻みで指数関数を計算するなら、一秒あたりの倍率に0.005をそのままかけるのではなくて、0.005乗するのではないでしょうか。
Posted by 通りすがり at 2010年04月23日 15:19
>>通りすがりさん

ご指摘ありがとうございます.まさか2年前の記事にコメントが付くとは思っていなかったのとこのブログは更新しなくなって久しいので気付くのが遅れました.

指摘の件ですが,確かに与えられた条件(1秒あたり4%身長が縮む)は比較的単純なので次の形の解析解を求めることが可能です.
y = 125*exp(-0.04*x)
この式を単純にプロットするだけならば単純に0.05刻みの値をxに与えれば良いだけです(ご指摘にある「0.005乗」というのはこのような解釈でよろしいでしょうか?).

しかし,今回の例はあくまでオイラー法の練習ですので,解析解を求めてしまっては意味がありません.与えられている条件は
x_0 = 125
dy = -0.04*y*dx
の2つだけです.この2つの条件のみを使って微分方程式の近似解を求めようというのがオイラー法と呼ばれる方法です.ですから,上記の条件をdx = 0.005として差分化した式をなんども繰り返すことになります.このとき,yの差分dyは上記の式から-0.04*y*0.005になります.

納得して頂けたでしょうか?「オイラー法 微分方程式」などのキーワードで検索するとより詳しい情報が得られると思います.
Posted by Rion778 at 2010年05月09日 18:03
またまた今更ですが、コメントありがとうございます!

オイラー法の概念的には正しいのですが、この例題は指数関数になっているので更新則の式が間違っているのではないでしょうかという意味でした。言葉足らずで誤解を招いてしまったようです。
年利200パーセントのサラ金が、毎月の増加率に換算すると、2^(1/12)になることを伝えたかったのです(汗


よって、この例では、

新y = 旧y - 旧y*(1-(0.04)^dX)

が正しいのではないかと思います。
Posted by 通りすがり at 2010年06月08日 12:10
>>通りすがりさん
コメントありがとうございます。

おそらく私が問題文を詳しく記載していないのが原因で勘違いをしておられると思うのですが、問題文そのものから直接に(つまり微分方程式を解析的に解くことなしに)指数関数を得ることはできないのです。そのため、オイラー法を適用する時点で指数関数が出現するというのは誤りです。

一度指数関数を導出するために必要な手順を説明させて頂きます。
この問題では
- 最初の身長は125cm(y0 = 125)
- 1秒当たりにすると4%縮む(dy/dt = -y*0.04)
という条件のみが与えられております。ここで注意しなければならないのは、1秒経過すると1秒前の身長に対して4%縮んでいるとは限らないという点です(この条件を満たそうとすると身長は離散的に変化しなければなりません)。条件が述べているのは「もし1秒前の時点で変化率がストップすれば1秒後の身長は1秒前の身長から4%縮んだものになる」という事で、実際には1秒の間にも変化率は時々刻々と変化するのだ、という点を考慮する必要があります。

まず変化率に関する条件を次のように変形します。
(1/y)dy = -0.04dt
t=0でy=y0、t=tでy=yとおいて両辺を積分しますと
ln(y/y0) = -0.04t
さらに両辺の指数をとれば
y/y0 = exp(-0.04t)
yに関して解くと
y = y0*exp(-0.04t)
という指数関数が得られます。この指数関数を得るためのステップを指して「解析解を得る」と呼んでおります。また、今回の例の場合ですと、得られる指数関数は通りすがりさんの例示されたような一般の指数関数ではなく、底がe(ネイピア数)である特別の指数関数になります。底がeである指数関数に対してはexp(...)という記述がよく使われます。

問題の条件によっては三角関数や対数関数が得られる場合もありますし、そもそも解析的に解くことが困難な場合もあります。解析的に解くことが困難な場合、指数関数などのような解析解を得ることを諦めて有限の時間のステップを設定して計算を繰り返す、すなわちオイラー法を行うことになります。

つまり、問題から指数関数が得られるという通りすがりさんの認識は正しいのですが、オイラー法というのはその指数関数(解析解)が得られないような問題に対して適用するものです。そして、今回の例題に対するオイラー法のプログラム中には指数関数は出現しません。なぜなら、yの変化量はdy/dt=-0.04yであり、yとその変化量の関係は比例だからです(もしここの関係が指数関数として与えられているなら当然指数関数を組み込むことになりますが、その場合与えられるプロットは全く違うものになります)。

また、年利200%のサラ金の場合ですが、これは金利が単利か複利かによって場合が異なります。単利の場合は時間と支払い残高の関係は「初期値*年利^時間(年)」なので通りすがりさんの式が正しいのですが、複利の場合はdy/dt=yとおいた上でこの微分方程式を解く必要があります。そしてこの場合は利子を計算する間隔を短くするほどに返済金額が増えていくのですが(利子が元金に組み入れられていきますので)、もしこの間隔を限りなく0に近づけることができたとすると、1年後の返済金額は2倍ではなくe倍になります。この「間隔を限りなく0に近づける」という操作が微分方程式を解析的に解くという操作に対応します。もし「間隔を0に近づけるが有限の値(例えば0.01)に留めておいてそのステップごとに返済額を計算する」だとすると、これがオイラー法に対応します。後者では原理的に指数関数が出現しえないことをご理解頂けると思います。
今回の例題は金利に例えると複利の話になりますので、その点で勘違いが生じたのかもしれません。
この点については以下のサイトで詳しく、私の説明などより分かりやすい解説がされておりますので、興味がおありでしたら一度ご覧下さい。
http://d.hatena.ne.jp/starbow/20070412/1176363152
Posted by Rion778 at 2010年10月10日 13:42
コメントを書く
お名前: [必須入力]

メールアドレス: [必須入力]

ホームページアドレス: [必須入力]

コメント: [必須入力]

認証コード: [必須入力]


※画像の中の文字を半角で入力してください。
この記事へのトラックバックURL(言及がない、関連がないと判断した場合削除することがあります)
http://blog.seesaa.jp/tb/95613472
※言及リンクのないトラックバックは受信されません。

この記事へのトラックバック
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。