静かなる名辞

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


【python】その矛盾した__eq__は・・・

私は疑問を持った

 pythonでは比較演算子==を使うと、内部的には__eq__メソッドが呼ばれる。

 ここから、素朴な疑問が生じる。比較演算子は二項演算子なので、2つのオブジェクトに対して適用される。

 どちらのオブジェクトの__eq__が呼ばれるのだろう? また、もし2つのオブジェクトの__eq__が違う値を返すとしたら、一体どんな事態が生じるのだろう? まずいことにならないのだろうか?

 ドキュメントにはちゃんとこのことが書いてある。

x==y は x.__eq__(y) を呼び出します

 3. データモデル — Python 3.6.5 ドキュメント

 これだけ。いまいち納得がいかない。

検証した

 こんなコードを書いてみた。

class Hoge:
    def __init__(self, name):
        self.name = name
    def __eq__(self, other):
        print("__eq__ of Hoge!  by", self.name)
        return True

class Fuga:
    def __init__(self, name):
        self.name = name
    def __eq__(self, other):
        print("__eq__ of Fuga!  by", self.name)
        return False

h1 = Hoge("h1")
h2 = Hoge("h2")
f1 = Fuga("f1")
print("h1 == h1")
print(h1 == h1)
print("h1 == h2")
print(h1 == h2)
print("h2 == h1")
print(h2 == h1)
print("h1 == f1")
print(h1 == f1)
print("f1 == h1")
print(f1 == h1)

 常にTrueを返す__eq__を持つHogeクラス、常にFalseを返す__eq__を持つFugaクラスを定義し、片っ端から比較している。どの__eq__がprintされたかもこれでわかるはずだ。

 結果は、こんなものだった。

h1 == h1
__eq__ of Hoge!  by h1
True
h1 == h2
__eq__ of Hoge!  by h1
True
h2 == h1
__eq__ of Hoge!  by h2
True
h1 == f1
__eq__ of Hoge!  by h1
True
f1 == h1
__eq__ of Fuga!  by f1
False

 とりあえず、興味深いのは二回printされたりはまったくしなかったこと。そして、常に左側のオブジェクトの__eq__メソッドが呼ばれているように見えること。

 要するにドキュメントの「x==y は x.__eq__(y) を呼び出します」は極めて正しく、それ以外の動作はまったくないということだ。

 意外な感じがする・・・。

まずい事態

 ということは、__eq__を好き勝手に拡張しているpython外部ライブラリとかだと、まずい事態が生じるのではないだろうか。

 たとえばnumpyのarrayをintと比較すると、次のように動作する。

>>> import numpy as np
>>> a = np.array([0,0,0,1,1,1,2,2,2])
>>> a == 0
array([ True,  True,  True, False, False, False, False, False, False])

 慣れない頃は不思議な感じがしていたが、これはnumpy.ndarrayクラスの__eq__がbooleanのndarrayを返すように書かれているからだ、と考えるとそれほどおかしくはない。

 ここまでは良い。だけど、このコードで0==aとするとintクラスの__eq__が呼ばれて違う結果が出てしまうのでは?

>>> 0 == a
array([ True,  True,  True, False, False, False, False, False, False])

 理解不能すぎる・・・。

 いやまてよ、pythonのintがndarrayに忖度しているだけなのでは? と思い、次のコードを実行してみた。

>>> class Hoge:
...     def __init__(self, name):
...         self.name = name
...     def __eq__(self, other):
...         print("__eq__ of Hoge!  by", self.name)
...         return True
... 
>>> h = Hoge("h")
>>> a == h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
__eq__ of Hoge!  by h
array([ True,  True,  True,  True,  True,  True,  True,  True,  True])
>>> h == a
__eq__ of Hoge!  by h
True

 とりあえず、a == hで__eq__ of Hoge! by hがいっぱい出てきたことは理解可能。numpy.ndarray.__eq__は各要素に対してotherの__eq__を呼ぶ感じに動作しているのだろう。

 h == aは普通にTrueになった。ということは、intとは違う動作になっている。intの方は、一体どんな挙動なのだろう?

 ググったら出てきた。

It works because int.__eq__() returns NotImplemented and when that happens it results in a call to other.__eq__(self) and that's what is returning True and False here.

python - Why is `int.__eq__(other)` a working comparison? - Stack Overflow

 NotImplementedはpythonの組み込み定数らしい。

3. 組み込み定数 — Python 3.6.7 ドキュメント

 intの__eq__は(int型、およびその他のintと直接比較できることになっている型以外との比較時には)NotImplementedを返す。これが返ると、otherの__eq__が呼ばれ、判定が行われる、とか。

 なるほど、実験してみよう!

>>> obj = 0
>>> obj.__eq__(h)
NotImplemented
>>> obj == h
__eq__ of Hoge!  by h
True

 すげえ、よくできてる・・・。つまり、intの忖度でもなんでもなく、この機能(NotImplementedが返ったときの挙動)をうまく使ってnumpy.ndarrayは実装されていたのだった*1

 逆に、自作オブジェクトでやるときは、想定していない型が来たときはNotImplementedを返すみたいな注意が必要という話でもあるが。

まとめ

 ま、私が心配する程度のことは、聡明なpython開発陣の皆さんはしっかり対処しておられるのでした。つーか、これを考えた人たちは本当に頭が良い。とてもよく出来たシステムだと思う。

 逆に、一瞬で終わりそうな比較なのに内部ではメソッドがバンバン呼ばれていて、そのコストは実はけっこう高いという現実も見てしまった感も、あるといえばある。動的型付け言語だから、仕方ないのだけど。

 とにかく、これで安心して__eq__を使えるようになったね!*2

*1:とすると、上の「__eq__ of Hoge! by h 」がいっぱい出てきたのも、おそらく一回np.int64の__eq__が呼ばれて、NotImplementedが返ってからHoge.__eq__が呼ばれる、という流れのような気がする。確証はない

*2:現実に実装する際どんな風にすれば良いのかについては、日本語圏だとこちらの記事が参考になりました:Pythonにおける同値性比較の実装

【python】MeanShiftのbandwidthを変えるとどうなるか実験してみた

 前回の記事ではMeanShiftクラスタリングを試してみました。

www.haya-programming.com

 このMeanShiftにはbandwidthというパラメータがあり、クラスタ数を決定する上で重要な役割を果たしているはずです。

 いまいち結果に納得がいかないというとき、bandwidthをいじって改善が見込めるのかどうか確認してみます。

プログラム

 例によってirisとwineで比較。簡単に書きました。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from sklearn.datasets import load_iris, load_wine
from sklearn.cluster import MeanShift, estimate_bandwidth
from sklearn.decomposition import PCA

