良いもの。悪いもの。

スピログラフに初めて触れたのは小学校に上がる前だった。当時はそれがスピログラフという名称であることを知らず、輪っかをぐるぐる回して綺麗な模様を描ける道具としてお気に入りの玩具になっていた。自分はマンデルブロ集合やジュリア集合のような数学的に表現できる図形や模様にともて魅力を感じるのだが、幼少のころのスピログラフがその切っ掛けだったのかもしれない。 長方形の定規に大きさの異なる3枚の歯車が付属しており、そして2つの円がくり抜かれている。さらに、矢印、五角形、半円などの型が子供心をくすぐったものだ。因みに、歯車と円に付いている歯数は刻印されており、3つの歯車はそれぞれ36、52、63、円の方は96、105の歯を持っている。 幼少のころ以来、スピログラフを見かけることはなかったのだが、最近になって100円ショップでも売られているという噂を聞き、少し探ってみた。近くの100円ショップでは扱っていなかったが、「デザイン定規」という名称でネットで78円で売っていたので思わず購入してしまった(当然だが送料のほうが高くついた…)。 スピログラフは一般的には内トロコイド(hypotrochoid; 内余擺線)といい、以下の式で表現できる(Wikipediaより)。 定円の半径 r c 、動円の半径 r m 、描画点の半径 r d とし、下図では、それぞれ、青、緑、赤の線で示している。そして、回転角を θとした軌跡がスピログラフの模様となる。 スピログラフでは半径よりも、円周の歯の数で示した方が分かりやすいだろう。円周と半径の比は変わらないのだから、歯の数の比がそのまま半径の比となる。例えば定円の歯が30で動円の歯が10ならば、 r c : r m =3:1ということだ。 さて、幼少のころはこんなに綺麗な模様が数式で表現できるなどと夢にも思わなかったわけだが、今はそれを理解し、コンピュータを使って自由に描くことができる。コンピュータであれば、動円が定円の外側に来る外トコロイド(epitrochoid; 外余擺線)も描くことができるし、描画点を動円の外に出すことも思いのままだ。 せっかくコンピュータを使って簡単に表示できるようになったというのに、プログラムを組むことが難しくては意味がない。しかし、Python+matplotlibであれば簡単に...

世の中のことをもっと知るにはどうしたら良いだろうと思うときがある。世の中の多くの事柄はログやデータに落とされる。Googleなどの検索サイトは良い例だろう。さて、そのログやデータをどうすれば良いのか? 多くの場合、視覚化が有効な手段となる。 まずは身の回りの日常的なデータやログを何とかしたい。ただ、日常のデータを視覚化するのに数十行以上のコードは書きたくない。まるで息をするかのごとく自然に視覚化を行いたいのだ。そのためには1～2行、長くて数行で済ませることが必要だ。そこでPythonとmatplotlibを使う。加えて、IPythonがあればなお良い。IPythonの導入については以前のブログ記事である IPythonの埋め込みプロットが素晴らしい を参考にして欲しい。 まずは事前にnumpyとmatplotlibをインポートしておく。できればscipyも。 >>> from numpy import * >>> from pylab import * 短いコードで視覚化を行うためには、Pythonの内包表記は必須だ。例えば、5, 2, 1, 5, 8をデータとするグラフを書きたいのならIPythonを使って以下の1行で実現できる。 >>> plot([5, 2, 1, 5, 8]) 数値が行ごとにinput.txtというファイルに書かれていた場合は以下の通り。 >>> plot([int(x) for x in file("input.txt")]) map関数を使えばもっとスマートに書ける。ファイル内の数値を文字列として配列で読み込み、それらの文字列をmap関数によりintで整数に変換する。 >>> plot(map(int, file("input.txt"))) さて、ログが2次元で書かれていた場合はどうするのか。例えば以下のデータがinput.txtに書かれていたとする。 2, 2.5 5, 6.2 6, 3.6 7, 6.3 10, 1.9 この場合、以下のようにすればプロットできる。ファイルから行ごとの文字列を読み出して、それを","を区切りとしてリス...

