Volo di notte

お勉強の成果メモや日常のこと

Google Optimization Toolsで組み合わせ最適化

巡回セールスマン問題やナップサック問題は、組み合わせ最適化と呼ばれる最適化問題の分野のひとつで扱われているとのこと。手法の全容は色々本が出ています。

今日から使える!組合せ最適化 離散問題ガイドブック (KS理工学専門書)

今日から使える!組合せ最適化 離散問題ガイドブック (KS理工学専門書)

実際の使い道として、生産工程の最適化などがあり、便利なソルバーとしてGoogleからGoogle Optimization Toolsが出ています。ちょっと使ってみました。

インストール

Python 3.X系の場合はpy3-ortools 6.3.4431 : Python Package Indexとして登録されているようなので…

pip install py3-ortools

でインストール。どうやら64bit系でしか動かない様なので、そこはご注意ください。

Python 2.7の場合は、ortools 6.3.4431 : Python Package Indexにあって、

pip install ortools

最新版はGithubにおいてあったりしますが、楽ちんな方法でやりました。

生産工程最適化問題

業界用語(?)で、ジョブショップ問題と呼ぶそうです。Google Optimization Toolsのページにもサンプルがあって、今回はこれを参考にしました。→The Job Shop Problem  |  Optimization  |  Google Developers. サンプルの中には若干間違いもあって、Feedback入れておきました。いつ反映されるかなー。

※ 生産工程の最適化についても制約条件や問題設定に応じて、オープンショップ問題、フローショップ問題、リスケジュール問題など色々種類があるそうです。

問題設定

こちらの問題設定をお借りして、3種類の機械(旋盤、ボール盤、フライス盤とかかな)の三種類を使う工程を3種類の部品(JOB)について、最適工程割り付けを行います。

機械と部品を便宜上、数字で呼ぶと扱いやすいので、

  1. 旋盤
  2. ボール盤
  3. フライス盤

としましょうか。工程は下記の通り(言うまでもなく、それぞれの部品の加工順序は守らなくてはなりません)

Job (部品) 工程1 工程2 工程3
部品0 旋盤 3分 ボール盤 2分 フライス盤 3分
部品1 旋盤 2分 フライス盤 3分 ボール盤 4分
部品2 ボール盤 3分 旋盤 2分 フライス盤 1分

これを、ソルバー用の書き方に書き換えると

# [m, t] : m 機械のID, t 加工時間(分)
jobs = [
    [[0,3],[1,2],[2,3]], # Job 0 / 部品0
    [[0,2],[2,3],[1,4]], # Job 1 / 部品1
    [[1,3],[0,2],[2,1]]  # Job 2 / 部品2
    ]

この問題を解かせます

結果のイメージ f:id:Chachay:20170805170527p:plain

ソルバーの宣言と問題の設定

横着して3次元の配列変数操作のためだけにnumpy使います。

# OR toolsのインポート
from ortools.constraint_solver import pywrapcp

# 変数操作を楽にするためNumpy
import numpy as np

# 視覚化用
import matplotlib.pyplot as plt

def main():
    # ソルバー定義
    solver = pywrapcp.Solver('jobshop')

    # 機械の種類
    # 0. 旋盤
    # 1. ボール盤
    # 2. フライス盤
    machines_count = 3
    all_machines = range(0, machines_count)

    # Jobの定義
    # [m, t] : m 機械のID, t 加工時間(分)
    jobs = np.array([
      [[0,3],[1,2],[2,3]], # Job 0 / 部品0
      [[0,2],[2,3],[1,4]], # Job 1 / 部品1
      [[1,3],[0,2],[2,1]]  # Job 2 / 部品2
      ])
    jobs_count = 3
    all_jobs = range(0, jobs_count)

    machines = jobs[:][:][0]
    processing_times = jobs[:][:][1]

    # 全てのJobを直列に実施した場合の時間 
    # すなわち最悪値
    horizon = 0
    for i in all_jobs:
        horizon += sum(processing_times[i])

変数と制約条件の定義

