静かなる名辞

pythonとプログラミングのこと


【python】sliceのちょっと深イイ(かもしれない)話

 リスト(じゃなくてもだけど)に次のようにアクセスするとき、内部的には__getitem__が呼ばれていることは、歴戦のpythonistaの皆さんには常識でしょう。

>>> lst = [1,2,3,4,5]
>>> lst[0]
1

 この様子を自作クラスで観察してみましょう。

>>> class Hoge:
...     def __getitem__(self, k):
...         print("__getitem__!")
...         return k
... 
>>> h = Hoge()
>>> h[0]
__getitem__!
0
>>> h["hogehoge~"]
__getitem__!
'hogehoge~'

 Hogeクラスは__getitem__が呼ばれると、「__getitem__!」とprintしてから__getitem__の引数をそのまま返します。上の結果から、実際に__getitem__が呼ばれていることがわかります。

 []は__getitem__の糖衣構文と言っても、まあ良いでしょう*1

 これで「ほーん、そうか」と納得しかけてしまいますが、「ちょっと待て、じゃあスライスはどうなるんだ・・・?」というのが疑問として浮かんできますね。

lst[0:5]

 一体何が渡るというんでしょう。「0:5」なんてオブジェクトはありませんから(そのまま書いても構文エラー)、スライスの場合は__getitem__とは違う仕組みで処理されるのでしょうか?

 そうはなっていません。

>>> h[0:5]
__getitem__!
slice(0, 5, None)
>>> type(h[0:5])
__getitem__!
<class 'slice'>

 おおお、sliceオブジェクトなんていうのが渡っている・・・。

 このsliceオブジェクトは、普通にpythonを書いている限り目にする機会はほとんどないと思います。

 それでも、公式ドキュメントにはしっかり載っています。

class slice(start, stop[, step])
range(start, stop, step) で指定されるインデクスの集合を表す、スライス (slice) オブジェクトを返します。引数 start および step はデフォルトでは None です。スライスオブジェクトは読み出し専用の属性 start、stop および step を持ち、これらは単に引数で使われた 値 (またはデフォルト値) を返します。これらの値には、その他のはっきりと した機能はありません。しかしながら、これらの値は Numerical Python および、その他のサードパーティによる拡張で利用されています。スライスオブジェクトは拡張されたインデクス指定構文が使われる際にも生成されます。例えば a[start:stop:step] や a[start:stop, i] です。この関数の代替となるイテレータを返す関数、itertools.islice() も参照してください。

2. 組み込み関数 — Python 3.6.5 ドキュメント

 なるほど~、こいつが渡ることで、あとは受けるオブジェクトの__getitem__が然るべき処理をしてくれればスライスが実現する仕組みになっているんですね。よくできてる・・・。

「いや、ちょっと待て。numpyで使うアレはどうなってるんだ」

 こういう奴のことですね。

a[:,0]

 上記sliceオブジェクトはstart, stop, stepしか持たないので、こういうスライスは手に負えなさそうです。

 なんてこった、今度こそ__getitem__では処理しきれなくて、なにか違う仕組みで処理されているのか・・・

 いません。

>>> h[:,0]
__getitem__!
(slice(None, None, None), 0)
>>> type(h[:,0])
__getitem__!
<class 'tuple'>

 __getitem__なのは間違いないみたいですが、 なぜかtupleが返ります。0要素目は空のslice, 二番目は入力がそのまま・・・? 一体どうなっているんだ?

 実はこうなっています。

>>> h[::,0:1:2,0:-1:-2,0,1,2,3]
__getitem__!
(slice(None, None, None), slice(0, 1, 2), slice(0, -1, -2), 0, 1, 2, 3)

 なるほど、カンマで区切られたものごとにtupleの要素になっているのか!

 ……と、書きましたが一応ちゃんと説明しておくと、そもそも「pythonのtuple」はカンマで区切られた要素によって成立する構文です。関数の引数リストなど、他の構文と被る場合はそちらが優先的に扱われますが。

>>> 1,2,3
(1, 2, 3)

 これについてはドキュメントに説明があります。

4. 組み込み型 — Python 3.6.5 ドキュメント


 要するに、単に添字の中にtupleを書いているだけ、とみなしてもまあ良いでしょう*2

 わかっちゃえばなんてことはないですね。とても自然な仕様です。あとは受け取った側で然るべき処理をするだけ。

 このように、スライスの裏ではsliceオブジェクトが暗躍しています。なんとなくこういうことを知っていると深イイと思えますね。また、なんとなくスライスの構文がよくわからなかった人も、この仕様が理解できれば自然に複雑なスライスを書けるようになることでしょう(その結果、難読コードが量産されてしまうかもしれないが・・・)。

余談

>>> lst = [0,1,2]
>>> lst[:,0]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: list indices must be integers or slices, not tuple

 上のことがわかっていれば理解できますが、そうでないと絶対理解不能なエラーメッセージです(tupleなんかないじゃんかって)。

追記:2018/06/21

 この話はドキュメントのどこに載っているんだろうとずっと思っていましたが、見つけました。

6. 式 (expression) — Python 3.6.5 ドキュメント

スライス表記に対する意味付けは、以下のようになります。プライマリの値評価結果は、以下に述べるようにしてスライスリストから生成されたキーによって (通常の添字表記と同じ __getitem__() メソッドを使って) インデクス指定できなければなりません。スライスリストに一つ以上のカンマが含まれている場合、キーは各スライス要素を値変換したものからなるタプルになります; それ以外の場合、単一のスライス要素自体を値変換したものがキーになります。一個の式であるスライス要素は、その式に変換されます。適切なスライスは、スライスオブジェクト (標準型の階層 参照) に変換され、その start, stop および step 属性は、それぞれ指定した下境界、上境界、およびとび幅 (stride) になります。式がない場所は None で置き換えられます。

 説明の内容はこの記事と基本的に同じですが、グダグダなこの記事と比べると、さすがによくまとまっています。公式は偉い。

*1:細かいことを言うと、lst[0] = "hoge"など代入する場合、del dct[0]などdel文を呼ぶ場合などは__getitem__とは異なったメソッドが呼ばれます。あくまでも単純に値を参照する場合の話です

*2:sliceオブジェクトに変換されるので単純な処理ではないが

複数の目的変数で回帰を行う方法

はじめに

 回帰分析を行う際、複数の目的変数に対して回帰をしたい場合があります。普通のモデルではできないのでちょっと面食らいますが、やり方は色々あるようです。

 目次

スポンサーリンク


目的変数の数だけ回帰モデルを作る方法

 単純に考えると、一つの目的変数を出力する回帰モデルを目的変数の数だけ用意してやれば、所要を達しそうです。

 python+sklearnを使えば、これに対応したモデルが最初から用意されています。

sklearn.multioutput.MultiOutputRegressor — scikit-learn 0.20.2 documentation

 コンストラクタには好きな回帰モデルを渡してあげることができます。それが目的変数の数だけコピーされ、内部で束ねられて回帰に使われます*1

複数の目的変数に対応したモデルを使う

 上の方法は単純ですが、回帰モデルの中には自然に複数の出力に対応しているものもあります。

 そういったモデルを使うことにどんなメリットがあるのか? というと、まず目的変数の数だけ回帰モデルを作るのと比べて無駄が減るので、計算コストがケチれる可能性があります(あくまでも「可能性」の話)。

 また、複数存在する目的変数の間に何らかの相関性があれば、それも踏まえて上手く学習することでモデルの性能が上がる可能性があります(こちらもあくまでも「可能性」)。

 そういった複数の目的変数に対応したモデルを幾つか紹介します。すべては網羅しきれないので、その点はご承知ください。