matplotlibで3次元プロットする際の雛形を示す。たったこれだけ。 from pylab import * import matplotlib.axes3d as p3 x = [2.36, 1.33, 1.10, 1.33, 2.36] y = [10.0, 8.00, 6.00, 4.00, 2.00] z = [1.10, 1.20, 1.44, 2.07, 4.30] ax = p3.Axes3D(figure()) ax.plot3D(x, y, z, color=(0.0, 0.0, 1,0)) ax.scatter3D(x, y, z, color=(1.0, 0.0, 0.0), marker="o", s=10) show() 以下はもう少しコメントと装飾を加えたソースコード。 from pylab import * import matplotlib.axes3d as p3 # データ. x = [2.36, 1.33, 1.10, 1.33, 2.36] y = [10.0, 8.00, 6.00, 4.00, 2.00] z = [1.10, 1.20, 1.44, 2.07, 4.30] ax = p3.Axes3D(figure()) # ラベル. ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") # 青線を引く. ax.plot3D(x, y, z, color=(0.0, 0.0, 1,0)) # サイズ10の赤丸で表示. ax.scatter3D(x, y, z, color=(1.0, 0.0, 0.0), marker="o", s=10) # 表示範囲. ax.set_xlim(1.0, 3.0) ax.set_ylim(1.0, 12.0) ax.set_zlim(1.0, 5.0) # DPI=100でPNGファイル(plot3d.png)を作成. savefig("plot3d.png", dpi=100) # 画面に表示. show()

matplotlibでは、polyfitを使えば簡単に多項式によるフィッティングを行える。図ではオリジナルのデータを赤でプロットし、それに対し、1次(黒)、2次(緑)、3次(青)の多項式によるラインを描いている。 以下、ソースコード。 from pylab import * x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] y = [1, 1, 2, 3, 5, 6, 8, 9, 8, 6] coeffs1 = polyfit(x, y, 1) coeffs2 = polyfit(x, y, 2) coeffs3 = polyfit(x, y, 3) besty1 = polyval(coeffs1, x) besty2 = polyval(coeffs2, x) besty3 = polyval(coeffs3, x) plot(x, y, "ro") plot(x, besty1, "-k", linewidth=2) plot(x, besty2, "-g", linewidth=2) plot(x, besty3, "-b", linewidth=2) grid(True) legend(["Data", "Order=1", "Order=2", "Order=3"], "best", numpoints=3) show()

matplotlibの線の色、線種、線幅の指定方法を記しておく。 線の色はplot上で、青なら"b"、緑なら"g"のように指定する。ここでは、白い線を見せるために背景を灰色にしている。 以下、図の作成に使ったソースコード。 from pylab import * axes(axisbg="#777777") # 背景を灰色に. x = arange(-20, 20, 0.3) plot(x+1, x, "b") # 青. plot(x+2, x, "g") # 緑. plot(x+3, x, "r") # 赤. plot(x+4, x, "c") # シアン. plot(x+5, x, "m") # マゼンタ. plot(x+6, x, "y") # 黄. plot(x+7, x, "k") # 黒. plot(x+8, x, "w") # 白. plot(x+9, x, color="#77ff77") # エメラルドグリーン. legend(["blue", "green", "red", "cyan", "magenta", "yellow", "black", "white", "#77ff77"], "lower right") xlim(2, 12) ylim(-5, 5) show() 次に、線種の指定方法を示す。 以下、ソースコード。 from pylab import * x = arange(-20, 20, 0.3) plot(x+1, x, "-") plot(x+2, x, "--") plot(x+3, x, "-.") plot(x+4, x, ":") plot(x+5, x, ...