ジョブショップ問題を、制約つき離散値組み合わせ最適化問題に変換するための作業をします

    # Jobをタスクの定義にばらします
    all_tasks = {}
    for i in all_jobs:
        for j in range(0, len(machines[i])):
            all_tasks[(i, j)] = solver.FixedDurationIntervalVar(0,
                                                              horizon,
                                                              processing_times[i][j],
                                                              False,
                                                              'Job_%i_%i' % (i, j))

    # 変数の設定と離接制約(機械は1度に1つのタスク)
    all_sequences = []
    all_machines_jobs = []
    for i in all_machines:
        machines_jobs = []
        for j in all_jobs:
            for k in range(0, len(machines[j])):
                if machines[j][k] == i:
                    machines_jobs.append(all_tasks[(j, k)])
        disj = solver.DisjunctiveConstraint(machines_jobs, 'machine %i' % i)
        all_sequences.append(disj.SequenceVar())
        solver.Add(disj)

それぞれの加工順序を守るための順序制約を追加

    # 作業順序の順序制約
    for i in all_jobs:
        for j in range(0, len(machines[i]) - 1):
            solver.Add(all_tasks[(i, j + 1)].StartsAfterEnd(all_tasks[(i, j)]))

最適化の目的関数設定

全行程の一番最後のタスクが終わる時間を評価関数の目的変数に。(全てのタスクの完了時刻の最大値 = 全行程の完了時刻)

 # 目的変数定義
  obj_var = solver.Max([all_tasks[(i, len(machines[i])-1)].EndExpr()
                        for i in all_jobs])

ソルバーのオプション設定

全行程の完了時刻が最小値になるものを探します

    # 探索的に目的変数が小さいものになるように…
    sequence_phase = solver.Phase([all_sequences[i] for i in all_machines],
                                solver.SEQUENCE_DEFAULT)
    vars_phase = solver.Phase([obj_var],
                            solver.CHOOSE_FIRST_UNBOUND,
                            solver.ASSIGN_MIN_VALUE)
    main_phase = solver.Compose([sequence_phase, vars_phase])

ソルバー

最適化実行

    # 解を保存する変数
    collector = solver.LastSolutionCollector()

    collector.Add(all_sequences)
    collector.AddObjective(obj_var)

    for i in all_machines:
        sequence = all_sequences[i];
        sequence_count = sequence.Size();
        for j in range(0, sequence_count):
            t = sequence.Interval(j)
            collector.Add(t.StartExpr().Var())
            collector.Add(t.EndExpr().Var())

    solver.Solve(main_phase, [objective_monitor, collector])

結果表示

ぐちゃぐちゃしてきてしまったので、最後のソースコード全てのところで…。

ソース

from __future__ import print_function

# OR toolsのインポート
from ortools.constraint_solver import pywrapcp

# 変数操作を楽にするためNumpy
import numpy as np

# 視覚化用
import matplotlib.pyplot as plt

