pureで統計
2014 アケオメ コトヨロ
去年頂いた『空飛ぶPython』を読んだ。Python３で頑張りましょって趣旨で書かれた本だ。 おいらは見よう見まねでPythonを使ってきたんで、こういう本は有り難い。炬燵に入って、 ぬくぬくしながら読みましたよ。
そして、むくむくとパソコンに向かうのでした。とは言え、いまだにデストリでデフォで Python3系をシステムスクリプトとして採用するって言う、勇気のある開発者はいないみたい だけどね。 そんな訳で、Python3が入ってる万次郎で試しつつ。本では、開発環境として、なんてったって アイドルってのを薦めていたけど、おいらはipythonさ。
[sakae@manjaro ~]$ ipython
Python 3.3.3 (default, Nov 26 2013, 13:47:45)
Type "copyright", "credits" or "license" for more information.
IPython 1.1.0 -- An enhanced Interactive Python.
? -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help -> Python's own help system.
object? -> Details about 'object', use 'object??' for extra details.
In [1]: a = ("Happy", "Python")
In [2]: a = ("Happy", "Ruby")
In [3]: a
Out[3]: ('Happy', 'Ruby')
あれ？ タプルは書き換え出来ないよって、本のあちこちで注意されてたけど、出来ちゃったよ。 これって得意の仕様変更が入ったの？ よう分からんわい。まてまて、幾らなんでも、もう少し 調べれ！
In [4]: a[1] = 'pure' --------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-4-07d8b85b7f99> in <module>() ----> 1 a[1] = 'pure' TypeError: 'tuple' object does not support item assignment
やっぱり、オイラーの早とちりだったわい。まだまだ修行が足りんな。でも何故タプルだけ 書き換え出来なくしたの？ 作者さんが関数族議員から脅されて、しぶしぶ実装した？
In [6]: dir(a) Out[6]: ['__add__', : '__rmul__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', 'count', 'index']
ここに出て来る、アンダーバーで固めた変数(?)が、Pythonの肝みたいだ。
In [7]: a.__str__ Out[7]: <method-wrapper '__str__' of tuple object at 0xb690340c>
gdbで開くと、オブジェクトそのものが見えるんだな。ちなみに、こやつはhaskellで言う showに相当するらしい。まあ、外枠は違っても、中は案外一緒っぽかったりして。
こんな暴言を吐くと、あの人が怒るかな。中身は一緒っぽかったりしても、それを取り巻く 応援団がいっぱい居るかが、言語の人気のバロメータだよと、かの本は言う。pythonは バッテリー。いろいろなライブラリィーが同梱されて出荷されてるから、余り余計な物は 入れなくても使えるよと。
便利そうなのの代表として、こんなのが出てた。
In [11]: import shelve
In [12]: lang = shelve.open("hoge")
In [13]: lang['ruby'] = ('matz', 'japan')
In [14]: lang['gauche'] = ('shiro', 'japan')
In [15]: lang.close()
Python流のデータ永続化。下請けにpickleモジュールを使って、簡単インターフェースを 実現してる。手軽に使えて便利そう。
便利だよー、だからみんなPythonで飛んでる世界を体験しよー。そして、disるのも忘れて いません。やっぱりあちらの人は一神教なんですな。 データ・サイエンスのプログラミング言語はRからPythonに置き換わる なんていう記事にも表れてきている。どうなんですかね？ 2013 Advent Calendarでも眺めてみれば、 今後の方向性が見えてくるかな。 で、おいらは去年に引き続き pure なんていう超マイナーなドイツ方言丸出しのHaskell もどきを使ってみる。なんたって、おいらは八百万神の国の人だもの、おおらか、おおらかに。
ああ、その前に、LispでOSを書く なんて刺激的なのを発見したぞ。これを試すには、gauche用の c-wrapperを入れれば良いのか。 後で、fedoraででも試してみるかな。
pure-csv
年が明け、還暦を過ぎて還暦プログラマの挑戦(Haskell に挑む→F#による言語造り) という状況。この方は病気のデパート風に伺えるけど、おいらも負けてはいないぞ。血圧高め という状況で、ずっと薬を飲んでる。
薬が切れそうなので、また貰ってこなければいけないんだけど、おいらがPythonで書いた 血圧記録・グラフかプログラム。一年毎の漬物にしてるんで、年初はデータが十分に溜まって いないのよ。（医者に行く時は、グラフを持って行くんだけど、そのグラフを有難がって コピーさせて下さいって言われる）
そんな訳で、今回は去年のまとめを持って行こう。一年の総括ですな。どんな総括すれば いいの？ そりゃ、医学 == 統計 でしょ。死亡率、罹患率、各種検査値からの異常発見、 これみんな統計技術ですよ。
最近患者が減ったなあ。じゃあ、患者を作ればいいんだな。誰でも分かる、腹回りなんてどうだ。85cm以上は 患者にしちゃえ。えと、腹回りと死亡率でも調べてみるか。カチャカチャカチャとEXCELを操作 する音。だめじゃな。そんじゃ、腹回りと 成人病の罹患率はどうだ。おおおー、優位な相関があるぞ。こりゃ使える。ってな具合。
なら、血圧データも統計処理したのを持ってけ。まあ、町医者に理解出来るかどうかは はなはだ疑問だけど。。おいらの知見の拡充にはなるな。プログラミングを通じてね。
Pythonからの出力は、いわゆるcsvファイルにしてる。
[sakae@pcbsd ~/src/pure]$ head -n 4 2013.csv 130101.04,118,73,57 130101.21,109,70,71 130102.04,128,82,60 130102.21,105,62,67
一行が４項のデータ。左から yymmdd.hhの測定日時、最高血圧、最低血圧、脈拍っていう並びだ。 これをpureで処理しよう。pureはマイナーゆえ、あれこれモジュールを同梱してない。 けど、有志がいろいろ作って提供してる。今回は、 Pure-CSV - Comma Separated Value Interface for the Pure Programming Language を使おう。
> using csv; > using namespace csv; > let p = open "2013.csv"; > let v = fgetr p; > close p; ()
普通にcsvファイルをオープンして、読み込んでからファイルをクローズする。超簡単っぽい。 どんな具合に読み込まれたか、確認。
> # v;
720
> take 4 v;
[{"130101.04","118","73","57"},{"130101.21","109","70","71"},{"130102.04","128","82","60"},{"130102.21","105","62","67"}]
朝に夜に血圧を測定してるんで、データは730あるはずだけど、５日分欠落してる。どこかへ 遊びに行ったとかで、データ無しなんだな。そんな事より、取り込んだデータが、文字列に なってるのが、いやらしいな。ちゃんと数値として読み込んで欲しいぞ。でも、測定日時が 数値になっちゃうと、困るな。だって、起床時と就寝時と分けて解析したいし。問題山積み。
マニュアルを精読すると、数値は数値らしく読み込むオプションが有るっぽい。そうすると、 日時のフォーマットを変えておかねば。 ついでに、csvファイルを指定してデータを読み込んでくれる関数を定義しておくか。
mycsv f = v with
d = dialect {quote_flag=>MINIMAL};
p = open (f, "r", d);
v = fgetr p;
z = close p;
end;
これでいいかな？
> let m = mycsv "2013.csv";
> take 2 m;
[{"130101:04",118,73,57},{"130101:21",109,70,71}]
うん、これでいいな。尚、オプションに指定した dialect って、方言って意味なのね。 csv.pureなんてファイルを見ると、EXCEL方言とか、RFC4180とかが登録されてた。
[sakae@pcbsd ~]$ rfc 4180
The Result:
4180 Common Format and MIME Type for Comma-Separated Values (CSV)
Files. Y. Shafranovich. October 2005. (Format: TXT=12931 bytes)
(Status: INFORMATIONAL)
RFCになったのは、以外に遅かったのね。CSVは適当に複雑なので、これを肴にRubyのあの人が 初Ruby本に解説してたな。RWH本にも有るぞ。HaskellのHackageにも山のように登録 されてる。こんなに有るのも迷惑だなあ。余り沢山選択子があると、どれも選ばれない ってのは、ビジネス界の常識ですだ。その点、pureは安心。ああ、話がそれたわい。
えと、少し統計関係の関数を定義しておくかな。
using math; // for sqrt
sum xs = foldl (+) 0 xs;
group n [] = [];
group n xs = (take n xs) : (group n (drop n xs));
avg xs = sum xs / #xs ;
ss xs = sum [(x - avg xs) ^ 2 | x = xs;] / ( #xs - 1) ;
sd xs = sqrt $ ss xs ;
median xs = if #xs mod 2 == 1 then
head $ drop (#xs div 2) $ sort (<) xs
else
(rs ! 0 + rs ! 1) / 2 with
rs = drop (#xs div 2 - 1) $ sort (<) xs
end;
sumは、Haskellにはデフォで付いているんだけど、pureには見当たらなかったので自前で定義。 ひょっとしてドイツ方言で、合計って定義されてるんですかね？
groupはHaskell方面で他の機能を実現してるかも知れないけど、リストを指定したエレメント 数で分割してくれるもの。ssは有名な分散ですな。各項を求める度に平均を求めてるけど、 これって無駄かな。まあ、一行野郎って事で許してくだせえ。
medianって、データを並べ替えした時の真ん中にある数値を取り出す演算。データ群の代表値 としてavgが普通に使われるけど、medeainの方が良いと仰る先生がいます。データ数が奇数と 偶数で場合分けしなきゃならないんだけど、pure君はoddもevenも持っていない。んな訳で、 さりげなく紛れ込ませておいたよ。
取り合えず、これだけあればいいかな。
> let am = [x | x = m; x ! 0 ! 7 == "0"] ; > let pm = filter (\x -> x ! 0 ! 7 == "2") m;
分析に先立ち、起床時と就寝時にデータを分けてみた。x!0 で、測定日時を選び出し、 そいつに!7を適用して、時の上位文字を指定してる。0なら起床時って訳。選び出すには、 条件を内包表記の中に書いてもいいし、ずばりfilter関数を使ってもよい。
> let hi = [x!1 | x = am; ]; > let low = [x!2 | x = am; ]; > let pls = [x!3 | x = am; ];
更に、起床時のデータを、最高血圧、最低血圧、脈拍に分解。
> avg hi; sd hi; median hi; 119.644444444444 7.1654713865796 120.0
算術平均、標準偏差、中央値を出してみた。無駄に下位桁まで表示しているな。 こういう時は、表示を制限しちゃおう。
> using system; > __show__ x::double = sprintf "%0.1f" x; > avg low; sd low; 77.7 4.1
上記は、一年分のデータのまとめだったけど、季節の変化で血圧って変動するの？ 24季節に分解してみる。
> map median $ group 15 hi; [126,126,125,124,119,122,121,117,119,117,117,116,112,113,120,117,115,118,118,122,120,123,128,124]
グラフでも書いてみないと、目がちかちかしちゃうな。もう少し荒くてもいいか。
> map median $ group 30 hi; [126.0,125.0,120.0,118.0,118.0,116.5,112.0,117.5,117.0,120.0,121.0,125.0]
うん、薬の飲んでいても、冬場は血圧高めなのね。よく分かるよ。
で、これだけじゃ、ちと味気ないので、相関とかに手を出してみる。最高血圧と最低血圧に 関係は有るのか、とかね。ぐぐったら、 高校数学C 統計学 なんてのが出てきた。 最近は高校でも統計学を教えるんですね。これって証券会社や日本競馬会がロビー活動して、 賭け事のイロハを若いうちに学んでおきましょって事だな。まあ、株のかけかたを研究して ノーベル賞を貰う御人が居るぐらいだから、世界的兆候なんでしょうな。
更にぐぐると、 平均と偏差､分散､相関の概念 なんてのが出てきて、このページの最後の式が簡単そうだったので、pureに落としてみる。
cor xs ys = (sum $ map tmc $ zip xs ys) / (#xs - 1) with
ax = avg xs ;
ay = avg ys ;
sx = sd xs ;
sy = sd ys ;
tmc z = ((z!0 - ax) / sx) * ((z!1 - ay) / sy) ;
end ;
２つのデータ系列を同時に扱うので、zipでまとめると、データの組（タプル）のリストに なるんで、それをmapに掛けて、それぞれの項を計算。最後は畳み込んで計にまとめるとな。 統計式にシグマが出てきた場合の定石だな。
> clear __show__ > cor hi low; 0.537233078618474 > cor hi pls; -0.0353964976740866
相関係数は、±1.0の範囲の数値になるんで、桁数表記をclearでデフォに戻した。それから 相関を取ってみた。最高血圧と最低血圧には正の相関有りって事だな。
pure-gplot
家でこうして遊んでいる分には、これで十分なのだけど、医者に結果を提出するには、数値を 持って行ってもしょうがない。ビジュアルじゃないとね。pureでグラフは書けるのか？
調べたら、gnuplotと連携出来るモジュールが提供されてた。パイプでgnuplotと接続して データを送り込むようだ。疎に結合してるっていいよね。unix的で。
昔おいらも、rubyのpopenを使って、gnuplotをコントロールし、pngファイルをCGIから 生成したっけなあ。その他、popenの先にgrepを接続して、よく使ったものだ。
使い方は簡単。
> using gplot;
> using namespace gplot;
> let gp = open "/usr/local/bin/gnuplot";
> plot gp hi ("style", "lines") ;
()
> output gp "png" "test.png" ;
gplot::output "png" "test.png"
> plotxy gp (hi, low) ("style", "lines") [] ;
> close gp;
いちいちコマンドを叩くって面倒。関数にしておくかな。
using gplot;
using system; // for time
let gp = gplot::open "/usr/local/bin/gnuplot";
lg xs = gplot::plot gp xs ("style", "lines") ;
// line graph with title and save png file
ltg xs ttl = ~a && b && lg xs && c with
a = gplot::title gp ttl;
s = str time + ".png";
b = gplot::output gp "png" s;
c = gplot::output gp "x11" "dmy";
end;
最初の関数は、データの確認用。後の方は、ファイルに落とす用。一応タイトルを 記入しないと、後で分からなくなるので、表題付き。
outputの切り替えが旨く行っていないため、一度しかファイルに落とせない。そして、これを 使ってしまうと、確認用の方も表示しなくなる。もう少し格闘しないと、使えるものに ならないな。ファイル名は、エポックタイムにしてある。
速くなーれ
相関係数を求める関数を定義したんだけど、結果が出るまで23秒も待たされるんだ。 気の短いオイラーには耐えられないな。
> tm n = cor (1..n) (1..n); > stats -m on
計算時間は計算サイズによって変わってくるはずなんで、簡単なベンチを書いて実行してみた。 (結果は見易いように、編集)
> tm 10; 0.001183s, 148 cells > tm 20; 0.00897s, 288 cells > tm 40; 0.059206s, 568 cells > tm 80; 0.401135s, 1128 cells > tm 160; 2.95929s, 2248 cells > tm 320; 23.3961s, 4488 cells
使用セルは、データサイズに対して比例して増えている。実行時間は、データサイズが２倍に なると、８倍の割合で増えている。O記法で書くと、O(n^3) だな。corを再掲すると
cor xs ys = (sum $ map tmc $ zip xs ys) / (#xs - 1) with
ax = avg xs ;
ay = avg ys ;
sx = sd xs ;
sy = sd ys ;
tmc z = ((z!0 - ax) / sx) * ((z!1 - ay) / sy) ;
end ;
withの中にあるのは、あくまで定義であって、評価はしていないんだな。withの前の 結果を計算してる所で、いちいち計算されてるのかな？ 検証の為に、axからsyまで 事前に計算しておくか。そしてtmcの定義もtopに出しておくか。
> let ax = 120; > let ay = 77; > let sx = 7.7; > let sy = 4.2; > tmc z = ((z!0 - ax) / sx) * ((z!1 - ay) / sy) ; > (sum $ map tmc $ zip hi low) / (#hi - 1); 0.542056399637085 0.007354s, 1804 cells
爆速になったぞ。思った通り、余計な計算をやってたんだ。どう書き換えればいいの？ 例を見てたら、withとwhenを併用してた。おいらも真似下。
coR xs ys = (sum $ map tmc $ zip xs ys) / (#xs - 1) with
tmc z = ((z!0 - ax) / sx) * ((z!1 - ay) / sy) ;
end when
ax = avg xs ;
ay = avg ys ;
sx = sd xs ;
sy = sd ys ;
end ;
早速実行
> coR hi low ; 0.537233078618474 0.102592s, 1809 cells > coR low pls; 0.0359148286634638 0.102687s, 1809 cells
withとwhenの使い所が分かったので、先にやったグラフをpngに落とす関数を修正。
ltG ttl xs = s when
a = gplot::title gp ttl;
s = str (time mod 86400) + ".png";
b = gplot::output gp "png" s;
e = lg xs ;
c = gplot::output gp "x11" "";
d = gplot::title gp "";
end;
結局、when節に書いた式は、その場で評価されるって事だ。ｃ式とｄ式は、gnuplotの状態を 元に戻す為に入れている。これで、上記の悩みは解消。思わぬ所で発見が有ったな。正月から 目出鯛いぞ。
ああ、使い方は
> lg $ map median $ group 15 hi ; () > ltG "hi 15 days priod" $ map median $ group 15 hi ; "26589L.png"
lgで、グラフの状態を確認しておき OK だったら、前の式をコンソールに呼び出して、 ちょいとタイトルを加えて実行ですかね。結果はpngファイルに送られてしまうので、表示は しないけどね。
更なる前進
とまあ、pureでも一応Pythonの統計とグラフを表示するpylabっぽいのが出来る（Rのrepl環境っぽいかな) ようになった。後は山程ある統計関数をpureに移植してくだけですかね。
上で、平気に標準偏差を求めてみたけど、そもそもデータの集合が正規分布になってるの？ そこら辺を確認しておかないといけないな。となると、次はヒストグラムの表示が欲しくなる。 どうやって、ランク分けしたらいいの？
ヒントはよそから貰ってくるのが一番。前にやった、 incanter あたりが良いですかね？ それともnumpy.histogram っつう事で、すでに自分のパソコン内に有るジャン。pythonからpureに移植しれ。 その前にどんな機能が有るか ヒストグラムの書き方２・・・Python あたりで確認かな。
まてまて、schemeチックに書いてあるRでもいいか。はて、Rのコードは入っていたかな？ それとも括弧まみれって事で、文字のヒストグラムを作ってみよう なんてのも楽しそう。
Hope
新年にあたり、今年１年、希望に満ちた年でありますようにってんで、 Hopeですよ。これ、haskellさんの親戚かな。
[sakae@pcbsd ~]$ hope >: dec fact : num -> num; >: --- fact 0 <= 1; >: --- fact n <= n*fact(n-1); >: fact 7; >> 5040 : num
>: uses lists, range; >: 1..7; >> [1, 2, 3, 4, 5, 6, 7] : list num >: product (1..7); >> 5040 : num
入れて、すぐに使えます。変数にしろ関数にしろ、ちゃんと型宣言しないと、定義させて もらえません。ruby人やPython人(pure人も)のためのHaskell人、養成ギブスと申せましょう。
fedup
新年になったので、Fedra 19 をupgradeしてみた。 Fedra Homeへ行ったら、やりかたが書いてあった。
yum install fedup fedup --network 20
19の環境で、上記を実行したら、1560個の20用rpmがDLされた。DLの一番最後がさっぱり 終わらないけど、じっくり待つ。そしてrebootしたら、実際にupgradeが適用されるのね。 何だか、毎月やってる某窓のあれみたいだな。これも永遠と待たされるけど、 お正月って事で、のんびり待ちましょう。 そんな訳で、無事に20になりました。
20になって、何か新しいのないかな。yum searchするのもいいんだけど、今回はyum list した結果をファイルに落とした。これって、バークレーDBか何かのdump結果になるのかな。 そして、そのカタログをgrepさ。こちらの方法の方がおいらに取っちゃストレス無いね。
[sakae@fedora ~]$ grep pure YUM-LIST emacs-pure.noarch 0.57-6.fc20 installed pure.i686 0.57-6.fc20 installed pure-devel.i686 0.57-6.fc20 installed pure-doc.noarch 0.57-6.fc20 installed emacs-pure-el.noarch 0.57-6.fc20 fedora ghc-pureMD5.i686 2.1.2.1-2.fc20 updates :
探し求めていたfedora用のpureがキターーだったので、早速入れておいたよ。
debianはどうする
って事で、新目のものにしてみた。
# unstable deb http://http.us.debian.org/debian/ unstable main contrib non-free deb-src http://http.us.debian.org/debian/ unstable main contrib non-free
apt-get update, apt-get dist-upgrade, apt-get autoremoveした。これで、fedora並みに 先端を走ってくれるかな。