PythonとShell Scriptを使えば研究がこんなにラクになる！

Excel芸人からの脱却

研究してるとこういうことありませんか？
「大量のアウトプットデータがcsvで手元にあって、それに対して色々な統計量を計算したりグラフを書きたい・・・」
みたいな状況

こういうときは多くの方はExcelを使うと思います。学習コスト低くていいですよねExcel。でも、こういう状況だとどうでしょう？

私の状況

※数値は実際とは異なります

各データディレクトリにEXPから始まるcsvデータが10個、トップディレクトリにtrue.csvというcsvがあります。

作業ディレクトリ

.
├── data1
│   ├── EXP0.csv
│   ├── EXP1.csv
│   ├── EXP2.csv
│   ├── EXP3.csv
│   ├── EXP4.csv
│   ├── EXP5.csv
│   ├── EXP6.csv
│   ├── EXP7.csv
│   ├── EXP8.csv
│   ├── EXP9.csv
├── data2
│   ├── EXP0.csv
│   ├── EXP1.csv
│   ├── EXP2.csv
│   ├── EXP3.csv
│   ├── EXP4.csv
│   ├── EXP5.csv
│   ├── EXP6.csv
│   ├── EXP7.csv
│   ├── EXP8.csv
│   ├── EXP9.csv
└── true.csv

データが取得したEXP*.csv(2,3行目はx,y座標)

data1/EXP0.csv

position,15.066053850884526,-0.9638563631665111,0.0
position,12.319995528632282,-2.4703297321884445,0.0
position,11.132448689202752,-3.712951287760663,0.0
position,13.148068482052047,-5.45596151293545,0.0
position,12.293560662818214,-5.953914266659009,0.0
position,11.601837637284774,-4.275343005129771,0.0
position,10.534288266317686,-2.3469451295032244,0.0
position,8.773279952240344,-7.19692299801392,0.0
position,6.9630858201847765,-7.043813715693978,0.0
position,5.170301830356207,-6.595090113534071,0.0
position,3.993313729379854,-5.989613016643561,0.0

正位置のx.y座標が書かれたtrue.csv

true.csv

16,-1
16,-3
16,-5
16,-7
16,-8
14,-8
12,-8
10,-8
8,-8
6,-8
4,-8

やりたいこと

各10個あるEXP*.csvのx,y座標（2,3行目）とtrue.csvのx,y座標のユークリッド距離（誤差）を計算して、最大・最小・平均・分散を計算したい
また誤差の累積度数分布（CDF）が書きたい（よく分からない人は、とりあえず複雑なグラフ書きたいと思ってもらえばOK）

これエクセルでやると

1. 各csvを開いてコピーしてエクセルに張り付けて
2. 累積度数分布計算のために必要なパラメータ計算して
3. グラフを書く
   となってかなりめんどくさいです。特にデータの数が増えるとその分だけ手間が増えます。dataディレクトリ2つならまだいいんですが実際は10近くあったり・・・

なのでShell Scriptとpythonを使ってやろう！というのが今回の趣旨

csvを結合するシェルスクリプト

**こういう出力が得たい（クリックで展開）**

1～11行目がEXP0.csvの内容 11~22行目がEXP1.csvの内容 ・・・

all.csv