def main():
    # ソルバー定義
    solver = pywrapcp.Solver('jobshop')

    # 機械の種類
    # 0. 旋盤
    # 1. ボール盤
    # 2. フライス盤
    machines_count = 3
    all_machines = range(0, machines_count)

    # Jobの定義
    # [m, t] : m 機械のID, t 加工時間(分)
    jobs = np.array([
      [[0,3],[1,2],[2,3]], # Job 0 / 部品0
      [[0,2],[2,3],[1,4]], # Job 1 / 部品1
      [[1,3],[0,2],[2,1]]  # Job 2 / 部品2
      ])
    jobs_count = 3
    all_jobs = range(0, jobs_count)

    machines = jobs[:,:,0].tolist()
    processing_times = jobs[:,:,1].tolist()
   
    # 全てのJobを直列に実施した場合の時間 
    # すなわち最悪値
    horizon = 0
    for i in all_jobs:
        horizon += sum(processing_times[i])

    # Jobをタスクの定義にばらします
    all_tasks = {}
    for i in all_jobs:
        for j in range(0, len(machines[i])):
            all_tasks[(i, j)] = solver.FixedDurationIntervalVar(0,
                                                              horizon,
                                                              processing_times[i][j],
                                                              False,
                                                              'Job_%i_%i' % (i, j))

    # 変数の設定と離接制約(機械は1度に1つのタスク)
    all_sequences = []
    all_machines_jobs = []
    for i in all_machines:
        machines_jobs = []
        for j in all_jobs:
            for k in range(0, len(machines[j])):
                if machines[j][k] == i:
                    machines_jobs.append(all_tasks[(j, k)])
        disj = solver.DisjunctiveConstraint(machines_jobs, 'machine %i' % i)
        all_sequences.append(disj.SequenceVar())
        solver.Add(disj)

    # 作業順序の順序制約
    for i in all_jobs:
        for j in range(0, len(machines[i]) - 1):
            solver.Add(all_tasks[(i, j + 1)].StartsAfterEnd(all_tasks[(i, j)]))
  
    # 目的変数定義
    obj_var = solver.Max([all_tasks[(i, len(machines[i])-1)].EndExpr()
                        for i in all_jobs])
    objective_monitor = solver.Minimize(obj_var, 1)
    
    # 探索的に目的変数が小さいものになるように…
    sequence_phase = solver.Phase([all_sequences[i] for i in all_machines],
                                solver.SEQUENCE_DEFAULT)
    vars_phase = solver.Phase([obj_var],
                            solver.CHOOSE_FIRST_UNBOUND,
                            solver.ASSIGN_MIN_VALUE)
    main_phase = solver.Compose([sequence_phase, vars_phase])
    
    # 解を保存する変数
    collector = solver.LastSolutionCollector()

    collector.Add(all_sequences)
    collector.AddObjective(obj_var)

    for i in all_machines:
        sequence = all_sequences[i];
        sequence_count = sequence.Size();
        for j in range(0, sequence_count):
            t = sequence.Interval(j)
            collector.Add(t.StartExpr().Var())
            collector.Add(t.EndExpr().Var())
    disp_col_width = 10

    # 最適化計算
    if solver.Solve(main_phase, [objective_monitor, collector]):
        print("\nOptimal Schedule Length:", collector.ObjectiveValue(0), "\n")
        sol_line = ""
        sol_line_tasks = ""
        print("Optimal Schedule", "\n")

        Result = np.zeros((3,3,3))
        
        for i in all_machines:
            seq = all_sequences[i]
            sol_line += "Machine " + str(i) + ": "
            sol_line_tasks += "Machine " + str(i) + ": "
            sequence = collector.ForwardSequence(0, seq)
            seq_size = len(sequence)

            for j in range(0, seq_size):
                t = seq.Interval(sequence[j]);
                sol_line_tasks +=  t.Name() + " " * (disp_col_width - len(t.Name()))

            for j in range(0, seq_size):
                t = seq.Interval(sequence[j]);
                sol_tmp = "[" + str(collector.Value(0, t.StartExpr().Var())) + ","
                sol_tmp += str(collector.Value(0, t.EndExpr().Var())) + "] "
                sol_line += sol_tmp + " " * (disp_col_width - len(sol_tmp)) 

            sol_line += "\n"
            sol_line_tasks += "\n"

            # 結果
            for j in range(0, seq_size):
                t = seq.Interval(sequence[j]);
                r, c = int(t.Name()[4]),int(t.Name()[6])
                Result[r,c,0]=i
                Result[r,c,1]=collector.Value(0, t.StartExpr().Var())
                Result[r,c,2]=collector.Value(0, t.EndExpr().Var())
            
        print(sol_line_tasks)
        print("Time Intervals for Tasks\n")
        print(sol_line)
        
        print(Result)
        
        # 結果可視化
        lglabels = ["Machine %i"%i for i in range(0,machines_count)]
        ylabels  = ["Job %i"%i for i in range(0,jobs_count)]
        colormap = {0:'#27d67e', 1:'#2780d6', 2:'#d62728'}
        color_mapper = lambda cps:[colormap[c] for c in cps]
        temp = np.array(processing_times).cumsum(axis=1)
        machines = np.array(machines)
        
        fig, ax = plt.subplots()

        for i in range(0,Result.shape[1]):
            ax.barh(all_jobs, Result[:,i,2]-Result[:,i,1], left=Result[:,i,1], color=color_mapper(Result[:,i,0]),align='center')
       
        ax.set_yticks(all_jobs)
        ax.set_yticklabels(ylabels)
        ax.invert_yaxis()
        ax.set_xlabel('Time')
        ax.set_title('Optimized Process')
        plt.legend(lglabels)
        
        fig, ax = plt.subplots()

        for i in range(0,temp.shape[1]):
            left_ind = temp[:,i-1] if i>0 else [0,0,0]
            ax.barh(all_jobs, temp[:,i], left=left_ind, color=color_mapper(machines[:,i]),align='center')
       
        ax.set_yticks(all_jobs)
        ax.set_yticklabels(ylabels)
        ax.invert_yaxis()
        ax.set_xlabel('Time')
        ax.set_title('Job Definition')
        plt.legend(lglabels)
        plt.show()