正準相関分析

 正準相関分析はこの手の話で出てくる代表的なモデルです。単純な手法ですが、けっこう奥深いといえば奥深いです。

 参考(過去に書いた記事):【python】正準相関分析(Canonical Correlation Analysis)を試してみる - 静かなる名辞

 これの良いところは、説明変数と目的変数*2のそれぞれでPCAみたく新たな軸を張り、次元削減を行ってくれることです。説明変数ン百次元、目的変数20次元みたいなケースだったとしても、次元削減の効果で「わかりやすい」結果が得られる可能性があります。つまり、現象を説明するモデルとしては非常に適しています。

 欠点は、回帰モデルとして考えると性能が高いと言えるかは微妙なこと、非線形への対応は基本的にはないことです。カーネルPCAみたくカーネル法で非線形対応させたモデルもありますが、良さげなライブラリが見当たらないのと、そこまでするなら他の手法を使いたいという気持ちがあるので紹介しません。

 sklearnのモデルはこれです。上に書いた通りカーネル正準相関の実装はありません。

sklearn.cross_decomposition.CCA — scikit-learn 0.20.2 documentation

 predictメソッドでXからYを予測できるので、普通に回帰に使えます。 

 入出力が割と線形なデータで、「説明」を重視したいときは使えると思います。

ランダムフォレスト回帰

 なぜかランダムフォレスト回帰は複数出力に対応しています。解説論文を見つけたので貼っておきます。興味のある方はどうぞ(私は読んでいません)。

 とにかく使いたければsklearnのRandomForestRegressorはそのまま使えます。目的変数も説明変数と同様に配列で入れてあげてください。

3.2.4.3.2. sklearn.ensemble.RandomForestRegressor — scikit-learn 0.20.2 documentation

多層パーセプトロン(ニューラルネットワーク回帰)

 ニューラルネットですからできて当然。複数出力にするためにやることといったら出力層ユニット数を増やすだけですから、一番シンプルかもしれません。これもsklearnのがそのまま使えます。

sklearn.neural_network.MLPRegressor — scikit-learn 0.20.2 documentation

まとめ

 複数の目的変数に対して回帰を行う場合について、2種類の方法を説明しました。

  • 単純に目的変数の数だけ回帰モデルを用意する方法
  • 複数の目的変数を出力できるモデルを用いる方法

 どちらが良いかは一概には言えません。データや目的に応じて、あるいは実際に走らせてみて評価指標や計算コストを勘案して考える必要があります。複数の目的変数に最初から対応したモデルの方が良いような気もしますが、そうとも言えないんじゃという話もあったりします。

 でもまあ、色々な選択肢があることは良いことです。いろいろ勘案して選べば良いでしょう。適当ですがこんな感じでシメます。

*1:ということだと思う・・実装読んでいないので断言しかねます

*2:回帰の記事なのでそう呼ぶが、ぶっちゃけ妥当ではない。CCAの枠組みではどっちがどっちでも大して構わないのだし

【python】リストの各要素に違う処理をする

問題設定

 想定しているのは、たとえばこんなシチュエーションです。

s = "hoge! 1234"
tmp = s.split()
lst = [tmp[0], int(tmp[1])]

 要するに、比較的短いリストだが性質の違うものが入っており、それぞれ違う処理をして返したいのです。

 それだけなら良いのですが、上の例だとs.split()の結果を一時変数に入れないとどうしようもないので(二回呼ぶと再計算されてしまう)、まどろっこしいことになります。ここでリスト内包を使おうとしても、まともな処理は書けません(ん、zipでlambdaと一緒に回せばできなくはないか)。

 (私が思いつく)解決策は2つあります。どちらも関数を使います。

解決策1:愚直に関数定義

 関数の引数に渡せばローカル変数として束縛されるので、できるという考え方です。

def f(x):
    return [x[0], int(x[1])]
s = "hoge! 1234"
lst = f(s.split())

 もういっそsを渡す関数にした方が素直ですが、趣旨からずれるのでNG。

 この方法は記述が短くならないという欠点があります。

解決策2:lambda

 lispのletと同じテクニックで一時変数を束縛しています。

s = "hoge! 1234"
lst = (lambda t:[t[0], int(t[1])])(s.split())

 とても短く書け、簡潔な記述です。可読性は人によって異なると思いますが、個人的には読みづらいと思います。

まとめ

 このテクニックはリスト内包表記の中でどうしても一時変数が必要になったときに使えます。覚えておいて損はないです。

【python】calendarモジュールの使い方

 calendarモジュールは標準ライブラリに入っていて、曜日や日付の計算にはけっこう便利なモジュールらしいです。

 でもあまり周知されていないので、使い方を(自分用に)メモっておきます。

 ドキュメントはここです。
8.2. calendar — 一般的なカレンダーに関する関数群 — Python 3.6.5 ドキュメント

 目次


スポンサーリンク



introduction

>>> import calendar
>>> print(calendar.month(2018, 4))
     April 2018
Mo Tu We Th Fr Sa Su
                   1
 2  3  4  5  6  7  8
 9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30

 テキストのカレンダーが出てきました。ちょっとほっこりします。

 日曜から週が始まるようにもできます。

>>> calendar.setfirstweekday(calendar.SUNDAY)
>>> print(calendar.month(2018, 4))
     April 2018
Su Mo Tu We Th Fr Sa
 1  2  3  4  5  6  7
 8  9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30

 オブジェクト指向でも関数型でもない、手続き型パラダイムって感じです>setfirstweekday()。ちょっとおっかない。このへんに依存する処理をしたいときは、import先やimport元で変なことにならないよう、一々確認してあげる必要があるということであります。

 まあ、ちゃんとクラスインスタンスを作れるので、そっちを使うようにすれば良いのですが。

calendarモジュールのクラス

 以下のクラスがあります。

  • Calendar
  • TextCalendar
  • HTMLCalendar
  • LocaleTextCalendar
  • LocaleHTMLCalendar

 なんとなく酷い気もしますが、大目に見ましょう。HTMLはどうでも良いので(いや、使いたいって人もいるだろうけど)、TextCalendarでintroductionと同じことをやってみます。

>>> from calendar import TextCalendar
>>> tcalendar = TextCalendar(firstweekday=6)
>>> print(tcalendar.formatmonth(2018, 4))
     April 2018
Su Mo Tu We Th Fr Sa
 1  2  3  4  5  6  7
 8  9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30

 firstweekday=6は、0が月曜日で6が日曜日という仕様なので、こうしております。

 ロケールも一応やってみます。

>>> from calendar import LocaleTextCalendar
>>> ltc = LocaleTextCalendar(firstweekday=6, locale="ja_JP.UTF-8")
>>> print(ltc.formatmonth(2018, 4, w=3))
     April 2018
Su Mo Tu We Th Fr Sa
 1  2  3  4  5  6  7
 8  9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30

>>> ltc = LocaleTextCalendar(firstweekday=6, locale="ja_JP.UTF-8")
>>> print(ltc.formatmonth(2018, 4, w=5))
          42018
 日   月   火   水   木   金   土
  1   2   3   4   5   6   7
  8   9  10  11  12  13  14
 15  16  17  18  19  20  21
 22  23  24  25  26  27  28
 29  30

 フォントの問題があるので、何をやってもあまり綺麗には見えません。半角英数と全角文字を同じ文字幅で表示するフォントがあれば、綺麗に見えることでしょう。

