読者です 読者をやめる 読者になる 読者になる

Kumpei Shiraishi's blog

物理、プログラミング、クライミングに関する雑記

DDSKKの統計をグラフにするコードを公開した

以前、sedとPython(pandas、matplotlib)を使って、DDSKKの統計をグラフにする手順をまとめた記事を書いた。

kumpeishiraishi.hatenablog.com

この記事に書いたような煩雑な手順を追わなくとも、Pythonとmatplotlibだけでグラフを描けるようにした。

github.com

詳しくはREADMEを参照。チンケなプログラムだが、あまり顧みられていない分野を自分で見つけてきて、それを埋めるプログラムを書いた初めての経験なので嬉しく思う(まあ、注目されていないということは、大して意味がないということではあるのだが)。皆さん、是非遊んでみてください。

SKKの統計をグラフにする

概要

SKKは、個人辞書に登録されている候補数の統計を記録することができる。この統計情報のグラフ化を試みる中で、以下のことを勉強することができた:

  • sedの使い方
  • pandasでcsvファイルを取り扱う方法
  • datetimeで日時データを取り扱う方法
  • matplotlib.datesで時系列データを表示する方法

この記事は、その個人的な記録である。グラフ化の元ネタはこのページ:~/.skk-record « 余談ですが……

結論

グラフ

Candidates in my SKK

コード・手順

  1. % cp ~/.skk-record ~/temp/rec.csvで、作業フォルダに統計情報をコピーする
  2. % sed -i '' -e 's/ 登録.*語数: /,/g' rec.csvで、日時と語数の列からなるcsvファイルに整形する(BSD sed)
  3. 以下のPythonコードを実行する:
import datetime as dt
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

# データフレーム
df = pd.read_csv('rec.csv',header=None)
df.columns = ['datetime','words']

# x軸時系列データのリストを整える
dates = [None]*df.shape[0]

for i in range(df.shape[0]):
    dstr = df.ix[i,'datetime']
    ddate = dt.datetime.strptime(dstr, '%a %b %d %H:%M:%S %Y')
    dates[i] = dt.datetime(ddate.year,ddate.month,ddate.day,ddate.hour,ddate.minute,ddate.second)
    # 当たり前だが、ここで時、分、秒を格納しなければ、もっと粗いグラフになる

# 以下、オブジェクト指向なmatplotlibの使い方をしている
# まず基本のFigureクラスのインスタンス。キャンバスだと思っておけばよい。
fig = plt.figure()

# 次にAxesクラスのインスタンス。白い背景と軸くらいに思っておこう。
ax = fig.add_subplot(1,1,1)

# axにプロットする。labelは凡例の名前を指定するため。
ax.plot(dates,df.ix[:,'words'],label='SKK jisyo candidates')

# 以下でグラフの体裁を整える
# 凡例の位置はlegend関数で指定できる
ax.legend(loc='upper left')

# axのx軸について、locatorとformatterを指定する
ax.xaxis.set_major_locator(mdates.DayLocator())
ax.xaxis.set_minor_locator(mdates.HourLocator(interval=12))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %d'))

# autofmt_xdate()で良い感じに回転させる
fig.autofmt_xdate(bottom=0.125,ha='center')

# 描画範囲を指定
start = min(dates)
end = max(dates)
ax.set_xlim(dt.date(start.year,start.month,start.day),dt.date(end.year,end.month,end.day+1))

# axにy軸のラベルを付加する
ax.set_ylabel('Number of words')

plt.savefig('a.png',dpi=300)

勉強したこと

sed

sedは置換、挿入、削除などのテキスト処理ができるコマンド。使い方は

を参照。ここでは今回使ったものだけを書いておく。

-eオプションは、その後に処理内容が来ることを表す。

-iオプションは、処理結果をファイルに上書きすることを表す。GNU sedとBSD sedで挙動が異なるので、後述する注意を参照。

特殊文字は、

  • .:任意の一文字
  • *:直前の文字が任意の個数続く文字列(0個も含む)
  • <space><space>*:連続するSPACEを表す
  • <tab>:TAB
  • .*:任意の文字列
  • \^:行頭
  • \$:行末