if __name__ == '__main__':
    main()

日本取引所グループのオプション理論価格等情報のPandas読み込み~整形まで

金融工学のお勉強に公開データを使ってブラックショールズ方程式のお勉強。お勉強開始にあたってデータ準備まで。使っている技術の詳細は以下の本で。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

まずは、いつものライブラリ読み込み。

%matplotlib inline
# %matplotlib notebook
import sqlite3
import pandas as pd
from matplotlib import pylab as plt
import numpy as np

from datetime import datetime
import dateutil

あとでcsv読み込み時に使う便利関数。

def strip(text):
    try:
        return text.strip()
    except AttributeError:
        return text

データ準備

オプション理論価格等情報 | 日本取引所グループから入手した2017年2月10日の終値データ。 色々入っているけど2017年3月限の原資産:日経225のオプション取引データだけ使う。

Pandasへデータ読み込み

  1. データ本体と別に配布されているヘッダー情報をもとに自分でヘッダーを起こして名前をつける
  2. 読み込み中にテキストデータの空白をトリムする
# ヘッダー名 (変数名順)
#     商品コード, 商品タイプ, 限月, 権利行使価格, 予備
#     プットオプション: 銘柄コード, 終値, 予備, 理論価格, ボラティリティ
#     コールオプション: 銘柄コード, 終値, 予備, 理論価格, ボラティリティ
#     原資産終値, 基準ボラティリティ
colName = ("CODE","TYPE","MATURITY","STRIKE", "RSV", 
           "PUT_CODE", "PUT_PRICE", "PUT_RSV", "PUT_TPRICE", "PUT_VOLATILITY",
           "CALL_CODE","CALL_PRICE","CALL_RSV","CALL_TPRICE","CALL_VOLATILITY",
           "F225_PRICE", "Base_VOL")

df = pd.read_csv('./ose20170210tp.csv',names=colName,
                 converters = {'CODE' : strip,
                               'TYPE' : strip})
df.head()

読み込んだデータの確認。

CODE TYPE MATURITY STRIKE RSV PUT_CODE PUT_PRICE PUT_RSV PUT_TPRICE PUT_VOLATILITY CALL_CODE CALL_PRICE CALL_RSV CALL_TPRICE CALL_VOLATILITY F225_PRICE Base_VOL
0 NK225E OOP 201703 11000.0 182031018 0.0 0.0 0.94 0.674892 192031018 8350.0 0.0 8372.22 0.792450 19378.93 0.1717
1 NK225E OOP 201703 11250.0 182031218 0.0 0.0 0.99 0.652465 192031218 0.0 0.0 8122.22 0.763183 19378.93 0.1717
2 NK225E OOP 201703 11500.0 182031518 0.0 0.0 1.03 0.627990 192031518 7850.0 0.0 7872.21 0.734523 19378.93 0.1717
3 NK225E OOP 201703 11750.0 182031718 1.0 0.0 1.07 0.603704 192031718 0.0 0.0 7622.20 0.706439 19378.93 0.1717
4 NK225E OOP 201703 12000.0 182032018 0.0 0.0 1.10 0.579896 192032018 0.0 0.0 7372.20 0.678904 19378.93 0.1717

2017年3月限の日経225オプションだけ抜き出し。ついでに要らない列を削除してスッキリ。

df = df.query("MATURITY == 201703 & CODE==\"NK225E\"")\
    .drop(['RSV','PUT_RSV','CALL_RSV','PUT_CODE','CALL_CODE','CODE','TYPE','MATURITY'], 1)
  • PUTとCALLが分かれてしまっているので、データを正規化。
  • CALL列がTRUEのときCALLのデータ、FLASEのとき、PUTのデータとする。
  • 行使価格も14000円以上、22000円未満に絞る