機能

 さて、存在するクラスのうち、

  • TextCalendar
  • HTMLCalendar
  • LocaleTextCalendar
  • LocaleHTMLCalendar

 これらはCalendarのサブクラスです。出力形式をformatするだけに存在しています。

 曜日や日付の計算に使いたいと思う機能は、

  • Calendar

 クラスに集約されています。

 そしてここには大したメソッド数はありません。なので、一つずつ紹介していきましょう。

iterweekdays()

曜日の数字を一週間分生成するイテレータを返します。イテレータから得られる最初の数字は firstweekday が返す数字と同じになります。

>>> from calendar import Calendar
>>> cl = Calendar()
>>> cl.iterweekdays()
<generator object Calendar.iterweekdays at 0x7fca5a3d3ca8>
>>> list(cl.iterweekdays())
[0, 1, 2, 3, 4, 5, 6]

 これだけ。

itermonthdates(year, month)

year 年 month (1–12) 月に対するイテレータを返します。 このイテレータはその月の全ての日、およびその月が始まる前の日とその月が終わった後の日のうち、週の欠けを埋めるために必要な日を (datetime.date オブジェクトとして) 返します。

 指定した年月の日を返しますが、週の中途半端なところから始まったり、中途半端なところで終わったりすると、前後の月も週が中途半端じゃなくなる範囲まで出してくれます。長いので整形した出力を見せます。

>>> cl = Calendar(firstweekday=6)
>>> list(cl.itermonthdates(2018, 4))
[datetime.date(2018, 4, 1), datetime.date(2018, 4, 2), datetime.date(2018, 4, 3), 
datetime.date(2018, 4, 4), datetime.date(2018, 4, 5), datetime.date(2018, 4, 6),
datetime.date(2018, 4, 7), datetime.date(2018, 4, 8), datetime.date(2018, 4, 9), 
datetime.date(2018, 4, 10), datetime.date(2018, 4, 11), datetime.date(2018, 4, 12), 
datetime.date(2018, 4, 13), datetime.date(2018, 4, 14), datetime.date(2018, 4, 15), 
datetime.date(2018, 4, 16), datetime.date(2018, 4, 17), datetime.date(2018, 4, 18), 
datetime.date(2018, 4, 19), datetime.date(2018, 4, 20), datetime.date(2018, 4, 21), 
datetime.date(2018, 4, 22), datetime.date(2018, 4, 23), datetime.date(2018, 4, 24), 
datetime.date(2018, 4, 25), datetime.date(2018, 4, 26), datetime.date(2018, 4, 27), 
datetime.date(2018, 4, 28), datetime.date(2018, 4, 29), datetime.date(2018, 4, 30), 
datetime.date(2018, 5, 1), datetime.date(2018, 5, 2), datetime.date(2018, 5, 3), 
datetime.date(2018, 5, 4), datetime.date(2018, 5, 5)]

 あんまり嬉しくないかも・・・。

itermonthdays2(year, month)

year 年 month 月に対する itermonthdates() と同じようなイテレータを返します。生成されるのは日付の数字と曜日を表す数字のタプルです。

 上とほぼ同じ。返り値の型だけ違います。

>>> list(cl.itermonthdays2(2018, 4))
[(1, 6), (2, 0), (3, 1), (4, 2), (5, 3), (6, 4), (7, 5), (8, 6), (9, 0), (10, 1), 
(11, 2), (12, 3), (13, 4), (14, 5), (15, 6), (16, 0), (17, 1), (18, 2), (19, 3), (20, 4), 
(21, 5), (22, 6), (23, 0), (24, 1), (25, 2), (26, 3), (27, 4), (28, 5), (29, 6), (30, 0), 
(0, 1), (0, 2), (0, 3), (0, 4), (0, 5)]

 前後月の日付は0で返されるようです。タプルの二番目の要素は例の曜日を表す数字。

 これを処理すると簡単そうで良いですね。

itermonthdays(year, month)

year 年 month 月に対する itermonthdates() と同じようなイテレータを返します。生成されるのは日付の数字だけです。

 なんで同じようなメソッドがいくつもあるんだろうか。

>>> list(cl.itermonthdays(2018, 4))
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 0, 0, 0, 0, 0]

 (命名が)投げやりだなー(棒)。

monthdatescalendar(year, month)

year 年 month 月の週のリストを返します。週は全て七つの datetime.date オブジェクトからなるリストです。

 週のリストを返すんだって。メソッド名からは想像できない機能で、ちょっとびっくりしています。

>>> list(cl.monthdatescalendar(2018, 4))
[[datetime.date(2018, 4, 1), datetime.date(2018, 4, 2), datetime.date(2018, 4, 3), datetime.date(2018, 4, 4), datetime.date(2018, 4, 5), datetime.date(2018, 4, 6), datetime.date(2018, 4, 7)],
[datetime.date(2018, 4, 8), datetime.date(2018, 4, 9), datetime.date(2018, 4, 10), datetime.date(2018, 4, 11), datetime.date(2018, 4, 12), datetime.date(2018, 4, 13), datetime.date(2018, 4, 14)],
[datetime.date(2018, 4, 15), datetime.date(2018, 4, 16), datetime.date(2018, 4, 17), datetime.date(2018, 4, 18), datetime.date(2018, 4, 19), datetime.date(2018, 4, 20), datetime.date(2018, 4, 21)],
[datetime.date(2018, 4, 22), datetime.date(2018, 4, 23), datetime.date(2018, 4, 24), datetime.date(2018, 4, 25), datetime.date(2018, 4, 26), datetime.date(2018, 4, 27), datetime.date(2018, 4, 28)],
[datetime.date(2018, 4, 29), datetime.date(2018, 4, 30), datetime.date(2018, 5, 1), datetime.date(2018, 5, 2), datetime.date(2018, 5, 3), datetime.date(2018, 5, 4), datetime.date(2018, 5, 5)]]

monthdays2calendar(year, month)

year 年 month 月の週のリストを返します。週は全て七つの日付の数字と曜日を表す数字のタプルからなるリストです。

>>> list(cl.monthdays2calendar(2018, 4))
[[(1, 6), (2, 0), (3, 1), (4, 2), (5, 3), (6, 4), (7, 5)], 
[(8, 6), (9, 0), (10, 1), (11, 2), (12, 3), (13, 4), (14, 5)], 
[(15, 6), (16, 0), (17, 1), (18, 2), (19, 3), (20, 4), (21, 5)], 
[(22, 6), (23, 0), (24, 1), (25, 2), (26, 3), (27, 4), (28, 5)], 
[(29, 6), (30, 0), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5)]]

 ・・・説明、要る?

monthdayscalendar(year, month)

year 年 month 月の週のリストを返します。週は全て七つの日付の数字からなるリストです。

>>> list(cl.monthdayscalendar(2018, 4))
[[1, 2, 3, 4, 5, 6, 7], 
[8, 9, 10, 11, 12, 13, 14], 
[15, 16, 17, 18, 19, 20, 21], 
[22, 23, 24, 25, 26, 27, 28], 
[29, 30, 0, 0, 0, 0, 0]]

 このモジュールの世界観は理解してしまえばとてもわかりやすいので、そういう意味では良いと思います。

yeardatescalendar(year, width=3)