コマンドの前に来る数字をアドレスと呼び、その行のみに対して処理を行う。例:'3d'ならば3行目のみを削除する。

置換処理はs(Substitution)コマンドで、's/処理前文字列/処理後文字列/g'のように用いる。g(Global)オプションを渡さない場合、各行の最初に登場した当該文字列にのみ処理を行う。区切り文字はsコマンドの後に来たものなら何でも良い。

さて、sedにはGNU系とBSD系があり、それぞれで挙動が違うことに注意。Mac OS XにはBSD sedが入っているが、Homebrewでgnu-sedを入れればGNU sedを使うこともできる。BSD sedを使う上で注意が必要なのは、

  • -iオプションを付けて上書きする場合
  • 改行\nを出力する場合

である。

-iオプションに関してBSD sedのマニュアルには

Edit files in-place, saving backups with the specified extension.

と書いてある。つまりBSD sedにおいては、-iオプションの後ろで指定した文字列をファイル名の末尾に付けて元ファイルのバックアップを取るのである。例:sed -i 'boke' -e 's/foo/bar/g' hoge.txtを実行したならば、hoge.txtに置換済みの内容(欲しい内容)が保存され、hoge.txtbokeに置換前の内容(バックアップ)が保存される。GNU sedと同じように、余計なファイルを生成せずに上書きしたいなら-i ''と空の文字列を指定すればよい。

改行出力についてはこのページを参考にした。結論として、例えば1行目に文字列を挿入して上書きするには、

sed -i '' 1s/\^/foobar\\$'\n'/ input.txt

を実行しよう(今回はpandasでヘッダーを付加したので、BSD sedでの改行出力は不要なのだが)。

pandas

pandasはPythonでデータ分析をする際の定番ライブラリ。インストールしたものの、碌に使っていなかったので、ちょうど良い機会だった。

以下のページを主に参照した:

ここでは、今回のコードを読むのに必要なこと等を書いておく。

csvデータを読み込むには、read_csv('ファイル名')とする。第1行をヘッダーとして認識するので、これが不要ならheader=Noneを与える。また、読み込み時にヘッダーを指定するにはnames='ABC'を与える。そうすれば、第1列をAなどと名付けられる。

read_csv()でデータフレームdfを生成した後から、ヘッダーを指定して列名を名付ける(上書き)にはdf.columns = ['columnA','columnB','columnC']とする。

shape(行数, 列数)を、describe()で平均、分散などの基本的な統計量を得ることができる。

各要素にアクセスするには、ix[]を用いる。例:df.ix[2,'columnB']で第3行、第2列の要素を得る。列は列名でも列番号でも指定できる。複数選択は:を使う。

ソートはsort()。引数に列名を渡して、昇順でソートする。降順はascending=False

行の追加には、まず追加するデータフレームadd = DataFrame([['foobar',1]])を用意する。列名は追加先のデータフレームdfと同じものを、indexは追加先の末尾に合うものを指定する(然もなくばデータフレームの列が拡張され、形が変わる)。そしてconcat([df,add],axis=0)が追加後のデータフレームである。axis=1とすれば列の追加になる。

datetime

datetimeはPythonで日時データを扱うモジュール。使い方は

がよくまとまっている。

文字列から日時にするには

tstr = 'Sat Apr 23 07:12:25 2016'
tdatetime = dt.datetime.strptime(tstr, '%a %b %d %H:%M:%S %Y')

とする。日時から文字列にするには

tdatetime = dt.datetime(year=2015, month=10, day=21, hour=16, minute=29)
tstr = tdatetime.strftime('%Y/%m/%d')

とする。

指定子を知るにはここを参照。現在時刻を取得する.now()などの使い方も載っている。

matplotlib.dates

詳しくは以下を参照:

matplotlibは、オブジェクト指向なプロットを提供する。そのようにプロットする際の流れはコードのコメント、またはここに挙げたページを見てほしい。ここでは、matplotlib.datesの簡単な使い方について書いておく。

手作業で日時を目盛tickに指定するのはかなり大変だ。しかしmatplotlib.datesを使えば、tickの位置を制御するlocatorと、tick毎につける名前の書式を指定するformatterを、簡単に指定できる。