def process(dataset, name):
    origin_bandwidth = estimate_bandwidth(dataset.data)
    rates = np.logspace(np.log10(0.2), np.log10(5), 11)
    fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(24,18))

    PCA_X = PCA().fit_transform(dataset.data)
    for target in range(3):
        axes[0,0].scatter(PCA_X[dataset.target==target, 0],
                        PCA_X[dataset.target==target, 1],
                        c=cm.Paired(target/3))
    axes[0,0].set_title("original label", fontsize=28)

    for r, ax in zip(rates, axes.ravel()[1:]):
        ms = MeanShift(bandwidth=r*origin_bandwidth, n_jobs=-1)
        y = ms.fit_predict(dataset.data)
        n_cluster = ms.cluster_centers_.shape[0]
        for target in range(n_cluster):
            ax.scatter(PCA_X[y==target, 0],
                       PCA_X[y==target, 1],
                       c=cm.Paired(target/n_cluster))
        ax.set_title("r:{0:.3f} b:{1:.3f}".format(
            r, origin_bandwidth), fontsize=28)
    fig.savefig(name+".png")

def main():
    iris = load_iris()
    wine = load_wine()

    process(iris, "iris")
    process(wine, "wine")

if __name__ == "__main__":
    main()

 bandwidthをsklearn.cluster.estimate_bandwidthの推定値(デフォルトで用いられる値)の1/5倍から5倍まで変化させ、結果をプロットします。

結果

 プロットされた結果を示します。

 結果の図の見方は、まずタイトルが

  • b

 sklearn.cluster.estimate_bandwidthによる推定値

  • r

 かけた比率

 という風に対応しており、あとは便宜的に2次元上に主成分分析で写像した散布図が、クラスタごとに色分けされて出ています。一枚目が本来のクラスに基づく色分け、r=1の図が推定値による色分けです。

 まずiris。

iris.png
iris.png
 きれいに元通りになるrは今回見た中にはありませんでした。クラスタ数的にはr=0.525とr=0.725の間くらいで3クラスタになりそうですが、この図を見るとそれでうまく元通りまとまるかは疑問です。

 次にwine。

wine.png
wine.png
 こちらもうまく元通りにはならないようです。そもそもデータが悪いという話はあると思います。

結論

 確かにクラスタ数は変わるが、クラスタリングの良し悪しが改善するかはなんともいえないですね。

 データをスケーリングしたり、もっと色々頑張ると改善は見込めるかもしれません。

【python】numpyでの等比数列の作り方

 ※2018年9月2日追記

 以前にこんな記事を書きましたが、改めて確認したところnumpyの機能で普通に書けました。

 その方法は追記に記したので、本文は読み飛ばしてください。

 内容そのものは記録として残しておきます。


 等比数列がほしくなった。作り方をメモしておく。

 目次

スポンサーリンク



比が決まっている場合

 たとえば、3倍ずつ増えていく等比数列がほしいとしよう。

 これは次のように簡単に作れる。

>>> import numpy as np
>>> rates = np.ones(10) * 3
>>> rates
array([3., 3., 3., 3., 3., 3., 3., 3., 3., 3.])
>>> rates ** np.arange(10)
array([1.0000e+00, 3.0000e+00, 9.0000e+00, 2.7000e+01, 8.1000e+01,
       2.4300e+02, 7.2900e+02, 2.1870e+03, 6.5610e+03, 1.9683e+04])

 できたけど、これは面白くもなんともないのだった。

範囲とnが決まっている場合

 たとえば「1から10までの範囲で10個ほしい!」というパターン。これは一筋縄ではいかなさそうである。

 幸い、np.logspaceという便利そうなものがある。

 numpy.logspace — NumPy v1.16 Manual

 これを使えばたぶんできるだろう。

>>> import numpy as np
>>> np.logspace(1,10,10)
array([1.e+01, 1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06, 1.e+07, 1.e+08,
       1.e+09, 1.e+10])

 期待通り動かないのだった。baseというオプションがあり、default=10で、その1,2,3,...,9,10乗が得られている訳である。

 ちょっと意図と違うので、こうしてみた。

>>> np.logspace(0, np.log10(10), 10)
array([ 1.        ,  1.29154967,  1.66810054,  2.15443469,  2.7825594 ,
        3.59381366,  4.64158883,  5.9948425 ,  7.74263683, 10.        ])

 成功した。冷静に考えると、この場合はこれで良い(log10(10)は1だ・・・)。

>>> np.logspace(0, 1, 10)
array([ 1.        ,  1.29154967,  1.66810054,  2.15443469,  2.7825594 ,
        3.59381366,  4.64158883,  5.9948425 ,  7.74263683, 10.        ])

 でもまあ、log10を使う方法だと他の範囲にも対応できる。

>>> np.logspace(0, np.log10(20), 10)
array([ 1.        ,  1.39495079,  1.94588772,  2.71441762,  3.78647901,
        5.2819519 ,  7.368063  , 10.27808533, 14.33742329, 20.        ])

 たとえば23から54まで10個取りたい場合は、どうしたら良いのだろう? こうする。

>>> np.logspace(np.log10(23), np.log10(54), 10)
array([23.        , 25.28791009, 27.80340855, 30.56913458, 33.60997942,
       36.95331033, 40.62921692, 44.67078193, 49.114379  , 54.        ])

 これで一般化できたので、いつでも使えるようになったと思う。

追記

 何もしなくても等比数列を生成するnumpy.geomspaceが用意されていて、普通にこちらでできました。

numpy.geomspace — NumPy v1.16 Manual

>>> import numpy as np
>>> np.geomspace(1, 10, 10)
array([ 1.        ,  1.29154967,  1.66810054,  2.15443469,  2.7825594 ,
        3.59381366,  4.64158883,  5.9948425 ,  7.74263683, 10.        ])
>>> np.geomspace(23, 54, 10)
array([23.        , 25.28791009, 27.80340855, 30.56913458, 33.60997942,
       36.95331033, 40.62921692, 44.67078193, 49.114379  , 54.        ])

 geometric progressionという英単語に馴染みがなくて、記事を書いたときには見つけられなかったのですが……。

 精進しようと思います。

【python】sklearnのMeanShiftクラスタリングを試してみる

はじめに

 MeanShiftはクラスタリングアルゴリズム。クラスタ数を自動で決定してくれるという長所がある。

 理論的には最急降下法で各クラスタの極大点を探していく感じらしいです。わかりやすい解説があったので、リンクを張っておきます(ただし私自身はすべては読み込めていない)。

Mean Shift Clustering

 このMeanShiftはsklearnに実装されているので、簡単に試してみることができます。

 sklearn.cluster.MeanShift — scikit-learn 0.20.1 documentation

 sklearnのトイデータで遊んでみましょう。

 目次

使い方

 いつものsklearnのモデルです。fitしてpredictするだけ。

 いつだったかFuzzy C-Meansをやったときは苦労しましたが、とりあえずそんな心配は要りません。

 となると気になるのはパラメータですが、指定しなくても

  • bandwidth

 勝手に推定される

  • seeds

 乱数のシードなので勝手に決まる。指定するときは[n_samples, n_features]が必要。

  • bin_seeding

 よくわからないけど、初期値の選び方みたいな。上のと関係がありそう。Trueにすると速くなるらしい。デフォルトのFalseの方がアルゴリズムとしては厳密なはず(よくわからんけど)。

  • min_bin_freq

 これも上のと関係がありそう。わかるようなわからないような感じ。

  • cluster_all

 すべての点をクラスタリングするかどうか。default=Trueなので敢えてFalseにする理由は・・・(高速化なんだろうな)。

  • n_jobs : int

 いつもの並列化数

 本質的な挙動に関わるのはbandwidthで、あとは高速化のために計算を端折るための引数がいっぱいあるだけっぽいです。

 そしてbandwidthも勝手に推定してくれるので、敢えて指定する必要性を感じません(推定の良し悪しがどうかという議論はありますが)。

 今回はn_jobs以外すべてデフォルトでやってみます。