指定された年のデータを整形に向く形で返します。返される値は月の並びのリストです。月の並びは最大で width ヶ月(デフォルトは3ヶ月)分です。各月は4ないし6週からなり、各週は1ないし7日からなります。各日は datetime.date オブジェクトです。

 以下のメソッドの出力は長いので省略。どんなものが返るのかは、実行すればすぐ理解できます。

yeardays2calendar(year, width=3)

指定された年のデータを整形に向く形で返します (yeardatescalendar() と同様です)。週のリストの中が日付の数字と曜日の数字のタプルになります。月の範囲外の部分の日付はゼロです。

yeardayscalendar(year, width=3)

指定された年のデータを整形に向く形で返します (yeardatescalendar() と同様です)。週のリストの中が日付の数字になります。月の範囲外の日付はゼロです。

 あー、疲れた。モジュールのリファレンスを書き写すような真似してもあまり意味なかったですね。

便利なusage

 そのうち(思いついたら)書きます。第n何曜日の計算とかに使えると思います。

まとめ

 ちょっと残念だけど、たまに使えるシチュエーションがありそうではある。

【python】pandasでデータを標準得点(z得点)に変換

 データの正規化(標準化)をpandasでもやってみる。

 正規化、標準化とは、データを分散1、平均0に変換する操作である。

スポンサーリンク




 自分で書いてもできるが、scipyの関数を使うと簡単にできる。

>>> import pandas as pd
>>> df = pd.DataFrame([[1,2,3,4,5,6],
                       [6,5,4,3,2,1],
                       [0,1,2,3,4,5],
                       [5,4,3,2,1,0]], columns=[*"ABCDEF"])
>>> df.apply(stats.zscore, axis=0)
          A         B         C         D         E         F
0 -0.784465 -0.632456  0.000000  1.414214  1.264911  1.176697
1  1.176697  1.264911  1.414214  0.000000 -0.632456 -0.784465
2 -1.176697 -1.264911 -1.414214  0.000000  0.632456  0.784465
3  0.784465  0.632456  0.000000 -1.414214 -1.264911 -1.176697
>>> df.apply(stats.zscore, axis=1)
         A        B        C        D        E        F
0 -1.46385 -0.87831 -0.29277  0.29277  0.87831  1.46385
1  1.46385  0.87831  0.29277 -0.29277 -0.87831 -1.46385
2 -1.46385 -0.87831 -0.29277  0.29277  0.87831  1.46385
3  1.46385  0.87831  0.29277 -0.29277 -0.87831 -1.46385

 axis=0だと列で計算した標準得点、axis=1で行で計算した標準得点になる。

【python】pandasでDataFrameの平均と標準偏差を計算する方法

列の平均と標準偏差を計算したい

 とても簡単にできます。

>>> import pandas as pd
>>> df = pd.DataFrame([[1,2,3,4,5,6],
                       [6,5,4,3,2,1],
                       [0,1,2,3,4,5],
                       [5,4,3,2,1,0]], columns=[*"ABCDEF"])
>>> df.mean()
A    3.0
B    3.0
C    3.0
D    3.0
E    3.0
F    3.0
dtype: float64
>>> df.std()
A    2.943920
B    1.825742
C    0.816497
D    0.816497
E    1.825742
F    2.943920
dtype: float64

 何も考える必要はないのだった。

 リファレンス:
pandas.DataFrame.mean — pandas 0.24.2 documentation
pandas.DataFrame.std — pandas 0.24.2 documentation

行の平均と標準偏差を計算したい

 「転置しとけば?」という天の声が聞こえたのを無視してやります。

 numpy配列のようにaxisを指定するだけなのでこれも簡単です。

>>> import pandas as pd
>>> df = pd.DataFrame([[1,2,3,4,5,6],
                       [6,5,4,3,2,1],
                       [0,1,2,3,4,5],
                       [5,4,3,2,1,0]], columns=[*"ABCDEF"])
>>> df.mean(axis=1)
0    3.5
1    3.5
2    2.5
3    2.5
dtype: float64
>>> df.std(axis=1)
0    1.870829
1    1.870829
2    1.870829
3    1.870829
dtype: float64

 よくできてますね。

特定の列・行だけ取り出してから計算する

 基本的なindexing操作と組み合わせて使うことで、特定の行・列だけに対して計算するということも可能です。

 A, Bに対してのみ出力させたい場合。

>>> df[["A", "B"]].mean()
A    3.0
B    3.0
dtype: float64

describeメソッドで全体の雰囲気を掴む

 describeメソッドを使うと様々な統計量を勝手に出してくれます。

>>> df.describe()  # 列ごとに
             A         B         C         D         E        F
count  4.00000  4.000000  4.000000  4.000000  4.000000  4.00000
mean   3.00000  3.000000  3.000000  3.000000  3.000000  3.00000
std    2.94392  1.825742  0.816497  0.816497  1.825742  2.94392
min    0.00000  1.000000  2.000000  2.000000  1.000000  0.00000
25%    0.75000  1.750000  2.750000  2.750000  1.750000  0.75000
50%    3.00000  3.000000  3.000000  3.000000  3.000000  3.00000
75%    5.25000  4.250000  3.250000  3.250000  4.250000  5.25000
max    6.00000  5.000000  4.000000  4.000000  5.000000  6.00000
>>> df.T.describe()  # describeで行ごとに処理したい場合は転置する
              0         1         2         3
count  6.000000  6.000000  6.000000  6.000000
mean   3.500000  3.500000  2.500000  2.500000
std    1.870829  1.870829  1.870829  1.870829
min    1.000000  1.000000  0.000000  0.000000
25%    2.250000  2.250000  1.250000  1.250000
50%    3.500000  3.500000  2.500000  2.500000
75%    4.750000  4.750000  3.750000  3.750000
max    6.000000  6.000000  5.000000  5.000000

 参考:
pandas.DataFrame.describe — pandas 0.24.2 documentation
pandasのdescribeで各列の要約統計量(平均、標準偏差など)を取得 | note.nkmk.me

【python】辞書で同じキーに複数の値を登録する

 ちょっとしたTips。

 辞書(dict)は通常、一つのキーには一つの値しか登録できない。代入しても上書きされる。

>>> d = {}
>>> d["hoge"] = 1
>>> d
{'hoge': 1}
>>> d["hoge"] = 2
>>> d
{'hoge': 2}

 こういうときどうすれば良いのかというと、値をリスト等にしておいて、そのリストにappendしていけば良い。どのように使いたいかにもよるのだが、大抵はこれで用が済む。

スポンサーリンク


 defaultdictを使うと面倒な初期化処理が省けて便利。

>>> from collections import defaultdict
>>> d = defaultdict(list)
>>> d["hoge"].append(1)
>>> d["hoge"].append(2)
>>> d
defaultdict(<class 'list'>, {'hoge': [1, 2]})
>>> d["hoge"]
[1, 2]

 別にlistじゃなくても、setだろうがdictだろうが何でも指定できる。tupleも指定はできるが、変更できないと何の役にも立たない。

>>> d = defaultdict(set) # setの場合
>>> d["hoge"].add(1) # setの場合はaddを使う
>>> d["hoge"].add(1)
>>> d["hoge"].add(2)
>>> d["hoge"]
{1, 2} # 重複しないことに注目
>>> d = defaultdict(dict) # dictの場合
>>> d["hoge"][1] = 1
>>> d["hoge"][2] = 2
>>> d["fuga"][-1] = -1
>>> d["fuga"][-2] = -2
>>> d
defaultdict(<class 'dict'>, {'hoge': {1: 1, 2: 2}, 'fuga': {-2: -2, -1: -1}}) # よくわからないけど何でもできそう

 参考:
8.3. collections — コンテナデータ型 — Python 3.6.5 ドキュメント

 mutableなコレクション型を辞書の値にしておく、という発想があればそれほど難しい話ではない。

【python】# coding: utf-8はもうやめる

 pythonのプログラムは先頭行(あるいはシェバンの次の二行目)でファイルの文字コードを指定することができます。エンコーディング宣言といいます。

 こんなのとか

# coding: UTF-8

 こういうのもありますね。これはemacsに自動認識させるための書式らしい*1

# -*- coding: utf-8 -*-

 これをずっと書いてたんだけど、PEP8を読んでいたらこんな記述に気づきました。

ASCII (Python 2) や UTF-8 (Python 3) を使用しているファイルにはエンコーディング宣言を入れるべきではありません。

はじめに — pep8-ja 1.0 ドキュメント

 えぇぇぇぇ!? と思ったんだけど、何回読み直しても「デフォルトエンコーディング使うならエンコーディング宣言は書くなよ!」と書いてあるようにしか読めない。

 デフォルトエンコーディングを使うなら不要なのは知っていたけど、コーディング規約で非推奨にされてたのですね・・・。

 ということで、「PEP8に準拠しろよ!」というのはpython使いの常識なので(本当か?)、そして私はpython3しか書かないので、個人的には今後コーディング宣言は使わないことにしました。シェバンも必要に迫られない限りは書かない人間なので、今後私のプログラムは一行目からimportで始まることに。・・・ちょっと寂しい気もするけど、すっきりはする。

 「入れるべきではない」とまで言い切っているのは少し不思議な感じはするけど、不要なものをわざわざ書けというよりは良いのかもしれませんね。

*1:そのくせemacs使いの僕は面倒くさくて上ので済ませてきたんだけど

map・filterとリスト内包表記はどちらを使うべきか?

はじめに

 pythonにはmap・filterという関数と、リスト内包表記という独自の記法があります。

 どちらを使っても同じようなことができますが、どちらを使うべきなのでしょうか?

 色々な視点から考えてみます。

 目次


スポンサーリンク


返り値の型

 map・filterにはmapオブジェクト、filterオブジェクトというジェネレータのようなものを返すという厄介な性質があります。ただしこれはpython3以降の仕様なので、python2を使っている方には当てはまりません(python3とpython2のソースコード互換性に効いてくるので、それはそれで難しい問題ではある)。

 リスト内包表記でジェネレータの返り値を望むのならジェネレータ式を使うことができますが、map・filterでリストを返したければlist()を使ってリストに変換するしかありません。

>>> [x*2 for x in range(10)]
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
>>> list(map(lambda x:x*2, range(10)))
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]

 list()は厄介で、何しろ6文字も余計なものが増えてしまいます。ソースコードが冗長になり、list()が増えすぎると可読性も損ないます。

 なお、star operatorを使って3文字で済ませる裏技もあります。これができるのは比較的新しいバージョンのpythonだけのはずですが・・・。

>>> [*map(lambda x:x*2, range(10))]
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]

 可読性は悪いです(というか初見殺し。そしてこの書き方は流行っていない)。

冗長さ

 上述の通り、map・filterはリストに変換してやる必要があるので、基本的に冗長です。lambdaが必要になるとなおさら、という気がします。再掲しますが、

>>> [x*2 for x in range(10)]
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
>>> list(map(lambda x:x*2, range(10)))
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]

 リスト内包表記の方は24文字、mapの方は34文字(どちらもスペース込み)ですから、勝ち目はなさそうです。

 しかし、「返り値はジェネレータのようなもので構わない(もしくは積極的にジェネレータを使いたい)」かつ「既存の関数を使い、lambdaを使わない」という条件だと、mapも輝き始めるようです。数字の文字列を一文字ずつintのリストに変換する例。

>>> (int(c) for c in "1234")
<generator object <genexpr> at 0x7fefabac23b8>
>>> map(int, "1234")
<map object at 0x7fefaba4b0f0>

 上は24文字、下は16文字です。内包表記の「for * in」がなくなるのと、呼び出しを省けるので短くなります。

 filterも事情は同様です。よって、冗長さについては「ケースバイケース」であり、より詳しく言うと「返り値がジェネレータで構わず、使う関数がすでに定義されている」場合はmap・filterが有利と言い切れます。逆に言うと、それ以外のケースでは内包表記の方が簡潔に書けるようです。

 たとえば、入力された数値n個をint型に変換する、といったコードでは、

a,b,c = map(int, input())

 のように、極めて簡潔な記述がmapを使うことで可能になります。

可読性

 可読性ははっきり言って、一長一短です。

 一般的なプログラミング言語の構文をベースに考えると、わかりやすいのはmap・filterの方で、何しろ普通の関数です。初心者でも安心して使えると言えます(無名関数と高階関数の概念さえ理解すれば)。ただし、現実的には増えまくるlist()のせいでコードがごちゃごちゃし、かえって可読性が下がります。また、mapとfilterはそれぞれ独立に使う必要があるので、なおさらコードがごちゃごちゃします。

 リスト内包表記は一見するとわかりづらいですが、慣れてしまえば簡潔で、書きやすく、読みやすいと言えます(本当か? 個人差は確実にあると思う)。また、一つの内包表記でfilterとmapを同時に適用できるので簡潔になります。

 多重の複雑なものになると、そもそも簡潔に書け、改行とインデントでかなりわかりやすく整理できる内包表記の方が相対的にかなり有利と感じています。

from pprint import pprint
pprint([[x*y for y in range(1,10)]
        for x in range(1,10)])

pprint(list(map(
    lambda x:list(map(
        lambda y:x*y,
        range(1,10))), 
    range(1,10))))

 かけ算九九の表を表示するプログラムですが、mapの方ははっきり言ってひでーと思います。この程度ならまだ読めなくはありませんが、条件が加わってfilterを追加するとかやり始めるともはや読めなくなります。

事故

 これはmap・filterの欠点というよりはジェネレータの欠点ですが、割と事故のもとになってくれます。

>>> result = map(lambda x:x*2, range(10))
>>> result[2]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'map' object is not subscriptable
>>> for x in result:
...     print(x)
... 
0
2
4
6
8
10
12
14
16
18
>>> for x in result:
...     print(x)
... 
>>> # 何も表示されない

 「んなもん、ちゃんと理解して使えば良いだろ」という意見も当然あると思います。それはそれで正しいと思いますが、理想論です。私は自分自身を含めて人間を100%信頼はしていません。事故る要素があれば、必ず事故るのです。

 そして、これが嫌という理由で私は普段、ほとんどリスト内包表記を使っています。*1

まとめ

 基本的には内包表記がずいぶん有利です。特別な理由がなければ内包表記で書けば良いと思います。

 map・filterはリストに既存の関数を適用するジェネレータを作るときに威力を発揮する可能性があります。あるいは、代入でunpackするときとか。

 というか、それ以外なさそうです。無名関数と組み合わせてまで使う必然性はないと思います。

*1:余談ですが、私はpythonのジェネレータというものはあまり好きではありません。mapやfilterやzip, rangeなどの返り値は、基本的にすべてlistでも別に良いと思います。これらの関数にはgeneratorオプションを追加してdefault=Falseとするか、python2よろしく対応するxrangeなどを作る、あるいはジェネレータを示す糖衣構文を新しく作って、その糖衣構文で囲むと値の評価時にジェネレータとして処理されるような仕組みを導入すれば良いと思います。これらについては、python4に期待です。