15.066053850884526,-0.9638563631665111
12.319995528632282,-2.4703297321884445
11.132448689202752,-3.712951287760663
13.148068482052047,-5.45596151293545
12.293560662818214,-5.953914266659009
11.601837637284774,-4.275343005129771
10.534288266317686,-2.3469451295032244
8.773279952240344,-7.19692299801392
6.9630858201847765,-7.043813715693978
5.170301830356207,-6.595090113534071
3.993313729379854,-5.989613016643561
16.33091125786493,-2.430017019017235
15.546685399770052,-3.1439349247372985
15.358443318203225,-5.188186868656099
15.803939574162833,-6.589544184050087
15.831574705831828,-7.954407158755075
13.870771718852922,-8.169534440618603
11.886055506642025,-7.634787110875791
10.135175263759995,-7.846760636016473
8.097965638406897,-7.180413132106358
6.116121345242821,-7.283352173524063
3.845316786655085,-7.55118438806831
15.556822841978402,-1.8508877393484549
15.822037588510105,-2.7397389516421837
16.179835414491677,-4.689521141099159
16.49082918770574,-6.354200310869064
16.680677753437532,-7.63659983475943
14.912406400154248,-8.388916230343355
10.399747107539138,-5.340893630088389
8.725053107040047,-5.864709379625428
8.103743077145955,-6.290356978450172
7.86285825348623,-6.398488476116505
7.834577847872636,-6.687439598791931
15.06455973931716,-4.212147043093617
14.746360403230558,-5.355163603871441
15.398591402920815,-5.656969294130488
15.937493635436809,-7.491544706983241
15.48357680580834,-8.32461598946261
13.589864833048253,-8.164541410322329
12.036900949008915,-7.733274046268607
10.253804184787986,-7.6566817883753275
8.300576277966238,-7.417723534180651
6.36111577489307,-6.987024081307401
4.321121024486674,-6.435155968298602
15.087854945648697,-2.498613532725733
14.846401151617915,-3.7751857824725614
15.844406157088924,-4.817787667912111
13.205637167460901,-7.396355249246657
13.440792741642062,-8.676351131889774
12.71654702940811,-9.30681499035674
11.23702991749847,-8.287467814802014
9.438639947749552,-8.117284937068618
8.293389155573697,-5.077281577003624
7.942443029430977,-5.787902935380523
6.070411685025602,-3.9373767255280763
14.655815826852917,-1.1996849880845688
12.222360856186253,-3.311190549291358
12.9190566007204,-5.111427067283211
13.423860472904536,-5.917070779647318
12.343745346110385,-7.41575220920902
11.711717430729434,-5.691923563575875
10.142301400871796,-6.094553186341865
8.575202445739784,-6.871644085709193
7.9158890921348934,-6.484863141786547
7.860011086423779,-6.75949188349868
5.119068459401129,-6.796963120860472
16.555956267741752,-0.9587462618382601
13.142137752799698,-4.66481470703684
12.656200636315637,-5.8295477268967995
12.093076489433395,-7.588061719303198
12.187481454059903,-7.470298460168054
10.36671026773049,-8.05266084311158
12.148675726347726,-7.599424744167443
10.493095045580684,-7.898743631267241
9.001118129428932,-8.061092368370923
6.908384930981552,-7.360068565793828
5.055933184996052,-7.6609498063357035
15.188931306952258,-1.549914581592393
15.737465449347583,-3.348924750158851
16.13769445228539,-4.863734469782087
13.48569507081511,-8.159926349123575
13.670303222915763,-9.039865302002957
13.928416896138017,-8.052003232101343
12.316424557066398,-8.23526384509547
9.985792092718828,-7.740003267438425
8.086746917731178,-8.014759911379112
6.4542676373641745,-8.105585582551974
4.37043612806629,-8.153947216386065
15.919736985155224,-1.9185185749033469
14.882752708851573,-2.496849360080571
15.907956404890063,-4.035278772220749
15.996947544866813,-5.549889428292126
16.189197308965333,-6.552017581046362
15.068266404624673,-7.6449547640093325
13.113295980390983,-7.934761441902208
11.331801356483977,-8.031430643610546
7.876342646837259,-5.1934801839702835
8.011799726393994,-5.8543327009229795
3.99957454885494,-7.9348479838165975
16.40942289466946,-1.4490245989348673
15.27992826354875,-2.720701597400545
15.523342741517212,-4.790080320428595
15.296252563511835,-5.679540988425728
15.877155497407925,-6.465048842768543
14.073641422550306,-7.699136846239562
12.46402527257389,-7.737845379876164
10.546937137947275,-7.818199674675334
9.60297097002121,-7.284815862221581
8.562833177448296,-7.080679570781472
6.620345660242456,-6.050833200782434

まずEXP*.csvが10個バラバラになっていると処理しにくいので統合するシェルスクリプトを書きます。とはいっても1行で十分。

ls data1/EXP*.csv | xargs cat | cut -d "," -f 2,3 > data1/all.csv

・・・ってシェルに慣れてないと意味不明ですよね！説明します！

ls data1/EXP*.csv

lsは現在のディレクトリにあるファイル（＆ディレクトリ）一覧を取得するコマンドです。
たぶん最も使うコマンドと行っても過言ではないでしょうｗ

またEXP*のようにワイルドカード使って書くことで、EXPから始まるファイル全てを対象にできます

ls data1/EXP*.csv
data1/EXP0.csv  data1/EXP2.csv  data1/EXP4.csv  data1/EXP6.csv  data1/EXP8.csv
data1/EXP1.csv  data1/EXP3.csv  data1/EXP5.csv  data1/EXP7.csv  data1/EXP9.csv