locatorでは分、時、日、週、月、年などの単位でインターバルを指定できる。例:matplotlib.dates.HourLocator(interval=6)は6時間毎にtickを記すlocatorである。他にも、毎週月曜日や、3年毎のイースターのようなtickを簡単に指定できる(公式ドキュメント参照)。

fotmatterでは、matplotlib.dates.DateFormatter()を使うことがほとんどだろう。例:matplotlib.dates.DateFormatter('%b-%d')はApr-24という形式で出力するformatterである。

これらをxaxis.set_major_locator()、xaxis.set_major_formatter()などに渡せばよい。

ちなみに、x軸のformatに日時を指定したが、多くの場合「2016/04/24」や「08:31:54」は横に長すぎて重なってしまう。autofmt_xdate()は、それを良い感じ回転、位置移動してくれる。デフォルトはbottom=0.2, rotation=30, ha='right'で、それぞれ下の余白、回転角、tickの相対位置を表す。また、今回はset_xlim()で描画範囲を指定した。指定しないと原点の下にx軸のtickが表示されず、ちょっと気持ち悪かったからだ。具体的には、x軸データの最大最小を取得し、両者の年月日のみを抜き出して、それを描画範囲とした。これならば、各日付の0時から24時まで描画されるので、欲しいtickになる。

matplotlibで定数関数を描く

やりたいこと

matplotlibで定数関数を描画したい。

ダメな例

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 1

x = np.arange(-1,1,0.1)
y = f(x)
plt.plot(x,y)

plt.show()

これだと

ValueError: x and y must have same first dimension

というエラーが出る。xyのデータ数が合っていない、ということだ。

上手くいく例

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 1

x = np.arange(-1,1,0.1)
y = [f(x[0])] * np.size(x)
plt.plot(x,y)
plt.show()

定数関数だと分かっているので、xの要素のどれかを指定してf(x)に渡し、xの配列長と同じ長さのリストをyに代入すれば良い。

ちなみにyの初期化は、こちらの方法の方がリスト内包表記よりも5.7倍くらい速かった。

matplotlibでグラフの場合分け

やりたいこと

matplotlibで、場合分けしたグラフを描きたい。

この記事では参考ページに倣って、