【python】mapで複数の引数を渡したいときはstarmapが便利

 pythonにはmapという関数があります。関数型プログラミングに役立ち、大変便利な関数です。

 しかし、これはデフォルトでは一つの引数を前提としています。

>>> list(map(lambda a,b: a+b, zip([1,2,3],[4,5,6])))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: <lambda>() missing 1 required positional argument: 'b'

 困ったものですね。こういうときは関数の側で工夫してあげるしかないのでしょうか?

>>> list(map(lambda x: x[0]+x[1], zip([1,2,3],[4,5,6])))
[5, 7, 9]

 内包表記使えよ、という考え方もあると思いますが、不問とします。

 こういうケースでうまく働くのが、itertools.starmapです。
10.1. itertools — 効率的なループ実行のためのイテレータ生成関数 — Python 3.6.5 ドキュメント

 starはunpackingのstar。こう書くともうだいたい察してしまった人もいるでしょうが、こう書けます。

>>> from itertools import starmap
>>> list(starmap(lambda a,b: a+b, zip([1,2,3],[4,5,6])))
[5, 7, 9]

 引数が格納されたtupleは、lambdaに渡されるときにunpackされる訳です。気が利いてますね。

 これと同様にデフォルト引数を渡せるitertools.starstarmapもあります・・・と言いたいところですが、ありませんでした。残念。

それでも内包表記で書きたい人へ

 まあ、これのためにimportするくらいなら・・・という気持ちはとてもよくわかりますし、私も普段はそうするクチなんですが。

>>> [a+b for a,b in zip([1,2,3],[4,5,6])]
[5, 7, 9]

 listとlambdaが端折れる分だけ短く書けているように見えますが、返り値がジェネレータで構わなくて、既存の関数を呼び出したいときはmapの方が短くなる(こともある)ようです。

追記

 この記事書いたときはすっかり忘れていたけど、そういえばmapはデフォルトで複数の引数を取れる。

>>> list(map(lambda a,b:a+b, [1,2,3],[4,5,6]))
[5, 7, 9]

 普通にこっちで良いですね。

 zipされたものを入力に渡したいときは工夫が必要だが、unpackしてもう一回zipしてまたunpackすれば良い。

>>> list(map(lambda a,b:a+b, *zip(*zip([1,2,3],[4,5,6]))))
[5, 7, 9]

【python】numpyで二次元配列を結合して三次元配列にする方法

 複数の二次元配列を結合して三次元配列に変換する方法について。

np.dstack

 そのものずばりのnp.dstackという関数がある。

numpy.dstack — NumPy v1.14 Manual

>>> a = np.array([[1,2,3],[4,5,6]])
>>> b = np.array([[7,8,9],[10,11,12]])
>>> a
array([[1, 2, 3],
       [4, 5, 6]])
>>> b
array([[ 7,  8,  9],
       [10, 11, 12]])
>>> np.dstack([a,b])
array([[[ 1,  7],
        [ 2,  8],
        [ 3,  9]],

       [[ 4, 10],
        [ 5, 11],
        [ 6, 12]]])

 なんか思ってたのと違う。理解するために色々やってみる。

>>> c = np.dstack([a,b])
>>> a.shape
(2, 3)
>>> c.shape
(2, 3, 2)
>>> c[:,0]
array([[ 1,  7],
       [ 4, 10]])
>>> c[:,1]
array([[ 2,  8],
       [ 5, 11]])
>>> c[:,2]
array([[ 3,  9],
       [ 6, 12]])
>>> c[:,0:,0] # 元の配列を取り出す。適当に書いてたらできた
array([[1, 2, 3],
       [4, 5, 6]])
>>> c[:,0:,1]
array([[ 7,  8,  9],
       [10, 11, 12]])

 わかるようなわからないような感じ。

 ドキュメントによると、

This is equivalent to concatenation along the third axis after 2-D arrays of shape (M,N) have been reshaped to (M,N,1) and 1-D arrays of shape (N,) have been reshaped to (1,N,1).

 つまりこれと同じということ。

>>> a_reshaped = a.reshape((2,3,1))
>>> b_reshaped = b.reshape((2,3,1))
>>> a_reshaped
array([[[1],
        [2],
        [3]],

       [[4],
        [5],
        [6]]])
>>> b_reshaped
array([[[ 7],
        [ 8],
        [ 9]],

       [[10],
        [11],
        [12]]])
>>> np.concatenate((a_reshaped, b_reshaped), axis=2)
array([[[ 1,  7],
        [ 2,  8],
        [ 3,  9]],

       [[ 4, 10],
        [ 5, 11],
        [ 6, 12]]])

 理解はできたが、ほしいものとは違うのだった。

単純なやり方

 普通はこっちがほしいと思うの。

>>> np.array([a,b])
array([[[ 1,  2,  3],
        [ 4,  5,  6]],

       [[ 7,  8,  9],
        [10, 11, 12]]])

 わかりやすい。上との対比で更にわかりやすくするため、np.concatenateで書いてみる。

>>> np.concatenate([[a], [b]], axis=0)
array([[[ 1,  2,  3],
        [ 4,  5,  6]],

       [[ 7,  8,  9],
        [10, 11, 12]]])

 つまりこっちは(M, N)の配列を(1, M, N)にしてaxis=0で結合したのと等価である。

追記

 こんな記事も書きました。numpyの配列の結合方法についてまとめています。併せて御覧ください。

www.haya-programming.com

【python】ランダムフォレストのチューニングにOOB誤り率を使う

 一般的な機械学習のアルゴリズムでは、パラメタチューニングにはグリッドサーチ・交差検証を組み合わせて使うのが割と普通だと思います。sklearnにはそれ専用のGridSearchCVというクラスまで用意されています。

 実際問題としては、GridSearchは良いとしても交差検証をやったのでは計算コストをたくさん食います。なので、他の方法で代替できるなら、そうしたいところです。

 そして、RandomForestはOOB誤り率という、いかにも強そうなスコアを計算できます。これの良いところは一回fitしただけで計算でき、交差検証と同じ(ような)結果が得られることです。

 なので、OOB誤り率を使ったパラメタチューニングを試してみたいと思います。

OOB誤り率とはなんぞ?

 OOB誤り率、OOBエラー、OOB error等と呼ばれます。これを理解するためにはランダムフォレストの学習過程を(少なくとも大雑把には)理解する必要があります。

 ランダムフォレストは木を一本作る度に、データをランダムサンプリングします(「ランダム」の所以です)。サンプリング方法はブートストラップサンプリングです。詳細はググって頂くとして、とにかくサンプリングの結果、各決定木に対して「訓練に使われなかったデータ」が存在することになります。逆に言えば、各データに対して「そのデータを使っていない決定木の集合」があります。このことを利用して、「そのデータを使っていない決定木」だけ利用して推定し、汎化性能を見ようというのがOOB誤り率のコンセプトです。

実験

 以下のようなコードを書きました。

# coding: UTF-8
import time
from itertools import product

from sklearn.datasets import load_digits
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.model_selection import  train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import precision_recall_fscore_support as prf

def oob_tune(clf, params, X, y):
    key_index = list(params.keys())
    
    result_dict = dict()
    for ps in product(*[params[k] for k in key_index]):
        clf.set_params(**dict(zip(key_index, ps)))
        clf.fit(X, y)
        result_dict[ps] = clf.oob_score_

    return dict(zip(key_index, 
                    sorted(result_dict.items(), 
                           key=lambda x:x[1],
                           reverse=True)[0][0]))
    