実験

 iris, wineで見てみる。せっかくなのでK-Meansと比較してみます。ついでに、入力をスケーリングすると結果が変わるかも見ます。

プログラム

 いろいろ手抜きをしています。が、とにかく結果は出ます。

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris, load_wine
from sklearn.cluster import MeanShift, KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

def main():
    iris = load_iris()
    wine = load_wine()

    kmeans = KMeans(n_clusters=3, n_jobs=-1)
    mean_shift = MeanShift(n_jobs=-1)
    s_kmeans = Pipeline([("scaler", StandardScaler()), 
                         ("kmeans", KMeans(n_clusters=3, n_jobs=-1))])
    s_mean_shift = Pipeline([("scaler", StandardScaler()), 
                             ("meanshift", MeanShift(n_jobs=-1))])


    # iris and wine
    pca = PCA(n_components=2)
    for dataset, dataset_name in zip([iris, wine], ["iris", "wine"]):
        fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(20,8))
        axes = axes.ravel()

        PCA_X = pca.fit_transform(dataset.data)
        origin_y = dataset.target
        km_y = kmeans.fit_predict(dataset.data)
        ms_y = mean_shift.fit_predict(dataset.data)
        s_km_y = s_kmeans.fit_predict(dataset.data)
        s_ms_y = s_mean_shift.fit_predict(dataset.data)
        n_clusters = [3, 3, mean_shift.cluster_centers_.shape[0],
                      3, s_mean_shift.named_steps.meanshift.cluster_centers_.shape[0]]

        for i, (y, name, n_cluster) in enumerate(
                zip([origin_y, km_y, ms_y, s_km_y, s_ms_y], 
                    ["original", "k-means", "mean-shift",
                     "scaling+k-means", "scaling+mean-shift"],
                    n_clusters)):

            for target in range(n_cluster):
                axes[i].scatter(PCA_X[y==target, 0],
                                PCA_X[y==target, 1],
                                c="rgb"[target])
            axes[i].set_title("{0} n_cluster:{1}".format(name, n_cluster))
        plt.savefig("{0}.png".format(dataset_name))

if __name__ == "__main__":
    main()

結果

 まずirisの結果から。

iris.png
iris.png

 MeanShiftは2クラスタと解釈しているようです。本来のデータとは異なりますが、敢えて人の目で見ると妥当な結果な気もします。この場合、スケーリングによる変化は微々たるものです。

 次に、wineの方。

wine.png
wine.png

 一見するとoriginalが悪すぎるように見えますが、PCAでむりやり二次元に写像しているためです。クラスタリング自体は写像前のオリジナルの空間で行っているため、影響はありません。

 この場合、合格と言って良いのは入力をスケーリングしたKMeansだけで、あとはダメダメです。特徴量の次元数が大きいから、うまく動いていないのでしょうか。ちょっと謎。

結論

 良いか悪いかの判断はつきかねますが、できることはわかりました。

 たぶんbandwidthを変えるとコロコロ結果が変わるのでしょう。どんな感じで変わるのかは、今後気が向いたときに検証しようと思います。

 →やりました。
www.haya-programming.com

【python】matplotlib.cmの使い方を説明しようと思う

テーマ

 matplotlib.cmが直接使うことはない謎の技術と思われがちなので、「普通に可愛い子なんですよ」って説明する。

cm (colormap) — Matplotlib 3.0.2 documentation

スポンサーリンク


cmとは何か

 とりあえずmatplotlib.cmの属性に何があるか見てみましょう。ちなみに、属性は組み込み関数dir()で見れます。