{
f(x) = \left\{ \begin{array}{ll}
x^2 & (0 \leq x \leq 1)\\
2-x & (1 < x \leq 2)\\
x^2 - 3x +2 & (2 < x \leq 3)\\
\end{array} \right.
}

を描く。

ダメな例

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    if 0<=x and x<=1:
        return x**2
    elif 1<x and x<=2:
        return 2-x
    elif 2<x and x<=3:
        return x**2-3*x+2

x = np.arange(0,3,0.001)
y = f(x)
plt.plot(x,y)

plt.show()

これを実行すると

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

というエラーになる。

関数f(x)の定義部分でxの大きさを判定しているのに、x = np.arange(0,3,0.001)というnumpyの配列を渡しているのが原因らしい。

上手くいく例

import numpy as np
import matplotlib.pyplot as plt

def f1(x):
    return x**2
def f2(x):
    return 2-x
def f3(x):
    return x**2-3*x+2

x1 = np.arange(0,1,0.001)
y1 = f1(x1)
plt.plot(x1,y1,"b")

x2 = np.arange(1,2,0.001)
y2 = f2(x2)
plt.plot(x2,y2,"b")

x3 = np.arange(2,3,0.001)
y3 = f3(x3)
plt.plot(x3,y3,"b")

plt.show()

明示的に別々の関数を定義すると上手くいく。複数の関数をmatplotlibで描画するには、単純にplotを複数回呼べばよい。グラフの色が別々だと格好悪いので、"b"を渡して全部青色にした。

そうして得られたグラフは以下。

Case in matplotlib

参考

3次元自己回避ランダムウォークをPythonで実装する

自己回避ランダムウォーク

自己回避ランダムウォーク(self avoiding random walk, SAW)とは、過去に占有した座標を回避しながら進むランダムウォークの一種である。3次元版のSAWは高分子のシミュレーションでよく用いられる。3次元では、以下のように漸近する。

{
\langle R^{2} (N) \rangle \sim N^{2\nu}
}

このとき{\nu \simeq 0.592}であることが知られている。

Rosenbluth-Rosenbluth最適化

さて、こんな3次元SAWをPythonでシミュレーションしてみようと思う。3次元直交座標を考え、原点から粒子がランダムウォークして動いていくとする。モノマー50個の長さのポリマーなら50回動くという要領で、動いた回数がポリマーの長さに対応するのである。

素朴な方法、次の格子をランダムに選び、それが自分の動いてきた履歴に含まれていたら終了、をそのまま実装したとする。ある程度の長さのポリマーを得るには、この方法では効率が悪すぎる。例えば長さ40のポリマーにまで到達する粒子は、わずか0.04%ほどである。そこで何らかの最適化が必要になる。

今回は、Rosenbluth-Rosenbluth最適化(Rosenbluth, Rosenbluth. J. Chem. Phys. 23, 356 (1955); http://dx.doi.org/10.1063/1.1741967)という方法を採る。この方法では次の格子をランダムに選ぶのではなく、6つの候補({x}正方向・負方向、{y}正方向・負方向、{z}正方向・負方向)から既に履歴に含まれているものを除外してランダムに選び、その代わりに重みをつける。

具体例を考えよう。{N-1}回の移動が終わり、{N}回目の移動を行う段階を考える。重み付けをする因子は{W(N)}である。6つの候補のうち、

  1. 履歴に含まれるものが1個(直前の位置)のみの場合、{W(N) = W(N-1)}
  2. 履歴に含まれるものが{m}個({m = 2, 3, 4, 5})の場合、{W(N) = \frac{m}{5}W(N-1)}
  3. 履歴に含まれるものが6個の場合、{W(N) = 0}として終了

以上の規則にしたがって、{W(N)}を決定する。但し、{W(0) = 1}である。

最後に、平均2乗変位{\langle R^{2} (N) \rangle}は以下の表式で求められる。

{ \displaystyle
\langle R^{2}(N) \rangle = \frac{\sum W(N)R^{2}(N)}{\sum W(N)}
}

但し、和は粒子番号についてとる。

Pythonコード

近日公開

SKKとispellの食い合わせが悪い

症状

DDSKKで書いた日本語と英語の混じった文章で、ispellによるスペルチェック候補を表示させようとすると、

Ispell and its process have different character map

というエラーとなり、表示できない。

原因と対策

~/.emacs.d/init.elに書いていた

(setq skk-isearch-start-mode 'latin)

が原因なので、

(setq skk-isearch-start-mode 'utf-8)

に書き換える。

コメント

インクリメンタルサーチ時はmigemoを使うので、SKKは引っ込んでいて欲しい。そんな場合、マニュアルではユーザ変数skk-isearch-start-modelatinにすべしとある。

これで欲しい動作は実現できるのだが、latin設定時に書いた日本語中のミススペルの正しい候補をispellで表示させようとすると、上に見たエラーを吐く。同様の症状はここなどで見られるが、何やら文字コード絡みの問題のようだ。

そこでutf-8に設定し直すと、当方の環境ではインクリメンタルサーチもスペルチェックも問題なしになった。

SKKは現在慣れようとしているところだが、未だに辞書ファイルの文字コードがEUCだったり、こういうマニュアルだったり、古式ゆかしきソフトであることよ。

日本で買うクライミング用品はなぜ高い?

現在私は、フィンランドに留学している。そして先日、ヨーロッパの有力なアウトドア用品通販サイトbergzeitでクライミング用品を購入する機会があった。金に汚いと思われるかもしれないが日本の定価を調べて計算したところ、今回は56%も安く購入することができた。

bergzeitはヨーロッパへの発送のみなので、日本から購入することは難しいが、最近は個人輸入で安く登山用品やクライミング用品を購入する人が増えているようである。また、海外旅行などの際に現地で安く買う人も多いだろう。こうした人々が共通して持っているのは、「日本の山道具は高い」という考えである。

この記事では、私が今回購入したクライミング用品に関して、ヨーロッパ価格と日本価格を比べてみたい。

なお、この記事で示した定価は2016年4月4日調べのものであり、使った変換レートは1ユーロ = 127円、1米ドル = 111円である。

価格表

品物 輸入代理店 bergzeit定価 輸入代理店定価 マージン率 備考
Petzl Sirocco アルテリア 9,599.93円 15,300円 59%
Petzl Grigri2 アルテリア 7,465.06円 12,000円 61%
Black Diamond Oz quickdraw ロストアロー 2,128.52円 2,400円 8% BD定価2,214.45円
DMM Phantom Screwgate carabinar ケーイーエム 1,488.44円 3,200円 115%
Wild Country Dyneema Sling 10mm 60cm モンベル 955.04円 1,300円 36%
  • 定価は全て税抜き価格である(bergzeitでは19%の付加価値税、日本では8%の消費税が加算される)。
  • マージン率とは、「(bergzeit定価と輸入代理店定価の差分)/(bergzeit定価)」で計算され、「ヨーロッパ価格の何割の下駄を履かせて販売しているか」を表す。
  • ただし、アメリカ企業であるBlack Diamondに関しては、備考に示したBlack Diamondウェブサイト上の定価を用いて計算した。

コメント

メーカーごとに見ていこう。

まずはフランスのクライミング、高所作業用品メーカーのPetzl。後述するBlack Diamondと並び立つ二大巨頭といった趣きがある。日本への輸入を請け負うのはアルテリアである。今回は軽量ヘルメットSiroccoとオートブロックビレイデバイスGrigri2を買ったが、どちらもヨーロッパ価格の60%くらいを上乗せして日本で販売しているようだ。

Black Diamondはクライミング用品メーカーの最大手と言っていいだろう。ヨセミテのクライミングがルーツなだけあって、ネイリングやビッグウォール関連のギアに関してはBDの独壇場である。このようなニッチなギアばかりではなく、カラビナやスリングを比較的低価格で販売しているので、ユーザーの裾野は広い。日本において、この低価格を実現しているのはロストアローだ。他にはクライミングシューズのScarpaやロープのBealなども扱っている。ロストアローの御大は山学同志会の坂下直枝、社員には倉上慶大、横山勝丘など有名クライマーの名も散見される。スポンサードクライマーに小山田大がおり、とにかくまぁ有力な代理店なのである。坂下直枝とBDの大ボスYvon Chouinardとのコネクションもあるようで、比較的良心的と言われる価格設定だ。実際、マージン率は抜きん出て低い。天晴れロストアローと言ったところか。

ウェールズのDMMを日本で取り扱うのは、ケーイーエムである。他にはロープで有名なedelweissも請け負っている。今回買ったPhantomに関して言えば、驚くなかれマージン率100%以上。ヨーロッパでクライミング用品を見ていると、DMMが飛び抜けて高いメーカーであるという印象はない。しかし日本では「DMM=かなり高い」という図式が出来上がっていると思う。その理由の一端はここにあると見て間違いないだろう。

最後に、ピークディストリクトの老舗Wild Countryである。輸入は彼のモンベル。モンベル本体は低価格で色々売っているという印象があるが、Wild Countryについても同様のようである。マージン率は、(ある意味異常な)ロストアローよりは高いが、アルテリア、ケーイーエムよりは低くなっている。

まとめ

Disclaimerとして、今回比較したクライミング用品は、各メーカーのラインアップのごく一部である。また、提示した「マージン率」という数字は、クライミング用品の値段を評価するのに多分に不適当である可能性はある。マージン率の分母に用いたbergzeitの定価に関しても、bergzeitがヨーロッパのほぼ全域をカバーする巨大な通販サイトであることを考慮すべきである。日欧のクライミング人口は桁違いで、それに伴う市場規模も全く異なる。さらに、世の中一般の輸入代理店がどれくらいの下駄を履かせているのかも、私はよく知らない。

しかし、そうは言っても、私にとっては驚きの数字だった。日本語でサポートを受けられるという利点はあるにせよ、これを考えると日本でギアを買うことは出来ないなあ。