def main():
    digits = load_digits()
    X_train, X_test, y_train, y_test = train_test_split(
        digits.data, digits.target)
    
    params = {"n_estimators":[50, 100, 500],
              "max_features":[4, 8, 12],
              "min_samples_leaf":[1, 2]}

    rfc = RFC(oob_score=False, n_jobs=-1)

    t1 = time.time()
    clf = GridSearchCV(rfc, params, cv=6, n_jobs=-1)
    clf.fit(X_train, y_train)
    t2 = time.time()
    best_params = clf.best_params_
    print("GridSearchCV") 
    print("time:{0:.1f}".format(t2-t1))
    print("best params")
    print(best_params)
    rfc = RFC(**best_params)
    rfc.fit(X_train, y_train)
    y_pred = rfc.predict(X_test)
    score = prf(y_test, y_pred, average="macro")
    print("p:{0:.3f} r:{0:.3f} f:{0:.3f}".format(*score))


    rfc = RFC(oob_score=True, n_jobs=-1)

    t1 = time.time()
    best_params = oob_tune(rfc, params, X_train, y_train)
    t2 = time.time()

    print("oob")
    print("time:{0:.1f}".format(t2-t1))
    print("best params")
    print(best_params)
    rfc = RFC(**best_params)
    rfc.fit(X_train, y_train)
    y_pred = rfc.predict(X_test)
    score = prf(y_test, y_pred, average="macro")
    print("p:{0:.3f} r:{0:.3f} f:{0:.3f}".format(*score))

if __name__ == "__main__":
    main()

 sklearnの機能でありそうな気がしましたが、見つからなかったので*1、OOB誤り率を使ったパラメタチューニングのプログラムは自分で書きました。

 同じパラメタ候補に対してのグリッドサーチで、GridSearchCVと比較しています。

 細かい説明は抜きで(気になる人はプログラムを追ってください。簡単なので)、結果を貼ります。

GridSearchCV
time:41.8
best params
{'n_estimators': 500, 'max_features': 4, 'min_samples_leaf': 1}
p:0.977 r:0.977 f:0.977
oob
time:13.1
best params
{'n_estimators': 500, 'max_features': 8, 'min_samples_leaf': 1}
p:0.979 r:0.979 f:0.979

 やはりOOB誤り率を使った方が速いです。最適パラメータは困ったことに両者で食い違っています。テストデータでのスコアはOOB誤り率で計算した最適パラメータを用いた方がびみょ~に高いですが、これくらいだとほとんど誤差かもしれません(1データか2データ程度の違い)。

まとめ

 交差検証でやるのと同程度の結果が得られる・・・と言い切るには微妙な部分がありますが、とにかくやればできます。

 自分でパラメタチューニングのプログラムを書くと、sklearnのインターフェースから逸脱するデメリットがあるので、実際にやるべきかどうかは微妙。scoringの関数と実質的にcross validateしないcv(ダミー)を渡せば、GridSearchCVを利用してもやれそうではありますが。ちょっと悩ましいところです*2

 とにかくさっさとパラメタチューニングしたいんだ! というときは使えるでしょう。

*1:本当に存在しないか、存在するけど私が見つけられなかったかのどちらか。

*2:それとも、私が見つけられてないだけで、sklearnで直接これができる方法があるのだろうか? もしそうなら、ご存知の方は教えていただきたい

【python】bitのリストを高速にintに変換する

やりたいこと

input:[0,1,0,0]
output:4

 極めて単純明快ですが、やるだけなら簡単なので速度を測ります。

 さらに、pure pythonでやると遅いことが目に見えているのでcythonで高速にしようというネタです。

pure pythonで書いたプログラム

 素晴らしいことに(?)、bitリストのintへの変換は既出ネタです。

arrays - Bits list to integer in Python - Stack Overflow

 書き方はだいたい数パターンなので、ここの回答の一つから引用します。

def mult_and_add(bit_list):
    output = 0
    for bit in bit_list:
        output = output * 2 + bit
    return output

def shifting(bitlist):
     out = 0
     for bit in bitlist:
         out = (out << 1) | bit
     return out

 リンク先は後者の方がじゃっかん遅いと主張しています。本当でしょうか?

 自分でも動かしてみます。

# coding: UTF-8

import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def mult_and_add(bit_list):
    output = 0
    for bit in bit_list:
        output = output * 2 + bit
    return output

def shifting(bitlist):
     out = 0
     for bit in bitlist:
         out = (out << 1) | bit
     return out

def bench(lst):
    t1 = time.time()
    mult_and_add(lst)
    t2 = time.time()
    res1 = t2-t1

    t1 = time.time()
    shifting(lst)
    t2 = time.time()
    res2 = t2-t1

    return np.array((res1, res2))

def main():
    result_lst = []
    for i in range(50):
        i = 10**i
        binlst = [int(x) for x in bin(i)[2:]]
        lst = []
        for j in range(100):
            lst.append(bench(binlst))
        result_lst.append(np.mean(lst, axis=0))
    df = pd.DataFrame(result_lst, columns=["muladd", "shift"])
    print(df)
    df.plot()
    plt.savefig("fig1.png")

if __name__ == "__main__":
    main()

ベンチマークの結果
ベンチマークの結果

 横軸は 10^nのnです。本当らしい。なんでだろう。

cythonを使う

 どうせpure pythonは遅いと思ったんです。だからcythonで書きました。

  • cython_bit2int.pyx
# coding: UTF-8
import array

def muladd_cy(lst):
    return _muladd(array.array("Q", lst))

def shift_cy(lst):
    return _shift(array.array("Q", lst))

cdef int _muladd(unsigned long[:] bitlist):
    cdef unsigned long int out = 0
    cdef int i = 0

    for i in range(len(bitlist)):
        out = out * 2 + bitlist[i]
    return out

cdef int _shift(unsigned long[:] bitlist):
    cdef unsigned long int out = 0
    cdef int i = 0

    for i in range(len(bitlist)):
        out = (out << 1) | bitlist[i]
    return out
  • bit2int.py
# coding: UTF-8

import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import pyximport; pyximport.install()
import cython_bit2int as cyb2i


def mult_and_add(bit_list):
    output = 0
    for bit in bit_list:
        output = output * 2 + bit
    return output

def shifting(bitlist):
     out = 0
     for bit in bitlist:
         out = (out << 1) | bit
     return out

def bench(lst):
    t1 = time.time()
    mult_and_add(lst)
    t2 = time.time()
    res1 = t2-t1

    t1 = time.time()
    shifting(lst)
    t2 = time.time()
    res2 = t2-t1

    t1 = time.time()
    cyb2i.muladd_cy(lst)
    t2 = time.time()
    res3 = t2-t1

    t1 = time.time()
    cyb2i.shift_cy(lst)
    t2 = time.time()
    res4 = t2-t1

    return np.array((res1, res2, res3, res4))

def main():
    result_lst = []
    for i in range(50):
        i = 10**i
        binlst = [int(x) for x in bin(i)[2:]]
        lst = []
        for j in range(100):
            lst.append(bench(binlst))
        result_lst.append(np.mean(lst, axis=0))
    df = pd.DataFrame(result_lst, columns=["muladd", "shift", 
                                           "muladd_cy", "shift_c"])
    print(df)
    df.plot()
    plt.savefig("fig2.png")