df_p = df[["STRIKE","PUT_PRICE","PUT_TPRICE", "PUT_VOLATILITY","F225_PRICE", "Base_VOL"]]\
    .rename(columns={'PUT_PRICE': 'OP_PRICE', 'PUT_TPRICE':'OP_TPRICE', 'PUT_VOLATILITY':'OP_VOL'})
df_p['CALL'] = False
df_c = df[["STRIKE","CALL_PRICE","CALL_TPRICE", "CALL_VOLATILITY","F225_PRICE", "Base_VOL"]]\
    .rename(columns={'CALL_PRICE': 'OP_PRICE', 'CALL_TPRICE':'OP_TPRICE', 'CALL_VOLATILITY':'OP_VOL'})
df_c['CALL'] = True
df = df_p.append(df_c).query("OP_PRICE > 1.0 & STRIKE < 22000 & STRIKE >= 14000")
del (df_p,df_c)
df.head()

CALLとPUTで共通の列名に変更。

STRIKE OP_PRICE OP_TPRICE OP_VOL F225_PRICE Base_VOL CALL
17 15250.0 2.0 1.53 0.318070 19378.93 0.1717 False
18 15500.0 2.0 2.00 0.306616 19378.93 0.1717 False
19 15750.0 3.0 2.56 0.294521 19378.93 0.1717 False
20 16000.0 3.0 3.00 0.279297 19378.93 0.1717 False
21 16250.0 5.0 5.00 0.275828 19378.93 0.1717 False

オプション価格情報をPUT, CALL別にプロット

cmap = plt.get_cmap('rainbow')
fig, ax = plt.subplots()

for i, (key, group) in enumerate(df.groupby(['CALL']), start=1):
    ax = group.plot(ax=ax, kind='Scatter', x='STRIKE', y='OP_PRICE',color=cmap(i / 2.0),label = "CALL" if key else "PUT")
    
plt.title(u"MAT 2017/03 / 2017/2/10 close")
fig.patch.set_facecolor('white') 
plt.show()

f:id:Chachay:20170212223312p:plain

理論価格とか基準ボラティリティって何だろうね。

参考:

chachay.hatenablog.com

Pandas Groupbyの練習 (チートシート?)

Pandasを使っているとGroupbyな処理をしたくなることが増えてきます。ドキュメントを読んだりしながらよく使ったりする機能の骨格をまとめました。手っ取り早く勉強するなら、本が簡単そうです。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

では、さっそくスタート。

pandasライブラリなどなどの読み込み

普段使っている感じのものをコピペ。matplotlibはinlineの方が好き。

%matplotlib inline
# %matplotlib notebook
import sqlite3
import pandas as pd
from matplotlib import pylab as plt
import numpy as np

from datetime import datetime
import dateutil

csvからpandasへのデータの読み込み

Climate Data Onlineから引っ張ってきた2016年1月の東京と京都に関する気象データをちょっと加工したものを使用。アメリカが出典なだけあって、単位がアメリカ風味。

df = pd.read_csv('./201601.csv',encoding="SHIFT-JIS")
df.head()
Station YEARMODA TEMP MAXT MINT WDSP MXSPD
0 Tokyo 20160101 45.1 53.8 37.2* 4.9 9.9
1 Tokyo 20160102 45.3 56.1 36.7 2.5 6.0
2 Tokyo 20160103 49.6 61.2 37.4 3.0 5.1
3 Tokyo 20160104 49.9 59.4 39.9 3.4 6.0
4 Tokyo 20160105 51.5 59.7 41.2 3.3 5.1
  1. Station : 観測地点。今回は東京か京都。
  2. YEARMODA : YYYYMMDDフォーマットの年月日データ
  3. TEMP : 1日の平均気温。単位 華氏[℉]
  4. MAXT : 1日の最高気温。単位 華氏[℉] *印付きには事情があるのだが、今はゴミ。
  5. MINT : 1日の最低気温。単位 華氏[℉] *印付きは事情があるのだが、今はゴミ。
  6. WDSP : 1日の平均風速。単位 1/10ノット
  7. MXSP : 1日の最大風速。単位 1/10ノット
df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 60 entries, 0 to 59
Data columns (total 7 columns):
Station     60 non-null object
YEARMODA    60 non-null int64
TEMP        60 non-null float64
MAXT        60 non-null object
MINT        60 non-null object
WDSP        60 non-null float64
MXSPD       60 non-null float64
dtypes: float64(3), int64(1), object(3)
memory usage: 3.4+ KB