>>> import matplotlib.cm as cm
>>> dir(cm)
['Accent', 'Accent_r', 'Blues', 'Blues_r', 'BrBG', 'BrBG_r', 'BuGn', 'BuGn_r', 'BuPu', 'BuPu_r', 'CMRmap', 'CMRmap_r', 'Dark2', 'Dark2_r', 'GnBu', 'GnBu_r', 'Greens', 'Greens_r', 'Greys', 'Greys_r', 'LUTSIZE', 'OrRd', 'OrRd_r', 'Oranges', 'Oranges_r', 'PRGn', 'PRGn_r', 'Paired', 'Paired_r', 'Pastel1', 'Pastel1_r', 'Pastel2', 'Pastel2_r', 'PiYG', 'PiYG_r', 'PuBu', 'PuBuGn', 'PuBuGn_r', 'PuBu_r', 'PuOr', 'PuOr_r', 'PuRd', 'PuRd_r', 'Purples', 'Purples_r', 'RdBu', 'RdBu_r', 'RdGy', 'RdGy_r', 'RdPu', 'RdPu_r', 'RdYlBu', 'RdYlBu_r', 'RdYlGn', 'RdYlGn_r', 'Reds', 'Reds_r', 'ScalarMappable', 'Set1', 'Set1_r', 'Set2', 'Set2_r', 'Set3', 'Set3_r', 'Spectral', 'Spectral_r', 'Vega10', 'Vega10_r', 'Vega20', 'Vega20_r', 'Vega20b', 'Vega20b_r', 'Vega20c', 'Vega20c_r', 'Wistia', 'Wistia_r', 'YlGn', 'YlGnBu', 'YlGnBu_r', 'YlGn_r', 'YlOrBr', 'YlOrBr_r', 'YlOrRd', 'YlOrRd_r', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_deprecation_datad', '_generate_cmap', '_reverse_cmap_spec', '_reverser', 'absolute_import', 'afmhot', 'afmhot_r', 'autumn', 'autumn_r', 'binary', 'binary_r', 'bone', 'bone_r', 'brg', 'brg_r', 'bwr', 'bwr_r', 'cbook', 'cmap_d', 'cmapname', 'cmaps_listed', 'colors', 'cool', 'cool_r', 'coolwarm', 'coolwarm_r', 'copper', 'copper_r', 'cubehelix', 'cubehelix_r', 'datad', 'division', 'flag', 'flag_r', 'get_cmap', 'gist_earth', 'gist_earth_r', 'gist_gray', 'gist_gray_r', 'gist_heat', 'gist_heat_r', 'gist_ncar', 'gist_ncar_r', 'gist_rainbow', 'gist_rainbow_r', 'gist_stern', 'gist_stern_r', 'gist_yarg', 'gist_yarg_r', 'gnuplot', 'gnuplot2', 'gnuplot2_r', 'gnuplot_r', 'gray', 'gray_r', 'hot', 'hot_r', 'hsv', 'hsv_r', 'inferno', 'inferno_r', 'jet', 'jet_r', 'ma', 'magma', 'magma_r', 'mpl', 'nipy_spectral', 'nipy_spectral_r', 'np', 'ocean', 'ocean_r', 'os', 'pink', 'pink_r', 'plasma', 'plasma_r', 'print_function', 'prism', 'prism_r', 'rainbow', 'rainbow_r', 'register_cmap', 'revcmap', 'seismic', 'seismic_r', 'six', 'spectral', 'spectral_r', 'spring', 'spring_r', 'summer', 'summer_r', 'tab10', 'tab10_r', 'tab20', 'tab20_r', 'tab20b', 'tab20b_r', 'tab20c', 'tab20c_r', 'terrain', 'terrain_r', 'unicode_literals', 'viridis', 'viridis_r', 'winter', 'winter_r']

 「アホ・・・」と思われた方もいるかもしれませんが、実際にこうなっています。すべてのカラーマップに対応する属性があります。

 ではカラーマップ自体と、その型を見ましょう。

>>> cm.Paired
<matplotlib.colors.ListedColormap object at 0x7fea0de2b2b0>
>>> type(cm.Paired)
<class 'matplotlib.colors.ListedColormap'>

 あまり参考にならなかったですね・・・でも実はこれのhelpには深いことが書いてあります。

>>> help(cm.Paired)
# 中略
 |  __call__(self, X, alpha=None, bytes=False)
 |      Parameters
 |      ----------
 |      X : scalar, ndarray
 |          The data value(s) to convert to RGBA.
 |          For floats, X should be in the interval ``[0.0, 1.0]`` to
 |          return the RGBA values ``X*100`` percent along the Colormap line.
 |          For integers, X should be in the interval ``[0, Colormap.N)`` to
 |          return RGBA values *indexed* from the Colormap with index ``X``.
 |      alpha : float, None
 |          Alpha must be a scalar between 0 and 1, or None.
 |      bytes : bool
 |          If False (default), the returned RGBA values will be floats in the
 |          interval ``[0, 1]`` otherwise they will be uint8s in the interval
 |          ``[0, 255]``.
 |      
 |      Returns
 |      -------
 |      Tuple of RGBA values if X is scalar, othewise an array of
 |      RGBA values with a shape of ``X.shape + (4, )``.

 __call__というのは、インスタンスを関数っぽく呼び出したときに呼ばれるメソッドです。そしてそれぞれのカラーマップは、クラスではなくオブジェクト(インスタンス)としてmatplotlib.cmの属性になっています。つまり呼ぶことができます。

 __call__自体の内容は、0.0から1.0までの値を引数に取り、RGBAのタプル等に変換して返すというものです。配列もnumpyのユニバーサル関数と同様に受け付けます。

 ということは? こうやって使えます。

>>> cm.Paired(0)
(0.6509803921568628, 0.807843137254902, 0.8901960784313725, 1.0)
>>> cm.Paired(1)
(0.12156862745098039, 0.47058823529411764, 0.7058823529411765, 1.0)
>>> cm.Paired([0,0.5,1])
array([[0.65098039, 0.80784314, 0.89019608, 1.        ],
       [0.99215686, 0.74901961, 0.43529412, 1.        ],
       [0.69411765, 0.34901961, 0.15686275, 1.        ]])

 もうお分かり頂けましたね。この機能は、その気になればmatplotlib.pyplotと組み合わせなくても使えるようなものになっています。独立性があるのです。

スポンサーリンク


実践編

 とりあえず乱数列を生成し、適当にplotして画像として保存してみます。

import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt

np.random.seed(0)
a = np.random.random((100,100))*100

plt.figure()
plt.imshow(a, cmap=cm.Paired)
plt.savefig("a.png")

a.png
a.png

 この結果自体は特に意味とかはないんですが。ここでどんな処理が走っているのかを追いかけて、matplotlibのplot機能を使わずに同じ画像を再現してみようと思います。

 とりあえず、入力の配列はcm.Pairedに渡ってRGBAに変換されています。それはそれで良いとして、入力の配列は(概ね)0~100のスケールになっているはずです。カラーマップのインスタンスの引数には0~1を指定してあげる必要があると、上で理解しました。誰がスケーリングしているのでしょう?

 愚直にmatplotlib.cmのドキュメントを読みに行きます。

cm (colormap) — Matplotlib 3.0.2 documentation

norm : matplotlib.colors.Normalize instance
The normalizing object which scales data, typically into the interval [0, 1]. If None, norm defaults to a colors.Normalize object which initializes its scaling based on the first data processed.

 わかりました。まあ、違うクラス(ScalarMappable)の引数っぽいんですが、とにかくこれを使っておけば良いのね。

matplotlib.colors.Normalize — Matplotlib 3.0.2 documentation

 まあ、これは適当に眺めておきます。難しいこと何も書いてないし。

 ちょっと色々やってみましょう。

>>> from matplotlib.colors import Normalize
>>> import matplotlib.cm as cm
>>> Normalize()([100.0, 30.0, 0.0])
masked_array(data=[1. , 0.3, 0. ],
             mask=False,
       fill_value=1e+20)
>>> cm.Paired(Normalize()([100.0, 30.0, 0.0]))
array([[0.69411765, 0.34901961, 0.15686275, 1.        ],
       [0.2       , 0.62745098, 0.17254902, 1.        ],
       [0.65098039, 0.80784314, 0.89019608, 1.        ]])
>>> sm = cm.ScalarMappable(norm=None, cmap=cm.Paired)
>>> sm.to_rgba([100.0, 30.0, 0.0])
array([[0.69411765, 0.34901961, 0.15686275, 1.        ],
       [0.2       , 0.62745098, 0.17254902, 1.        ],
       [0.65098039, 0.80784314, 0.89019608, 1.        ]])

 まあ、使えそうという感触があります。

 では、上と同じ画像をmatplotlib.pyplotに依存しないで生成し、画像ファイルとして保存してみましょう。pltの代わりにPILを使います。

import numpy as np
import matplotlib.cm as cm
from PIL import Image

np.random.seed(0)
a = np.random.random((100,100))*100

sm = cm.ScalarMappable(norm=None, cmap=cm.Paired)
im_array = sm.to_rgba(a, bytes=True)
Image.fromarray(im_array).save("b.png")

b.png
b.png

 できました。このように使えます(枠とか画像サイズとかは勘弁してください・・・)。

まとめ

 この記事みたいなことを実際にやる必要があるか? というと微妙ですが、matplotlib.cmの実装自体はわかりやすいものだ、ということを理解してもらえたと思います。

 皆さんも、matplotlib.cmのいろいろな可愛がり方を工夫してあげてはいかがでしょうか*1

*1:たとえばこういう風に使っている人がいて、かっこいいなーと思ったりする訳です: matplotlibで色をグラデーション的に選択 - Qiita

【python】複数のin演算子を一つにまとめる方法

はじめに

 こういう状況を考える。

>>> s = "hoge! fuga! piyo!"
>>> if "hoge" in s and "fuga" in s and "piyo" in s:
...     print("fizz!")
... 
fizz!

 文字列の中に部分文字列が含まれているかを判定する、という状況で、ただ判定したい条件が複数ある。

 壮絶にまどろっこしい。できればこんな風に書きたい。

>>> if ("hoge", "fuga", "piyo") in s:

 ただ、こんな文法はpythonにはない(というかstrの__contains__メソッドあたりがこういう引数を想定していない)。

 これもだめ。

>>> if "hoge" or "fuga" or "piyo" in s:

 なんとなくいけそうな気もするけど、実は"hoge" or "fuga" or "piyo"が先に評価されて

>>> "hoge" or "fuga" or "piyo"
'hoge'

 となる。pythonの論理演算子は、基本的に受け取ったオブジェクトを返すので、こういう挙動になってしまう。

Python の or と and 演算子の罠 - Qiita

 最終的な結果は"hoge" in sなのでまったく正しくない。

 なので、これに近づける方法を考える。

スポンサーリンク


all()とany()

 どちらも組み込み関数である。all()は論理積、any()は論理和に対応する。

 実行例はこんな感じ。

>>> all([False,False])
False
>>> all([False,True])
False
>>> all([True,True])
True
>>> any([False,False])
False
>>> any([False,True])
True

 これを踏まえ、冒頭のコードを以下のように書き換えてみる。

>>> s = "hoge! fuga! piyo!"
>>> if all(x in s for x in ("hoge", "fuga", "piyo")):
...     print("fizz!")
... 
fizz!

 まだちょっとダルい。mapでも書いてみる。

>>> s = "hoge! fuga! piyo!"
>>> ins_f = lambda x:x in s
>>> if all(map(ins_f, ("hoge", "fuga", "piyo"))):
...     print("fizz!")
... 
fizz!

 これはエレガントな感じがする。一行増えているけど。一行増やしたくなければ黒魔術を使える。

>>> s = "hoge! fuga! piyo!"
>>> if all(map(s.__contains__, ("hoge", "fuga", "piyo"))):
...     print("fizz!")
... 
fizz!

 anyもまったく同様にやるだけなので、説明は省略。

注意

 ドキュメントを読むとわかるのだが、all(), any()は機能的にはとても単純なようだ。

組み込み関数 — Python 3.7.4 ドキュメント

 それは別に良いのだが、先に引数をすべて渡す必要があるせいで、リストを引数に渡した場合、短絡評価されない(というか短絡評価する前に計算してしまっているので無意味)という問題がある。

 all(), any()を使うときは引数をジェネレータ式などにしてあげた方が良いだろう。map objectもイテレータで、同様の性質があるので、問題ない。

www.haya-programming.com

まとめ

 この記事の方法を使うと、複数のin演算子が入るような場合にまとめて簡潔に書くことができる。

 あまり使う機会は多くないかもしれないが、参考にしてほしい。

【python】numpyで乱数のseedを設定する方法

 「seed(種)」とか「random state」とか呼ばれる奴の設定方法。これを設定することで、乱数の処理に再現性を与えることができる。

方法

 np.random.seed()を呼ぶと、とりあえずseedが引数でリセットされる。

https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html

 やってみる。

>>> import numpy as np
>>> np.random.seed(0)  # seedを設定
>>> np.random.random((2,3))  # 一様乱数の配列を作ってみる
array([[0.5488135 , 0.71518937, 0.60276338],
       [0.54488318, 0.4236548 , 0.64589411]])
>>> np.random.random((2,3))  # 同じにはならない
array([[0.43758721, 0.891773  , 0.96366276],
       [0.38344152, 0.79172504, 0.52889492]])
>>> np.random.seed(0)  # もう一度設定する
>>> np.random.random((2,3))  # 最初と同じになる
array([[0.5488135 , 0.71518937, 0.60276338],
       [0.54488318, 0.4236548 , 0.64589411]])

 まあ、理解はできる。

スポンサーリンク


 ところで、np.random.seed()のドキュメントにはこう書いてあった。

This method is called when RandomState is initialized. It can be called again to re-seed the generator. For details, see RandomState.

 RandomStateというクラスが内部では動いているということなのだろうか。

https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html

Parameters:
seed : {None, int, array_like}, optional

Random seed used to initialize the pseudo-random number generator. Can be any integer between 0 and 2**32 - 1 inclusive, an array (or other sequence) of such integers, or None (the default). If seed is None, then RandomState will try to read data from /dev/urandom (or the Windows analogue) if available or seed from the clock otherwise.

 よくわからないが、乱数の状態のクラスなので(直訳)、たぶんインスタンス化してよろしく使えるのだと思う。

>>> rand_state = np.random.RandomState(0)
>>> rand_state.random_sample((2,3))  # なぜかrandomメソッドがなくてrandom_sampleメソッドが同じ機能を持っている・・・
array([[0.5488135 , 0.71518937, 0.60276338],
       [0.54488318, 0.4236548 , 0.64589411]])
>>> rand_state.random_sample((2,3))
array([[0.43758721, 0.891773  , 0.96366276],
       [0.38344152, 0.79172504, 0.52889492]])

 できた。それでえっと、何に便利かなこれは・・・。

 まあ、常識的な用途ではnp.random.seed()があれば困らないと思う。わざわざオブジェクト化して取り扱いたいような状況は、たぶんあるにはあるだろうけど、ライブラリとかアルゴリズムのコアなコード書いてるとき以外はほぼないだろう。

まとめ

 簡単にできる。

【python】operator.itemgetterを使うべきか否か問題

はじめに

 この記事を開いた人の大半は「itemgetter? なにそれ」という反応でしょう。
 (いや、検索で来た人はそうでもないかもしれないけど)

 itemgetterは以下のように使えるものです。

>>> lst = list(zip([1,3,5,6,7,1,4], [3,4,1,0,8,5,2]))  # 特に値に意味はない
>>> sorted(lst, key=lambda x:x[0])
[(1, 3), (1, 5), (3, 4), (4, 2), (5, 1), (6, 0), (7, 8)]
>>> sorted(lst, key=lambda x:x[1])
[(6, 0), (5, 1), (4, 2), (1, 3), (3, 4), (1, 5), (7, 8)]
>>> from operator import itemgetter
>>> sorted(lst, key=itemgetter(0))
[(1, 3), (1, 5), (3, 4), (4, 2), (5, 1), (6, 0), (7, 8)]
>>> sorted(lst, key=itemgetter(1))
[(6, 0), (5, 1), (4, 2), (1, 3), (3, 4), (1, 5), (7, 8)]

 つまり、itemgetter(0)はlambda x:x[0]と等価です。

 参考:
 10.3. operator — 関数形式の標準演算子 — Python 3.6.5 ドキュメント

 では、どちらを使うべきなのでしょうか?



 目次

可読性とタイプ数

可読性

 可読性は正直なんとも言い難いです。

 とりあえずlambda式はlambda式の概念さえ理解していればわかるともいえますし、それでもわかりづらいという人もいるでしょう。また、見た目はグロいです。

 itemgetterは見た目はすっきりしています。これはとても大切なことで、lambdaまみれで入り組んだコードはもううんざりです。ただ、itemgetter自体を知らない人は恐らく多いので、初見では「helpを見る」「ググる」という作業が発生します。また、何をやっているか明示的にわかりづらいので、lambdaの方がマシという考え方もあります。

タイプ数

 タイプ数的には、はじめにを見て頂ければわかる通り互角です。itemgetterという途方もなく長い識別子のせいで、lambda式に勝てません。

 一応、asで別名にしてしまうという反則技(?)はあります。

from operator import itemgetter as get  # もうちょっと良い名前を考えて

 こうするとタイプ数的にはlambdaに勝てそうです。ただ、こうするのは如何なものか? と相成ります。混乱を招きそうなので、あまりやりたくない手です。

結論

 可読性・タイプ数では決められない。

速度差

 まさかこんなのに速度差なんてねーだろ、そう思っていた時期が私にもありました。

 都合によりminで比較。こうすると常に全数線形探索されるはずなので、再現性等が良いです。

>>> import timeit
>>> lst = list(zip(range(10000), range(10000)))
>>> timeit.timeit(lambda :min(lst, key=lambda x:x[0]), number=1000)
1.33255122898845
>>> timeit.timeit(lambda :min(lst, key=itemgetter(0)), number=1000)
0.6424821490072645

 倍違うのだった。えぇ・・・

 怪しいので、key関数を外部に定義してみます。

>>> key1 = lambda x:x[0]
>>> key2 = itemgetter(0)
>>> timeit.timeit(lambda :min(lst, key=key1), number=1000)
1.32674505902105
>>> timeit.timeit(lambda :min(lst, key=key2), number=1000)
0.6525059329869691

 変わらず。

結論

 速いので、特にこだわりがなければitemgetterを使いましょう。

 私はわざわざimportするのがダルいのと、個人的にlambdaがそれなりに好きという理由でlambdaを使い続けると思います。競プロとかのときはitemgetter有利なんじゃないでしょうか。

【python】np.matrixの速度を測る

 numpyで行列演算を行う方法としては、普通のnumpy配列に行列演算系の関数を適用していく方法と、あまり知られていないがnp.matrix型やnp.mat型を使う方法がある。

 速度が違ったりするのだろうか? 仮に違うと困る(というか場面によって適切な方を選ぶ必要が生じてくる)ので、こんなコードを書いてみた。

import timeit

import numpy as np 

def main():
    a = np.random.random((100, 100))
    b = np.random.random((100, 100))
    a_m = np.matrix(a)
    b_m = np.matrix(b)
    a_mt = np.mat(a)
    b_mt = np.mat(b)
   
    print("add")
    print(np.allclose(a+b, np.array(a_m+b_m), np.array(a_mt+b_mt)))
    print(timeit.timeit(lambda : a+b, number=100000))
    print(timeit.timeit(lambda : a_m+b_m, number=10000))    
    print(timeit.timeit(lambda : a_mt+b_mt, number=10000))    

    print("mul")
    print(np.allclose(a*b, 
                      np.array(np.multiply(a_m, b_m)),
                      np.array(np.multiply(a_mt, b_mt))))
    print(timeit.timeit(lambda : a+b, number=10000))
    print(timeit.timeit(lambda : np.multiply(a_m, b_m), number=10000))        
    print(timeit.timeit(lambda : np.multiply(a_mt, b_mt), number=10000))        
    
    print("dot")
    print(np.allclose(np.dot(a,b),
                      np.array(a_m*b_m),
                      np.array(a_mt*b_mt)))
    print(timeit.timeit(lambda : np.dot(a,b), number=10000))
    print(timeit.timeit(lambda : a_m*b_m, number=10000))
    print(timeit.timeit(lambda : a_mt*b_mt, number=10000))

if __name__ == "__main__":
    main()

 100*100の行列の加算、要素積、行列積の計算を、np.arrayを使う方法とnp.matrix、np.matを使う方法で書いて、速度を比較している(ついでに計算結果が同じになるかも確認しているけど、これは念のためにやっているだけで、無論Trueになる)。

 結果は以下の通り。

add
True
0.6425584320095368
0.1307114879891742
0.1383463980164379
mul
True
0.07879749499261379
0.13494228100171313
0.1337542689871043
dot
True
0.8895792970142793
0.9789612090098672
1.0360321149928495

 変わるのか(困惑)。一瞬計測誤差かと思って何回も回したけど、傾向は変わらず。add, mulはnp.arrayの方が速い。dotはほとんど速度差はないようだが。

 何かの間違いだろうと思って検索したら、こんなスタックオーバーフローの質問が出てきてしまった。

python - numpy np.array versus np.matrix (performance) - Stack Overflow

 np.arrayの方が速い、と言われてしまっている。こうなると認めるしかないようだ。

 こうなると、それなりに速度重視の場面ではnp.matrixは使えないですねぇ。かといってnp.arrayで行列演算は苦行(Cとかフォートランでforループ回して書くのに比べればずっとマシだろうけど)。

 嫌な現実を見せられてしまった感がある。

 2018年11月15日追記:
 そうこうしているうちに、numpy 1.15.0から「no longer recommended」になってしまいました。

numpy.matrix — NumPy v1.15 Manual

 配列型をいろいろな関数と組み合わせて行列演算してくれということなのでしょう。

【python】numpy.meshgridの基本的な使い方まとめ

はじめに

 numpyのmeshgridの使い方があまり理解できなくて、なんとなくコピペで動かしていたので、どのようなものなのかまとめておくことにしました。

 目次

とりあえずmeshgridを作ってみる

 meshgrid自体はこのように何の変哲もないものです。

>>> import numpy as np
>>> a = np.arange(0,3,0.5)
>>> b = np.arange(0,10,2)
>>> a
array([0. , 0.5, 1. , 1.5, 2. , 2.5])
>>> b
array([0, 2, 4, 6, 8])
>>> X, Y = np.meshgrid(a, b)
>>> X
array([[0. , 0.5, 1. , 1.5, 2. , 2.5],
       [0. , 0.5, 1. , 1.5, 2. , 2.5],
       [0. , 0.5, 1. , 1.5, 2. , 2.5],
       [0. , 0.5, 1. , 1.5, 2. , 2.5],
       [0. , 0.5, 1. , 1.5, 2. , 2.5]])
>>> Y
array([[0, 0, 0, 0, 0, 0],
       [2, 2, 2, 2, 2, 2],
       [4, 4, 4, 4, 4, 4],
       [6, 6, 6, 6, 6, 6],
       [8, 8, 8, 8, 8, 8]])

 つまり(0,0)の値を(X[0,0],Y[0,0])として表現できる。ここまでは常識的な話。

計算する

 numpy配列なので、配列全体でそのまま計算することができます。numpyでいうところの[1, 2, 3] + [4, 5, 6]が[5, 7, 9]になるのと同じです(この表記はリストのリテラルなのでできませんが)。

>>> Z = np.sin(X+Y)
>>> Z
array([[ 0.        ,  0.47942554,  0.84147098,  0.99749499,  0.90929743,
         0.59847214],
       [ 0.90929743,  0.59847214,  0.14112001, -0.35078323, -0.7568025 ,
        -0.97753012],
       [-0.7568025 , -0.97753012, -0.95892427, -0.70554033, -0.2794155 ,
         0.21511999],
       [-0.2794155 ,  0.21511999,  0.6569866 ,  0.93799998,  0.98935825,
         0.79848711],
       [ 0.98935825,  0.79848711,  0.41211849, -0.07515112, -0.54402111,
        -0.87969576]])

 ま、これはmeshgridの機能というより、numpyの機能。

plotしてみる

 ぶっちゃけ3Dのplot以外でmeshgrid使うことってあるんですか・・・?

>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>> fig = plt.figure()
>>> ax = Axes3D(fig)
>>> ax.plot_wireframe(X, Y, Z)
<mpl_toolkits.mplot3d.art3d.Line3DCollection object at 0x7f04d0367860>
>>> plt.show()

meshgridを使って三次元プロット
meshgridを使って三次元プロット

 XとYを足してsinに通すという謎の関数にしてしまったので見た目は気持ち悪いですが、それはさておきうまくプロットすることができています。matplotlibには同じような3次元データのプロット用の関数がたくさんあり、だいたいmeshgridを受け付けるので、いろいろな方法でプロットすることができます。

xyz座標の配列に変換する

 meshgridからx,y,zの座標列に変換したい、ということもあるでしょう。あまり難しく考える必要はなくて、Xをflattenすればx座標列が、Yをflattenすれば……というように座標を得ることができ、あとは適当に結合すれば任意の形にできます。

>>> x = X.ravel()
>>> y = Y.ravel()
>>> z = Z.ravel()
>>> x.shape
(30,)
>>> y.shape
(30,)
>>> z.shape
(30,)
>>> x
array([0. , 0.5, 1. , 1.5, 2. , 2.5, 0. , 0.5, 1. , 1.5, 2. , 2.5, 0. ,
       0.5, 1. , 1.5, 2. , 2.5, 0. , 0.5, 1. , 1.5, 2. , 2.5, 0. , 0.5,
       1. , 1.5, 2. , 2.5])
>>> y
array([0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 6, 6, 6, 6,
       6, 6, 8, 8, 8, 8, 8, 8])
>>> z
array([ 0.        ,  0.47942554,  0.84147098,  0.99749499,  0.90929743,
        0.59847214,  0.90929743,  0.59847214,  0.14112001, -0.35078323,
       -0.7568025 , -0.97753012, -0.7568025 , -0.97753012, -0.95892427,
       -0.70554033, -0.2794155 ,  0.21511999, -0.2794155 ,  0.21511999,
        0.6569866 ,  0.93799998,  0.98935825,  0.79848711,  0.98935825,
        0.79848711,  0.41211849, -0.07515112, -0.54402111, -0.87969576])

 こうすると色々な関数に渡せそうで便利。結合については、こちらの記事を参照してください。

www.haya-programming.com

 たとえばxyzで(n_samples, 3)にしたければnp.stack([x, y, z], axis=1)などでよさそうです。

逆にxyz座標からmeshgridを生成する

 x, y, zの座標データが与えられていて、プロットするためにmeshgridに変換したいというシチュエーションもあります。

 この場合は元のデータが都合よくグリッド状になっていないことが多いので、補完という処理をした上で変換することが可能です。

 詳細はこちらの記事を御覧ください。

www.haya-programming.com

まとめ

 いろいろ動かしてみると、それなりにやっていることが単純なのがわかって、親近感が湧いてきました。とにかく三次元プロットでは必需品なので、meshgridは使いこなそうということです。

【python】sklearnのclass_weightの挙動

はじめに

 先に断っておくと、class_weightの挙動はモデルによって異なる可能性が十分ある。今回はsklearn.svm.SVCとsklearn.ensemble.RandomForestClassifierのドキュメントを参照して、一応基本的に共通する部分を抜き出した。

 class_weightを調整する必要が出てきたときは、自分が使うモデルで確認してください。

 参考:
3.2.4.3.1. sklearn.ensemble.RandomForestClassifier — scikit-learn 0.20.1 documentation
sklearn.svm.SVC — scikit-learn 0.20.1 documentation

解説

 与えられる引数はNone(デフォルト)、辞書、"balanced"という文字列の3種類というのが基本だと思う*1

 Noneを渡した場合、class_weightはすべてのクラスに対して1であると仮定される。これは悩むような話ではない。

 辞書を渡す場合、キーはクラスラベル、値は重みになる。すべてのクラスラベルと対応する重みを辞書の要素に与えなくても動くようだが、まあ与えた方が無難。詳しく検証はしていない。

 重みは大きいほどそのクラスを重視する意味になる。なので、たとえばクラス0に100件、クラス1に1000件のデータがあったとしたら、何もしないとクラス1が重視されすぎてしまう懸念がある。なので、クラス1を1/10してやるか、クラス0を10倍してやれば同じ比率になるだろう、という気がする。

{0:1000/100,1:1}
{0:1, 1:100/1000}

 直感的にはこんな感じで良い。

 "balanced"を指定すると似たような演算をしてくれる。ただし、ちょっとやっていることが異なる。公式によると、この演算が行われるらしい。

n_samples / (n_classes * np.bincount(y))

 n_samplesは全サンプル数の合計(恐らくfit時のサンプルである)。n_classesはクラス数だが、np.bincountというのが見慣れない処理である。でもやることは簡単。

>>> import numpy as np
>>> np.bincount([0,0,0,1,1,1,1,2,2,2,2,3,4,5,6,5,4,4,4,3,2,1])
array([3, 5, 5, 2, 4, 2, 1])

 つまりクラスラベルごとに出現回数を数えている(index=0の要素が0の出現回数・・・という仕組み)。これの逆数を取って、あとは定数をかけているだけの式と解釈できる。そこそこ無難そうな式ではある。

まとめ

 普通は意識しないと思います・・・。

 クラスごとのサンプル数の偏りが無視できそうな程度に小さかったり、無視できそうな程度にサンプルが大量にある場合はNone(デフォルトのまま)、

 偏っている場合でも特殊なケース以外は"balanced"で良いと思う。

 これをパラメタチューニングしてどうこうとかは、考えない方が良いんじゃないかなぁ・・・。

*1:ランダムフォレストは複数値への分類に対応しているのでlistも受け取るようだ。また、ランダムフォレストはブートストラップの絡みで“balanced_subsample”というのも指定できるようだ。これらは特殊だと思うので説明しない(ドキュメントに書いてあるし)

【python】sklearnで「何もしない」モデルがほしい

 sklearnで「何もしない」モデルがあると、チョー便利。個人的にはそう思う。

 どうやって使うかというと、具体的には以前の記事で書いたFeatureUnionと組み合わせて使う。

 参考(以前の記事):【python】複数の特徴をまとめるFeatureUnion - 静かなる名辞

 たとえば、100次元の特徴量があったとき、「入力をそのまま出力する(何もしない)モデル」と「何らかの加工をして出力するモデル」を組み合わせたモデルに入力すれば、100次元+加工した特徴量の次元の新たな特徴量とすることができる。これが有効な場面はそれなりに多そう。

 ・・・なのだが、sklearnにはそれに当てはまるようなモデルは用意されてなさそうだし、ググっても出てこない。自分で書いている人はいた。

from sklearn.base import BaseEstimator, TransformerMixin

class IdentityTransformer(BaseEstimator, TransformerMixin):
    def __init__(self):
        pass
    
    def fit(self, input_array, y=None):
        return self
    
    def transform(self, input_array, y=None):
        return input_array*1

 引用元:sklearn Identity-transformer – Laurent H. – Medium

 原文通りだけど*1は要らないと思う。とにかく、同じようなことを考える人はいるものだ。簡単に書けるのもわかったけど、そのために自分でクラスを定義するのもなんだかなぁ、という気がする。もっと色々なメソッドに対応させようとすると面倒臭さが増していきそうなのも、あまりよくない。

 そこで改めてドキュメントを漁ったら、使えそうなのがあった。FunctionTransformerなるもの。

sklearn.preprocessing.FunctionTransformer — scikit-learn 0.20.1 documentation

 任意の関数を渡してTransformerを生成できるというなかなかのスグレモノである。こいつは関数が渡されないとき、入力をそのまま出力するとドキュメントに書いてある。使えそう。簡単に確認する。

>>> from sklearn.datasets import load_iris
>>> from sklearn.preprocessing import FunctionTransformer
>>> iris = load_iris()
>>> identity_transformer = FunctionTransformer()
>>> iris.data[:10]
array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2],
       [5.4, 3.9, 1.7, 0.4],
       [4.6, 3.4, 1.4, 0.3],
       [5. , 3.4, 1.5, 0.2],
       [4.4, 2.9, 1.4, 0.2],
       [4.9, 3.1, 1.5, 0.1]])
