Saturday, January 23, 2016

matplotlib入門

これまでグラフを描く時にはRを使ったりpythonでmatplotlibでやったりと色々かじっていたのですが、2ヶ月くらい前にRで2万枚くらいのグラフを描いていた結果、この言語で続けるのはもう無理だという気分になりました。だって変数のスコープがないんだもん。

それでmatplotlibに移ることにしたのだけれど、Rでやってた基本的なことするのにいちいちstackoverflowの記事を読むのが大変だったので、ごく簡単なところをここにメモすることにしました。

情報としては、matplotlib.org本家があるのと、Wes McKinneyの"Python for Data Analysis"を持っています。本の方はMatplotlibに限らずnumpy, pandasなど色々幅広いです。

使い始めてすぐ気になるのは、plt.plotでできることと、subplotのax.plotでできることが似ている一方、違っている点で、結局matplotlibのソースを見て納得しました。pltの方はpyplotモジュールの関数で、axの方はaxesの下の_subplots.AxesSubplotあたりのクラスの関数なのですね。要するに別のもの。

1. 数字のカンマ区切り(1000ごと)。みんなグラフの数字はカンマで区切りましょう!excelでも「数値」を選べばそうなります。カンマをつけてくれないと、いちいちゼロの数を数えるのが面倒です!私は以下のようにFuncFormatterを使ったcommaformを渡すようにしています。

2. 色の変え方。系列ごとにblue, redなどで変える方法はあちこちに出てます。数字によって連続的に変える場合は以下のようにnormalizationと、colormapを使います。

3. colorbarで色の意味を示すのも以下のようにします。例では2の累乗で表示していますがこの例だとデフォルトは10の累乗になります。コメントアウトするとわかります。colorbarの書き方は以下のようにする方法と、コメントアウトしてあるplt.colorbar(sc)
の方法で設定する方法があるようです。

4.  最近もっとも驚いたのは、グラフの上にマウスを持ってきて、'l'(エル)とか'L'と打つと画面上でy軸x軸がログスケールに変わることです。みんな知ってるのかなあ。それはそうとして以下の例ではコメントアウトしましたが、ax.get_xaxis().set_major_locatorでbase=2でx軸をlogに変えることができます。

5. matplotlibはスタイルシートをサポートするようになりました。冒頭コメントアウトしてありますが、plt.style.availableで使えるスタイルシートの一覧が見られます。自分のスタイルシートは
$HOME/.matplotlib/stylelib/foobar.mplstyle
のようなファイルに入れておくと読んでくれます。

 6. matplotlibをいじってみたい気になったので、ソースをgit cloneしてdevelopモードで動くようにしました。ダウンロードしたソースに対してpython setup.py developとすると
/usr/local/lib/python3.5/site-packages/matplotlib.egg-link
というファイルが出来、そこにcloneしたソースへのリンクが書かれています。このローカルコピーで動くことになります。





以下が例。背景が黒いのはdark_backgroundスタイルシートのせい。



Friday, January 22, 2016

壁にぶつかるボールが円周率を作る

2015年2月号の数学セミナーの時枝先生という方の「ちいさい数理みつけた」という記事の中に、壁にぶつかる2つの物体の話が出てくる。

最初、mは静止しており、Mがぶつかってくる。左端には壁があって、床はツルツルしている。mは弾き飛ばされて壁とMに何度もぶつかる。その回数が問題。

質量m, Mの間に $\frac{M}{m} = 100^d,  (d=0,1,2...) $の関係がある時に、衝突の回数が円周率のd桁になるという不思議な話だ。
この記事にはそのことしか書いていないので、自分で確かめよう。

久しぶりに高校の物理を思い出して、m, Mそれぞれの速度をu, vとし、エネルギー保存の法則と、運動量保存の法則を解けば良いことに気づく。衝突後のu,vをそれぞれu',v'とすると、
$mu + Mv = mu' + Mv'$
$\frac{1}{2} mu^2 + \frac{1}{2} Mv^2 = \frac{1}{2} mu'^2 + \frac{1}{2} Mv'^2$
壁に向かう方の符号を正としている。
これを解くと、以下のようになる。

$ \begin{pmatrix}
u' \\
v'
\end {pmatrix}
= \frac {1}{m+M}
\begin{pmatrix}
2 M & m - M \\
M-m & 2 m
\end {pmatrix}
\begin{pmatrix}
u \\
v
\end {pmatrix}
$

mの速さは壁にぶつかっても変わらず、符号だけ変わる。m,Mが衝突すると、mは壁にもぶつかってふたたびMに衝突する。

問題なのは停止の条件で、Mが押し返されて右向きに動き始めてからも、Mの速さよりmが上まわればまた後ろから衝突する。mの方がMより遅くなっても、一回だけ壁にぶつかる。

回数をカウントするプログラムを書いたのが以下。$r = \frac{M}{m} $ とし、$ v = 1$ とした。

結果、以下のように円周率に従っている。

0 3
1 31
2 314
3 3141
4 31415
5 314159
6 3141592
7 31415926
どこにも円周率らしきものがないのに、結果として円周率になっているのが不思議である。

停止の条件が自分ではよくわからず、時枝先生の記事にはロシア人らしき人の名前が引用されているだけだったので、元ネタの所在も分からなかったのだけれど、「衝突 円周率」で調べているうちにNY Timesの記事が見つかって、そのアニメを見て停止の条件がわかった。というかこのアニメ、前に見たような気がする。


-----
(2/6/2016)
2つ後の記事に書いたがブログにtexで数式が書けるようになったので、式のところを直した。

Saturday, January 16, 2016

山の高さの測り方

太陽系最大の火山は火星のOlympus Monsで、こちらを読むと以下のように書いてあります。
Olympus Mons is a shield volcano 624 km (374 mi) in diameter (approximately the same size as the state of Arizona), 25 km (16 mi) high, and is rimmed by a 6 km (4 mi) high scarp.
 高さが25kmあって、広がりは624km。だいたいアリゾナ州と同じ大きさというのですが、その例えだと大きいのか小さいのかわからないですね。しかし日本の本州の幅が200kmだということを考えれば、尋常でない大きさだとわかります。高さが25kmというのも想像を絶しています。

この記事のもうちょっと下を見ると、地球における最大の火山が紹介されていて、それはハワイのMauna Loaだというのです。
To compare, the largest volcano on Earth is Mauna Loa. Mauna Loa is a shield volcano 10 km (6.3 mi) high and 120 km (75 mi) across. 
 高さが10kmというわけですが、ちょっと考えればこれはエベレストより高いことがわかる。なぜそうなのかを考えると、通常我々が知っている山の高さは海抜であって、火星には海がないので、測り方が違っているのであろうということにすぐ気づくわけです。

Mauna Loaが10kmは海底からからだろうなと思って探すと、こういうのが見つかります。理屈としては、この火山は海の底から続いているので、底から測るのが正しいというわけです。で、その数字は以下のように17km。
Thus, the total relief of Mauna Loa, from its true base to its summit, is about 17,170 m (56,000 ft).
10kmという話はどうなってしまったのでしょう。不思議なものだと思って、さらにここの「世界で一番高い山」という記事ですが、

1. エベレスト 海抜8850m
2. Mauna Kea 海底から10km
3. Chimborazo 地球の中心から6384km

と書いてある。 Mauna Loaの立場はどうなってしまったのだろうと思いますが、地球の中心から、っていうのも面白かった。