pandasでのdatetime型への変換; 文字列→浮動小数点変換

ちょっとした前処理をします。MAXT, MINTがたまに入っている*のせいでobject扱いに。floatに直そう。ついでにYEARMODAも変換。

# YEARMODAをdatetime型にパース
df['YEARMODA'] = df['YEARMODA'].apply(lambda x: dateutil.parser.parse(str(x)))
# *を除去してfloat型に。
df['MAXT'] = df['MAXT'].apply(lambda x: float(str.replace(x, '*', '')))
df['MINT'] = df['MINT'].apply(lambda x: float(str.replace(x, '*', '')))
df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 60 entries, 0 to 59
Data columns (total 7 columns):
Station     60 non-null object
YEARMODA    60 non-null datetime64[ns]
TEMP        60 non-null float64
MAXT        60 non-null float64
MINT        60 non-null float64
WDSP        60 non-null float64
MXSPD       60 non-null float64
dtypes: datetime64[ns](1), float64(5), object(1)
memory usage: 3.4+ KB
df.head()
Station YEARMODA TEMP MAXT MINT WDSP MXSPD
0 Tokyo 2016-01-01 45.1 53.8 37.2 4.9 9.9
1 Tokyo 2016-01-02 45.3 56.1 36.7 2.5 6.0
2 Tokyo 2016-01-03 49.6 61.2 37.4 3.0 5.1
3 Tokyo 2016-01-04 49.9 59.4 39.9 3.4 6.0
4 Tokyo 2016-01-05 51.5 59.7 41.2 3.3 5.1

Pandasの行方向計算

華氏から摂氏への変換

[C] = ([F] − 32) * 5/ 9で変換

df["TEMP"] = (df["TEMP"]-32)*5/9
df["MAXT"] = (df["MAXT"]-32)*5/9
df["MINT"] = (df["MINT"]-32)*5/9

風速:0.1ノットからm/sへの変換

0.1 ノット = 0.0514444[m/s]だそうです

df["WDSP"] *= 0.051444
df["MXSPD"] *= 0.051444
df.head()
Station YEARMODA TEMP MAXT MINT WDSP MXSPD
0 Tokyo 2016-01-01 7.277778 12.111111 2.888889 0.252076 0.509296
1 Tokyo 2016-01-02 7.388889 13.388889 2.611111 0.128610 0.308664
2 Tokyo 2016-01-03 9.777778 16.222222 3.000000 0.154332 0.262364
3 Tokyo 2016-01-04 9.944444 15.222222 4.388889 0.174910 0.308664
4 Tokyo 2016-01-05 10.833333 15.388889 5.111111 0.169765 0.262364

groupbyとグループ毎での統計量

観測地点ごとにグルーピングして統計情報を出力をする練習。東京、京都それぞれの2016年1月の平均気温の平均、2016年1月の最高気温・最低気温を出してみましょう。

個別に出力する

統計量を個別に出力する場合です。あとでまとめて出す方法やります。

平均気温の平均

df.groupby("Station")["TEMP"].mean()
Station
Kyoto    5.881226
Tokyo    6.035842
Name: TEMP, dtype: float64

最高気温

df.groupby("Station")["MAXT"].max()
Station
Kyoto    16.000000
Tokyo    16.222222
Name: MAXT, dtype: float64

最低気温

df.groupby("Station")["MINT"].min()
Station
Kyoto   -4.111111
Tokyo   -2.611111
Name: MINT, dtype: float64

参考:Pandas組み込みのGroup用統計関数

上で使った基本的なもの以外にも、Essential Basic Functionality — pandas 0.19.2 documentationにある通り、いろいろな関数が組み込まれています。抜粋して翻訳。

