眠たくて頭痛いけど, せっかく伊庭先生がMCMCについて話してくださるというのに参加しないという選択枝は無い.
こないだ社内勉強会で川合さんによる統計熱力学でのモンテカルロの話を聞いたばかりだしな.
というわけで会社へ.
次の魔方陣の個数を拡張アンサンブルモンテカルロ法を用いて求める話もとても面白かった.
他の話も聞きたかったけど残念ながらここで帰宅.
インテルソフトウェア・カンファレンス2011に行ってきた.
無料なのに軽食が出て, それがえらくおいしかった. というのはおいといて.
Cilk Plusの配列構文はおもろいな.
C[0:n] = A[0:n] + B[0:n];
if (A[:] > 0) { B[:] += C[:]; }
みたいなの. 関係ないけど最初"Click Plus"と読んでて「なんでシリックって発音してるのかな」と思ってたのは内緒.
parallel_forはさりげなくVC2010でもppl.hでさりげなく対応してたりするし今度ちゃんと使ってみよう.
コンパイラにさまざまな情報を与えられる#pragmaや__assume文. VCのはいまいち(書いてもコード生成に殆ど影響なさげ)だけど, iccのは相当強力らしい.
__assume_aligned(p, 32); // pが32byte alignmentされている __assume(n % 8 == 0); #pragma loop count avg(32); // ループの平均回数が32回
VTuneの使い方のデモもほどよくディープでよかった. float→intキャストが遅かった時代のasm { fld x; fistp ix; }はよく使ったなあ.
期待以上のカンファレンスだった.
初めてのSymPy
昨日と同じことを今度はSymPyでもやってみる. 一部うまくいかない(分からない)ところがあったけど一応答えは出た. SymPy歴30分なのでよりよい方法があれば教えてください.
SymPyはPythonのライブラリなのでPythonをインストールしてから入れる. pythonを起動して対話的に
from sympy import *
まず使いたい変数を指定する. Symbolかvarを使う. varだとまとめて定義できて便利そう.
x = Symbol('x')
t = Symbol('t')
a = Symbol('a')
b = Symbol('b')
c = Symbol('c')
d = Symbol('d')
e = Symbol('e')
# または
var('x t a b c d e')
Pythonでは巾乗は^ではなく**なのに気を付ける.
f=integrate((exp(x)-(a+b*x+c*x**2+d*x**3+e*x**4))**2,(x,0,t))
solveに渡す式は等式ではなく関数なのに注意. 関数の零点を求めるということ. 最初ちょっと悩んだ.
ans=solve([diff(f,a),diff(f,b),diff(f,c),diff(f,d),diff(f,e)],[a,b,c,d,e])
答えは辞書で返ってきてるのでans[a]とかで参照する. t = 0での極限を出してみる.
>>> limit(ans[a], t, 0) 1 >>> limit(ans[b], t, 0) 1 >>> limit(ans[c], t, 0) 1/2
うまく行ってる感じ. ところがans[d]やans[e]はエラーになった.
>>> limit(ans[d], t, 0)
File "<stdin>", line 1, in <module>
File "C:\p\Python\lib\site-packages\sympy\series\limits.py", line 116, in limit
(...)
raise ValueError('(%s).series(%s, n=8) gave no terms.' % (f, x))
ValueError: ((-4*_w**4*exp(_w) + 5*_w**4 + 84*_w**3*exp(_w) + 96*_w**3 - 780*_w*
*2*exp(_w) + 840*_w**2 + 3720*_w*exp(_w) + 3840*_w - 7560*exp(_w) + 7560)/_w**8)
.series(_w, n=8) gave no terms.
テイラー展開か何かが最大n = 8で決め打ちなのかな. 何かで変更できるのかちょっと調べたけど分からなかったのでとりあえず進める. 近似値をえるにはevalf()を使う. が, ans[a]に対してどうevalすればよいのか分からなかったので代わりにlimitを使った. 分母が0にならないところの値なので大丈夫だろう.
>>> limit(ans[a],t,log(2)/4096).evalf(50) 1.0000000000000000000000045896459498312533391085579 >>> limit(ans[e],t,log(2)/4096).evalf(50) 0.041670192358256038810322743017020389259233306763901
さて, もとに戻って上記limitのエラーの部分を見てみる. 実装コードをすぐ見られるのは利点やね. sympy/series/gruntz.pyの508行目.
def calculate_series(e, x, skip_abs=False, logx=None):
""" Calculates at least one term of the series of "e" in "x".
This is a place that fails most often, so it is in its own function.
"""
f = e
for n in [1, 2, 4, 6, 8]:
series = f.nseries(x, n=n, logx=logx)
if not series.has(O):
# The series expansion is locally exact.
return series
series = series.removeO()
if series:
if (not skip_abs) or series.has(x):
break
else:
raise ValueError('(%s).series(%s, n=8) gave no terms.' % (f, x))
return series
うむ, どうみても8までしかやってない. しかも数値固定. 試しにforループのところを10まで増やしてみる. お, 今度はうまく行った. なんでパラメータで変更できるようになってないのだろう.
>>> limit(ans[d],t,0) 1/6 >>> limit(ans[e],t,0) 1/24
全般的に速度はMaximaよりもっさりという感じ. Mapleは流石にそれらよりずっと速い(10倍は違うと思う). でもPythonなのは何かと使いやすいなあ. 最後に全体のコード
from sympy import *
var('x t a b c d e')
f=integrate((exp(x)-(a+b*x+c*x**2+d*x**3+e*x**4))**2,(x,0,t))
ans=solve([diff(f,a),diff(f,b),diff(f,c),diff(f,d),diff(f,e)],[a,b,c,d,e])
limit(ans[e],t,log(2)/4096).evalf(50)
limit(ans[e],t,0) ## コードをいじらないとエラー
初めてのwxMaxima
こないだのx86/x64最適化勉強会の私の発表[高速な倍精度指数関数expの実装]のおまけのところでやってた数式処理ってどうするのという質問があったので書いてみる. 私はMapleを使ってたのだけどMaximaはフリーソフトで同じことができる. なお私はMaxima歴30分なのであしからず.
今回やりたいことはt = log(2) / 4096が与えられたときに
f:=∫0t(exp(x) - (a + bx + cx2 + dx3 + ex4))2dx
を最小とするようなa, b, c, d, eは何か?という問題(expの計算ではb=1固定でやった)を手で計算せずにコンピュータにさせること. まずは積分させてみよう.
integrate((exp(x)-(a+b*x+c*x^2+d*x^3+e*x^4))^2,x,0,t);
と入力してみる. 関数をxについて区間[0, t]で積分せよということ. おっと入力の終わりは[Enter]ではなくctrl-[Enter]なので注意. すると
Is t positive, negative, or zero?
と聞いてくる. へぇ. positiveと打つと無事計算結果が出た. 最初からtが正を仮定しておくには
assume(t > 0);
としておけばよいようだ. 関数の定義は := を使う.
f(a,b,c,d,e):=integrate((exp(x)-(a+b*x+c*x^2+d*x^3+e*x^4))^2,x,0,t);
fの最小値を求めるにはfをa, b, c, d, eで微分してそれらを0と置いた連立方程式を解けばよい. 微分はdiffを使う. たとえばaで微分するには
diff(f(a,b,c,d,e),a);
でOK. 連立方程式を解くにはsolveを使う. 複数の式は[]でくくる. たとえば
solve([2*x+3*y=13,4*x+y=11],[x,y]);
とするとx = 2, y = 3と表示される. 今回欲しい値はtの関数なのでそれをans(t)と定義してみよう.
ans(t):=solve([diff(f(a,b,c,d,e),a)=0,
diff(f(a,b,c,d,e),b)=0,
diff(f(a,b,c,d,e),c)=0,
diff(f(a,b,c,d,e),d)=0,
diff(f(a,b,c,d,e),e)=0],[a,b,c,d,e]);
えらく長い式が出力された. 正しく計算できてるかt → 0での値を見てみよう. 極限値を取るにはlimitを使う.
limit(ans(t),t,0); [[a=1,b=1,c=1/2,d=1/6,e=1/24]]
ちゃんとexp(x)をテイラー展開したときの最初の4項の係数がでている. よさそうだ.
数値の定義は:を使う.
t : log(2)/4096
こうしてからans(t)とすると厳密解が表示される. 今回欲しいのは小数の値なのでfloatを使ってみる.
float(ans(t)); [[a=-7482132.582755132, b=2.646673732478273*10^11, c=-7.036349605814977*10^15, d=6.479223200922161*10^19, e=-1.912918493656289*10^23]]
あれ, 値があきらかにおかしい. 実はMapleでも同じことが起こって少し悩んでいた. 精度を大きく取らないと誤差が大きくてよい近似値にならないようだ. 内部演算の精度を変更するにはfpprecに値を設定すればよい. 任意精度の計算にはbfloatを使う.
fpprec : 50;bfloat(ans(t)); [[a=1.0000000000000000000000045896015893185326124785262b0, b=9.9999999999999999918636653907610524181757630578724b-1, c=5.0000000000003365597465341327601927431869990904259b-1, d=1.6666666613631816947174247945213363211134669284956b-1, e=4.1670192348702792353452213368172890252007837484349b-2]]
今度はちゃんと出た. 今回の操作を全て書いてみよう.
assume(t>0);
f(a,b,c,d,e) := integrate((exp(x) - (a + b*x + c*x^2 + d*x^3 + e*x^4))^2, x, 0, t);
ans(t):=solve([diff(f(a,b,c,d,e),a)=0,
diff(f(a,b,c,d,e),b)=0,
diff(f(a,b,c,d,e),c)=0,
diff(f(a,b,c,d,e),d)=0,
diff(f(a,b,c,d,e),e)=0],[a,b,c,d,e]);
fpprec : 50;bfloat(ans(log(2)/4096));
これだけでよい. すばらしい. 今検索するとはじめてのMaximaという膨大な(1000ページ超!)Maximaと数学の解説書を書いてる方がいらっしゃった. ジャンルが多岐に渡っていて目次を見ているだけでため息が出る.
わたしのモチベーションが下がったというのもあるのだが（まあ、それが唯一の原因であるのだけど）って書いてるやん. 開催メリットが低くなっちゃったらどうしようもない. いくら開催コストを下げるアイデアを出したとしても, さまざまな理由で却下されるだろう. ここは初心に戻って最新kernelのソースのどこかを読んでモチベーションを回復させるか, あるいは回復するまでしばらく放置でええんじゃないかな.