| （パイプ）

| はシェル特有の記法でパイプと呼ばれます。
a | bのように書くことでaの標準出力をbの標準入力に受け渡すことが出来ます。
超便利！！

【パイプ】Linuxのコマンドで使う縦線「|」の意味と使い方 | UX MILK

xargs cat

まずcatについて説明します。cat file1とすることでファイルの内容を標準出力に出力します。またcat file1 file2とすることで「file1の最終行にfile2の内容を出力」することができます。つまり複数のファイルを行方向に連結できます（ちなみに列方向はpasteコマンド）。

次にxargsコマンドですが、これは受け取った入力をコマンドの引数にすることができます。つまり

echo "file1 file2" | xargs cat 

これは以下と同じ意味になります

cat file1 file2

xargs コマンド――コマンドラインを作成して実行する

ここまでの出力結果はこうなります。各11行あるEXP0~9までのcsvデータが行方向に連結して110行の出力が得られています。

ls data1/EXP*.csv | xargs cat
position,15.066053850884526,-0.9638563631665111,0.0
position,12.319995528632282,-2.4703297321884445,0.0
position,11.132448689202752,-3.712951287760663,0.0
position,13.148068482052047,-5.45596151293545,0.0
・・・（長いので省略）
position,8.562833177448296,-7.080679570781472,0.0
position,6.620345660242456,-6.050833200782434,0.0

cut -d "," -f 2,3

次にcutコマンドです。これは区切り文字でフィールドを分けてくれるコマンド。

• -d で区切り文字を指定
• -f でどのフィールドを取得するか（注：1から始まる）

今回は2列目と3列目が欲しいので -f 2,3 の指定でいいですね。
これでここまで加工できました！

ls data1/EXP*.csv | xargs cat | cut -d "," -f 2,3
15.066053850884526,-0.9638563631665111
12.319995528632282,-2.4703297321884445
11.132448689202752,-3.712951287760663
13.148068482052047,-5.45596151293545
・・・（長いので省略）
8.562833177448296,-7.080679570781472
6.620345660242456,-6.050833200782434

> data1/all.csv

ここまでで画面上に出力できましたが、これをファイルに保存したいですよね？
> 書き出したいファイル名 とすれば画面上に出力したデータをファイルにデータを書き出せます！（上書き処理。無い場合は作成）
ちなみに>>とするとファイルの更新です。

ただし、これだと画面上に出力されなくなるので、書き込む内容を確認しながらファイルに書き込みたい場合は| tee ファイル名とすると画面に書き出しつつ同時にファイルに保存できます。（tee -aとすると上書きを防げる）

• リダイレクトの使い方
• 【 tee 】 標準入力を標準出力とファイルに出力する | 日経 xTECH（クロステック）

csvからグラフや統計量を計算

統計量の計算

さて、ここまででcsvが出来ました！
ここから平均や分散を計算するのはawkコマンドでも出来るのですが、さすがに大変になってくるのでpythonの出番です！

このように

• テキスト(csvなど)処理はシェルスクリプト
• 実際の計算はpython

というように「言語の得意分野」に合わせて分けるのがコツかなと思います。
統計処理はこんな感じで書きました。numpyライブラリ使うのがコツですね！

# coding:utf-8
import numpy as np
import sys
from collections import OrderedDict

def print_statistics(dir_name, errors, d_num, r_num):
    print("-----" + dir_name + " all_xy" + "-----")
    print("Average error")
    print(np.mean(errors))
    print("Variance")
    print(np.var(errors))
    print("Max error")
    print(np.max(errors))
    print("Min error")
    print(np.min(errors))
    errors_sorted = np.sort(errors)
    # python3のroundは端数が0.5の時に偶数方向へ丸めるので注意 例)22.5 -> 22
    percent90 = round((d_num*r_num/10)*9) - 1
    print("90% error")
    print(errors_sorted[percent90])

# @Output
# {PNUM100: [
#   [x0,y0],
#   [x1,y1],
#   ...
# ],
# PMUM500: [
# ...
# ],
# ...
# }
def load_all_csvs(args):
    all_xy_s = OrderedDict() # note: 辞書順のループは入力した順番を必ずしも保持しないのでOrderedDictを使う
    for i in range(len(args)):
        if i==0:
            continue # args[0]は自分のファイル名
        all_xy_s[args[i]]  = np.loadtxt(args[i] + "/all.csv", delimiter = ",")
    return all_xy_s