関数 機能
count Non-Nullなデータの数
sum 合計値
mean 平均値
mad 平均絶対偏差
median メジアン(中央値)を出力
min 最小値
max 最大値
mode モード(最頻値)
abs 絶対値
prod 総乗
std 標準偏差. 自由度を引数にとれるがデフォルトは1.
var 不偏分散
sem 平均値の標準誤差(cf : 推計統計学)
skew 歪度
kurt 尖度
quantile 分位数またはクォンタイル。デフォルトでは中央値を出力する。
.quantile(.25)で25パーセンタイルを出力。
[参考]pandas.quantile — pandas
cumsum 累積和を出力(デフォルトで行方向)
cumprod 累積乗を出力 (デフォルトで行方向)
cummax 累積最大値。axis指定が必要。
cummin 累積最小値。累積最大値と似てる。

まとめて出力する

個別に統計量を出力させてもよいのですが、まとめてバッっと出したい場合は。.agg関数を使うと素敵になる。

aggregations = {
    'TEMP': {
        '1月の平均気温の平均': 'mean'
    },
    'MAXT': {
        '1月の最高気温' : 'max'
    },
    'MINT': {
        '1月の最低気温' : 'min'
    }
}
df.groupby("Station").agg(aggregations)
TEMP MAXT MINT
1月の平均気温の平均 1月の最高気温 1月の最低気温
Station
Kyoto 5.881226 16.000000 -4.111111
Tokyo 6.035842 16.222222 -2.611111

Groupbyとmatplotlib

Pandasでグループ毎に分けてプロットする方法。連続してプロットすると重ね描きになるのはmatplotlibと同様。

cmap = plt.get_cmap('rainbow')
fig, ax = plt.subplots()

#df['YEARMODA'] = df.YEARMODA.astype(np.int64)

for i, (key, group) in enumerate(df.groupby(['Station']), start=1):
    ax = group.plot(ax=ax, x='YEARMODA', y='TEMP',color=cmap(i / 4.0),label = key+"_avg")
    
for i, (key, group) in enumerate(df.groupby(['Station']), start=1):
    ax = group.plot(ax=ax, x='YEARMODA', y='MAXT',color=cmap(2*(i+1) / 4.0),label = key+"_max")
    
plt.title(u"2016/01 Temp")
fig.patch.set_facecolor('white') 
plt.show()

f:id:Chachay:20170212151243p:plain

Groupレベルの処理をDataFrameに戻す

Group別の処理結果を使って、元のデータを処理したい時のお話。

DataFrameにGroup別最大値を追加

DataFrameに、観測地別2016年1月最高気温をMAXT_Monthという列名で追加

df['MAXT_Month'] = df.loc[df.groupby('Station')['MAXT'].transform('idxmax'), 'MAXT'].values

## 別の方法
# df = df.join(df.groupby('Station')['MAXT'].max(), on='Station', rsuffix='_Month')

df.head()
Station YEARMODA TEMP MAXT MINT WDSP MXSPD MAXT_Month
0 Tokyo 2016-01-01 7.277778 12.111111 2.888889 0.252076 0.509296 16.222222
1 Tokyo 2016-01-02 7.388889 13.388889 2.611111 0.128610 0.308664 16.222222
2 Tokyo 2016-01-03 9.777778 16.222222 3.000000 0.154332 0.262364 16.222222
3 Tokyo 2016-01-04 9.944444 15.222222 4.388889 0.174910 0.308664 16.222222
4 Tokyo 2016-01-05 10.833333 15.388889 5.111111 0.169765 0.262364 16.222222

Group別に他の列の条件に応じて代入値を決める

例えば、観測地別2016年1月 最高気温日の最大風速をMXSPD_of_DAY_MAXTという列名で追加

df['MXSPD_of_DAY_MAXT'] = df.loc[df.groupby('Station')['MAXT'].transform('idxmax'), 'MXSPD'].values
df.head()
Station YEARMODA TEMP MAXT MINT WDSP MXSPD MAXT_Month MXSPD_of_DAY_MAXT
0 Tokyo 2016-01-01 7.277778 12.111111 2.888889 0.252076 0.509296 16.222222 0.262364
1 Tokyo 2016-01-02 7.388889 13.388889 2.611111 0.128610 0.308664 16.222222 0.262364
2 Tokyo 2016-01-03 9.777778 16.222222 3.000000 0.154332 0.262364 16.222222 0.262364
3 Tokyo 2016-01-04 9.944444 15.222222 4.388889 0.174910 0.308664 16.222222 0.262364
4 Tokyo 2016-01-05 10.833333 15.388889 5.111111 0.169765 0.262364 16.222222 0.262364