if __name__ == "__main__":
    main()

 array使うのはフェアな比較じゃないって? まあ大目に見て下さい。

ベンチマークの結果2
ベンチマークの結果2

 桁数が少ない領域では速度差が少ない(というかヘタしたら負けてる)ので微妙かも。恐らくlistからarrayへの変換が重いのでしょう。こっちだとshiftの方が速くなるのは理屈通りでした。

【python】順列・組み合わせを計算する方法

 Pythonで、順列(Permutation)と組み合わせ(Combination)がほしくなるときがある。また、順列・組み合わせの数がほしくなることもある。

 順列・組み合わせそのものはitertoolsで、その数はscipyで出せる。計算方法についてまとめておく。

 目次

スポンサーリンク


順列・組み合わせそのものがほしい場合

 つまり"abc"から2個取り出す→[["a","b"], ["b", "c"], ["a", "c"]]という結果を期待している場合。

 これは標準ライブラリのitertoolsでできる。

順列の場合

 itertools.permutationsを使う。

>>> from itertools import permutations
>>> list(permutations("abc", 2))
[('a', 'b'), ('a', 'c'), ('b', 'a'), ('b', 'c'), ('c', 'a'), ('c', 'b')]

組み合わせの場合

 itertools.combinationsを使う。

>>> from itertools import combinations
>>> list(combinations("abc", 2))
[('a', 'b'), ('a', 'c'), ('b', 'c')]

 極めて容易。

順列・組み合わせの数だけほしい場合

 順列・組み合わせの数がほしいときもある。つまりnPrとnCrである。数字がほしいだけで、結果そのものは要らないケース。

 困ったことに、直接これを計算してくれる関数はmathモジュールなどには用意されていないようだ。実現する方法は何通りかある。

itertoolsの関数の結果をlen()する

 安直。パフォーマンス上はよくないだろうし、無駄っぽい。が、一応できる。

>>> from itertools import permutations
>>> from itertools import combinations
>>> len(list(permutations("abc", 2)))
6
>>> len(list(combinations("abc",2)))
3

 まあ、他の方法を使った方が良さげ。

自分で実装する

 数学的定義通り書くことができる。

>>> def P(n, r):
...     return math.factorial(n)//math.factorial(n-r)
... 
>>> def C(n, r):
...     return P(n, r)//math.factorial(r)
... 
>>> P(3, 2)
6
>>> C(3, 2)
3

 もうちょっと工夫して、無駄な階乗の計算を端折ることも考えられる。ただ、あまりごちゃごちゃ書くのも大変だし、実はscipyにあるので車輪の再発明をする必要すらない。

 nPrとかnCrを自分で実装している人は、素直にライブラリを使いましょう。

ライブラリを使う(おすすめ)

 scipyにある。どちらもscipy.specialから使える。

scipy.special.perm — SciPy v1.3.0 Reference Guide

scipy.special.comb — SciPy v1.3.0 Reference Guide


 注意点としては、exact=Trueを指定しないとfloatで近似値が返されること。この機能は近似値の方が便利/それで十分という用途では活用すれば良いと思う(ちなみにデフォルトはFalse。なので正確な値がほしいときは一々指定してやる必要がある)。

 あと、使い方はいまいちよくわからないけどndarrayも渡せるし、組み合わせの方はrepetitionオプションで重複ありの組み合わせも計算できる。

>>> from scipy.special import perm, comb
>>> perm(3, 2, exact=True)
6
>>> comb(3, 2, exact=True)
3

 あっさりしたものですね。

まとめ

 順列・組み合わせそのものがほしいときはitertools、順列・組み合わせの数がほしいときはscipyを使おう。
(余談だが、scipyは色々あって便利なのに、ドキュメントの検索のしづらさとgoogleの検索順位の低さで損してる気がする。)

【python】sklearnのtolってなんだ?

 sklearnの公式ドキュメントをよく読む方なら、色々なモジュールに"tol"というオプションがあることに気づいていると思います。

 たとえばSVCだと、こんな風に書いてあります。他のモジュールも似たり寄ったりですが。

tol : float, optional (default=1e-3)

Tolerance for stopping criterion.

 出典:
sklearn.svm.SVC — scikit-learn 0.20.1 documentation

 よくわからないし、なんとなく重要そうじゃないし、デフォルトから変える必要もないでしょ、ということで無視されがちなパラメタですが、冷静に考えたら何なのかまったくわからない。実は重要だったりすると大変ですね。ということで、調べました。

 まずtoleranceという単語とcriterionという単語に対応する日本語がぱっと出てきません(そこからかよ)。仕方がないのでググると、「耐性」や「公差」「許容誤差」のような意味が出てきます。なんとなく雰囲気は伝わりました。

toleranceの意味・使い方 - 英和辞典 Weblio辞書

 ちなみに後者のcriterionは基準という意味でした。これは雰囲気でもなんでもなく、そのまんま基準です。

 よって、「Tolerance for stopping criterion.」を訳すと「打ち切るための許容誤差の基準」となり、なんとなくわかるようでわかりません。

 仕方ないので色々キーワードを変えて検索していると、CrossValidatedの質問を見つけました。

machine learning - What exactly is tol (tolerance) used as stopping criteria in sklearn models? - Cross Validated

tol will change depending on the objective function being minimized and the algorithm they use to find the minimum, and thus will depend on the model you are fitting. There is no universal tolerance to scikit .
超訳:目的関数とそれを最小化するアルゴリズムに依存するから、モデルによってちげーよ。sklearnには普遍的なtolなどない

 あっ、そうですか・・・。

 これで終わってしまうのも寂しいので、SVMでtolを変えて結果が変わるか実験してみます。

# coding: UTF-8

import time

from sklearn.datasets import load_digits
from sklearn.svm import SVC
from sklearn.model_selection import cross_validate, StratifiedKFold as SKF

def main():
    digits = load_digits()
    
    svm1 = SVC(C=5, gamma=0.001, tol=1)
    svm2 = SVC(C=5, gamma=0.001, tol=0.5)
    svm3 = SVC(C=5, gamma=0.001, tol=0.001)
    svm4 = SVC(C=5, gamma=0.001, tol=0.00001)

    skf = SKF(random_state=0)

    scoring = {"p": "precision_macro",
               "r": "recall_macro",
               "f":"f1_macro"}

    skf = SKF(n_splits=5, shuffle=True, random_state=0)

    for svm, tol in zip([svm1, svm2, svm3, svm4], [1, 0.5, 0.001, 0.00001]):
        t1 = time.time()
        scores = cross_validate(svm, digits.data, digits.target,
                                cv=skf, scoring=scoring)
        t2 = time.time()
        print("tol:{0:5} time:{1:8.3f} p:{2:.3f} r:{3:.3f} f:{4:.3f}".format(
            tol,t2-t1,
            scores["test_p"].mean(),
            scores["test_r"].mean(),
            scores["test_f"].mean()))

if __name__ == "__main__":
    main()

 結果は、

tol:    1 time:   0.661 p:0.990 r:0.989 f:0.989
tol:  0.5 time:   0.833 p:0.991 r:0.990 f:0.991
tol:0.001 time:   1.193 p:0.992 r:0.992 f:0.992
tol:1e-05 time:   1.226 p:0.992 r:0.992 f:0.992

 性能への影響はごく僅かですが、トータルの処理時間は若干(数10%程度)削れるようです。少しでも軽くしたい、というときは検討する(性能への悪影響が抑えられる範囲でtolを上げる)価値はありますね。