all_xy_s = load_all_csvs(args=sys.argv)
true_xy = np.loadtxt("true.csv", delimiter = ",")

d_num = 10 # データ数
r_num = 11 # 行数
true_xy = [true_xy[i%r_num] for i in range(len(true_xy)*d_num)] # change len(true_xy) into 110 for all.csv
for dir_name, all_xy in all_xy_s.items():
    errors = [np.linalg.norm(true_xy[i] - all_xy[i]) for i in range(len(all_xy))]
    print_statistics(dir_name=dir_name, errors=errors, d_num=d_num, r_num=r_num)

※ 説明のため動けばいいというコードになっていますが、適宜普通に実行するコードはmainで囲むなどすべきです

グラフの描画

また、csvからグラフを書くのはpythonのmatplotlibが便利です！
特に公式のsampleが優秀で、作成できるグラフを一覧で表示してくれています。グラフをクリックするとサンプルコードまで出てくるので本当に便利・・・！！
（今回は少し難しいCDFというグラフを書こうとしたので stack overflowのお世話になりましたが笑）

Thumbnail gallery — Matplotlib 2.0.2 documentation

# coding:utf-8
from collections import OrderedDict
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
import sys

def load_all_csvs(args):
    all_xy_s = OrderedDict() # note: 辞書順のループは入力した順番を必ずしも保持しないのでOrderedDictを使う
    for i in range(len(args)):
        if i==0:
            continue # args[0]は自分のファイル名
        all_xy_s[args[i]]  = np.loadtxt(args[i] + "/all.csv", delimiter = ",")
    return all_xy_s

args = sys.argv

# データと真の値を読み込み
all_xy_s = load_all_csvs(args=sys.argv)
true_xy = np.loadtxt("true.csv", delimiter = ",")

d_num = 10 # データ数
r_num = 11 # 行数
true_xy = [true_xy[i%r_num] for i in range(len(true_xy)*d_num)] # change len(true_xy) into 110 for all.csv

for dir_name, all_xy in all_xy_s.items():
    errors = [np.linalg.norm(true_xy[i] - all_xy[i]) for i in range(len(all_xy))] # ユークリッド距離の計算

    errors_sorted = np.sort(errors)
    percent90 = (d_num*r_num//10)*9 - 1 
    print("-----" + dir_name + "-----")
    print("90% error")
    print(errors_sorted[percent90])

    #0から8まで80個にすれば0.1刻みになる
    hist_fig, hist_ax = plt.subplots()
    y, bins, patches = hist_ax.hist(errors, bins=80, range=(0,8), density=True, cumulative=True, histtype="step", label="hist")
    plt.close(hist_fig) # y, bins があればヒストグラム自体は必要ないのでここで消す(もう少しいいやり方がありそう)

    # 散布図の描画
    plt.plot(bins[:-1], y)
    plt.scatter(bins[:-1], y, s=10, label=dir_name)
# plt.title("title")
plt.xlabel("Error [m]")
plt.ylabel("Cumulative frequency [%]")
plt.legend()
plt.show()

python cdf_plot.py data1 data2 とすると、こんなグラフが！！

これでdataが増えてもコマンド一発でいつでもグラフ描ける環境が手に入りました！

GitHubに置いてます！

今回のコードはGitHubに置いているので参考にしてみてください！（とりあえずサクッと確認するためのコードなので保守性とかは考慮してないですので、実際に使う場合には修正した方がいいと思います。とはいえ今思うと流石にmainで囲むくらいすればよかった気もしますね（　＾ω＾）・・・）
記事では触れませんでしたが（コマンドをいちいち覚えておくのが大変なので）シェルスクリプトにしてすぐに実行できるように工夫されています。
GitHub - koyo-miyamura/research_csv_analyze: 研究用に作成した測位csvデータからCDFや統計データを出力するスクリプト

シェルスクリプトは面白い！

シェルスクリプトには他にもgrepやawkやsedなどが便利なので、機会があれば勉強してみると面白いと思います！

シェルスクリプトについてもっと知りたい方は、筆者がシェルスクリプト勉強するきっかけにもなった
まんがでわかるLinux　シス管系女子
がオススメです！

追記

7/29更新のQiita公式デイリーランキング載りました！

【毎日自動更新】Qiitaのデイリーストックランキング！ウィークリーもあるよ