>>> identity_transformer.fit_transform(iris.data[:10])
array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2],
       [5.4, 3.9, 1.7, 0.4],
       [4.6, 3.4, 1.4, 0.3],
       [5. , 3.4, 1.5, 0.2],
       [4.4, 2.9, 1.4, 0.2],
       [4.9, 3.1, 1.5, 0.1]])

 使えた。

 本来の使い方ではないと思うけど、これがベター・・・なのかなぁ。もっと良いやり方を知っている方がいたら、教えてください。

【python】複数の特徴をまとめるFeatureUnion

 単一の入力データから、複数の処理方法で幾つもの異なる特徴量が得られる・・・というシチュエーションがある。

 この場合、「どれが最善か」という観点でどれか一つを選ぶこともできるけど、そうすると他の特徴量の情報は捨ててしまうことになる。総合的な性能では他に一歩譲るが、有用な情報が含まれている特徴量がある・・・というような場合は、ちょっと困る。

 こういう状況で役に立つのがFeatureUnion。特徴抽出や次元削減などのモデルを複数まとめることができる。

 結果はConcatenateされる。Concatenateというのがわかりづらい人もいると思うけど、たとえば手法1で10次元、手法2で20次元の特徴量ベクトルが得られたら、これをそのまま横に繋げて30次元のベクトルとして扱うということ。