一昔前は2次元しか扱えなかったが、最近では3次元プロットも利用できるmatplotlib。 Cookbook/Matplotlib/mplot3D あたりを読めばなんとなく使い方はわかる。利用するには、以下のようにモジュールを読み込んでおく。 import pylab as p import matplotlib.axes3d as p3 ネットで調べてみたら scipy+matplotlibによる３元連立常微分方程式の軌道 について書かれているサイトを発見。ふむふむ、参考になる。 Gnuplot.py は知らなかったな。でも最終更新が2003年と古い。numpyではなく古いNumericを使っているっぽい。やはり、matplotlibがいいなぁ。でも、plot3Dだとうまくいかないようなことが書かれている。そうなのか…と思って試してみた。例示されているコード中のscatter3Dをplot3Dに変更するだけ。見た目はまったく問題ないように見えるけど、何か問題があるのかな? バージョンの違い? どちらにしろ自分としてはこれで十分かな。 余談だが、現時点で最新のscipy-0.5.2とnumpy-1.0.3の組み合わせでは正常にコンパイルできない。 svc co http://svn.scipy.org/svn/numpy/trunk numpy svc co http://svn.scipy.org/svn/scipy/trunk scipy 上記のようにsvcで最新のバージョンを取ってくるか、numpy-1.0.2を利用しなければならない。

「 カオス―新しい科学をつくる 」は15年も前に読んだ本なのだが、非常に面白かった。以前からマンデルブロ集合やローレンツアトラクタなどに興味を持っていたのだが、素人にもわかりやすく記述されており、世界にはなんて奇妙で不可思議な現象があるのだろうと素直に感じた。事実は小説より奇なりとはよく言ったものだ。こんなこと、事実を目の当たりにしなければ、到底信じられない。 例えば、x(next)=rx(1-x) という式がある。これは一見何の変哲もない式に見えるが、rのパラメータを変化させていき、xがどの値に落ち着くか観察してみるといい。とても奇妙な振る舞いが見られる。xの初期値を0.4としたとき、r=2.8ならばxの値は1つに落ち着くが、r=3.2だと2つ、r=3.5だと4つ、r=3.7だと無数になる。そうかと思えば、r=3.84だと3つに収束する。なんとも不思議だ。こういった予測できない振る舞いがカオスなのである。 上記のグラフはPython+matplotlibで描いたが、そのコードを書いておく。 from pylab import * def calc(r, x, n): values = [] f = lambda r, x: r * x * (1.0 - x) for i in xrange(n): x = f(r, x) values.append(x) return values if __name__ == "__main__": n = 1000 x = 0.4 for r in arange(2.8, 4.0, 0.001): plot([r for i in xrange(n-n//10)], calc(r, x, n)[n//10:], "b.", markersize=1) axis([2.8, 4.0, 0.0, 1.0]) text(2.9, 0.2, r"$x_\rm{next} = rx(1-x)$", fontsize=30) xlabel(r"$r$", fontsize=20) ylabel(r...

matplotlib はPythonで利用できるグラフ作成モジュールで、高機能で使いやすい。とても便利なので是非ともインストールしておこう。 以下のようにモジュールを読み込み、matplotlibを利用できるようにする。 from pylab import * 直線を描く。 plot([1, 2, 3]) show() 続けてグラフを描くとき、以前のグラフを残したままだと同じ図にプロットされてしまう。以前の図の内容を消去したい場合は、 clf() とすればよい。 グラフの範囲をaxisで指定し、x, yに対応したリストでポイントに赤丸を打ち、赤い線で結ぶ("ro-")。 plot([1, 2, 3, 4, 5], [5, 10, 20, 50, 100], "ro-") axis([0, 6, 0, 100]) show() axisの代わりに、 xlim(0, 6) ylim(0, 100) でも可。 二次曲線を緑の線("g")で描く。タイトルもつける。 x = arange(0.0, 10.0, 0.01) title(r"Quadratic Curve", fontsize=20) plot(x, x**2, "g") show() 正弦曲線を描く。タイトルとラベルをつける。タイトルやラベルにはTeXの形式で数式を書くこともできる。 t = arange(0.0, 2.0, 0.01) s = sin(pi * t) plot(t, s) title(r"Sine Curve", fontsize=20) xlabel(r"$\theta$", fontsize=15) ylabel(r"$\rm{sin}(\theta)$", fontsize=15) show() 正弦曲線のグラフ内に二次曲線のグラフを埋め込む。さらに、EPS形式でグラフをファイルに書き出す。 t = arange(0.0, 2.0, 0.01) s = sin(pi * t) plot(t, s) title(r"Sine Curv...