sklearn.pipeline.FeatureUnion — scikit-learn 0.20.1 documentation

 ちなみに、こいつはsklearn.pipeline以下に存在する。Pipelineの兄弟みたいな扱い。引数の渡し方とかもほとんど同じである。

 簡単に試してみよう。digitsの分類を行うことにする。PCA+GaussianNB, LDA+GNB, FeatureUnion(PCA, LDA)+GNBの3パターンでスコアを見比べる。

import warnings
warnings.filterwarnings('ignore')

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.naive_bayes import GaussianNB
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.model_selection import cross_validate, StratifiedKFold

def main():
    digits = load_digits()
    
    pca = PCA(n_components=30)
    lda = LDA()
    gnb = GaussianNB()
    
    pca_gnb = Pipeline([("pca", pca), ("gnb", gnb)])
    lda_gnb = Pipeline([("lda", lda), ("gnb", gnb)])
    pca_lda_gnb = Pipeline([("reduction", FeatureUnion([("pca", pca),
                                                        ("lda", lda)])),
                            ("gnb", gnb)])

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

    for name, model in zip(["pca_gnb", "lda_gnb", "pca_lda_gnb"], 
                           [pca_gnb, lda_gnb, pca_lda_gnb]):

        skf = StratifiedKFold(shuffle=True, random_state=0)
        scores = cross_validate(model, digits.data, digits.target,
                                cv=skf, scoring=scoring)
        
        p = scores["test_p"].mean()
        r = scores["test_r"].mean()
        f = scores["test_f"].mean()
        print(name)
        print("precision:{0:.3f} recall:{1:.3f} f1:{2:.3f}".format(p,r,f))

if __name__ == "__main__":
    main()

 結果は、

pca_gnb
precision:0.947 recall:0.944 f1:0.945
lda_gnb
precision:0.955 recall:0.953 f1:0.953
pca_lda_gnb
precision:0.959 recall:0.957 f1:0.957

 ちょっと微妙だけど、誤差ではないみたい。このように比較的手軽に性能を改善できることがわかる(効くかどうかはケースバイケースだけど)。

有意水準5%の論文が100本あったら

 この記事は思いついたままに書いたポエム。


 有意水準5%とは、その判断(主張)の妥当性が95%である、ということを意味する。

 よって、有意水準5%で検定したら、100回に5回は第1種の過誤を犯す。

 有意水準5%の論文が100本あったら、(いちおうすべての論文が正しいプロセスを踏んでいると仮定しても)うち5本は間違っている。


 恐らく現実の論文はそんなに酷いことにはなっていないと思う(ただし、100本あったらそもそもプロセスが正しくないものは一定数入ってくるだろうけど)。これがどういうことなのかというと、

  • そもそも5%なんて甘い有意水準は使っていない(これはあると思うけど、とりあえず無視することにする)
  • 最初からある程度有力な仮説を立てて検証しているので、95%に「仮説の妥当性」がかかってくると考えられる(ちょっと異論もあるかもしれないが、「ある仮説を妥当だと思って検証し、けっきょく妥当ではなかったという事象の可能性」を考えるとけっきょく効いてくると思われる)
  • 論文に発表した以外にも色々実験をしたりして、妥当性を判断している。一回検定しただけ、というのは考えづらい(逆に言えば、再現が難しい話(世の中から取ってきた統計量をそのまま使うような奴)だと5%は割とそのまま5%かもしれないので、注意が必要)

 恐らくこのような事情が絡んでくるので、有意水準5%の妥当性は実際には98%とか99%とか、それぐらいには信頼できるように(個人的には)思える*1

 

*1:だからって98%とか99%でも割ときつい水準だと思うので、これ以上緩める理由はないと思うが

【python】1つおきにリスト・文字列などから抽出する

 スライスの基本的な話なんだけど、意外と知らない人が多いと思うので。

 スライスで1つおきに取り出すには、こうする。

>>> "hogehoge~"[::2]
'hghg~'

 スライスで指定できるのはstart, stop, stepであり、上のように指定するとstart, stopはNoneでstepが2になる。

 参考:【python】sliceのちょっと深イイ(かもしれない)話 - 静かなる名辞


 ではstepを1にすると? これは元の文字列が返る。

>>> "hogehoge~"[::1]
'hogehoge~'

 たまにこれをreverseに悪用(?)する人がいる。[::-1]と指定するとそういう挙動になるので。

>>> "hogehoge~"[::-1]
'~egohegoh'

 でもこれは可読性が悪い、というか完全に初見殺しなので、やめた方が良いと思う。reversedというものがあるので、そっちを使おう。

 ・・・と思って書いてみたけど、腹立たしいことに空気を読まない子(iteratorで返してくれるアレ)なので、こんな真似をする必要がある。

>>> "".join(reversed("hogehoge~"))
'~egohegoh'

 どっちがマシかなぁ・・・。

 話をもとに戻して、[::-1]でこうなるということは、これもできる。1つおきに取り出して逆転させる。

>>> "hogehoge~"[::-2]
'~ghgh'

 便利、なのだろうか(どんなときに?)。