【python】rangeではin演算子が使える。速度は微妙かも

はじめに

　今日コードを書いていて、rangeでもinが使えることに気づきました。

>>> 10 in range(20)
True

　ドキュメントを見るとシーケンス型としての機能は一通り備えているようです。

range オブジェクトは collections.abc.Sequence ABC を実装し、包含判定、要素インデックス検索、スライシングのような機能を提供し、負のインデックスをサポートします (シーケンス型 — list, tuple, range を参照):

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

　ちなみにfloatでもなんとなくTrueになりましたが、あくまでも離散値の包含で比較される雰囲気です。

>>> 10.0 in range(20)
True
>>> 10.1 in range(20)
False

測ってみる

　こういうものがあると、速度が気になります。

>>> r = range(1000)
>>> l = list(range(1000))
>>> s = set(l)
>>> import timeit
>>> timeit.timeit(lambda : 500 in r)
0.2108146829996258
>>> timeit.timeit(lambda : 500 in l)
8.18245309899794
>>> timeit.timeit(lambda : 500 in s)
0.13419233300010092
>>> timeit.timeit(lambda : 0 <= 500 < 1000)
0.12016064700219431

　同じ長さのlistよりは圧倒的に速いものの、setに負けるという結果に。恐らくシーケンスに展開して線形探索をする訳ではないものの、内部処理がそこそこ複雑なのでしょう。また、単に値の区間だけ確認したいのなら、不等式による比較が最速のようです（setが異様に速いと言うべきか・・・）。

使いどころ

　rangeということはstepも入れられるので、複雑な条件のときにはmodで書くより可読性が良いかもしれません。

>>> [x in r for x in range(10)]
[True, False, False, True, False, False, True, False, False, True]

　他に積極的に使う理由は思いつきません。

xyzの点データを内挿してmeshgridにしmatplotlibでプロットする

はじめに

　pythonでmatplotlibを使って作図するとき、三次元のデータでpcolormeshとかcontourでやるような等高線プロットを作りたいんだけど、手持ちのデータはxyzが紐付いた点のバラバラな離散データだけ……ということがままあります。

　散布図ならそれでも良いのですが、等高線などをプロットしようとするとmeshgrid的な形式のデータが要ると思います。困ります。

　なので、内挿してなんとかしてみます。

方針

　matplotlibのドキュメントを読んでいると、それらしいものを見つけました。

mlab — Matplotlib 3.0.2 documentation

matplotlib.mlab.griddata(x, y, z, xi, yi, interp='nn')[source]
Deprecated since version 2.2: The griddata function was deprecated in Matplotlib 2.2 and will be removed in 3.1. Use scipy.interpolate.griddata instead.

　scipyのscipy.interpolate.griddataを使え、ということらしいです。

scipy.interpolate.griddata — SciPy v0.18.1 Reference Guide

　ちょっと使い方がわかりづらいですが、仮に

• xy：shape=(n_samples, 2)のxy座標のデータ
• z：shape=(n_samples,)のz座標のデータ
• X, Y：xy空間のmeshgrid

　という変数を置くとすると

Z = griddata(xy, z, (X, Y))

　でZ（zのmeshgrid）が得られる、ということのようです。詳細はドキュメントのサンプルを見て確認してみてください。

実験

　二次元の正規分布を作り、真値と内挿で得られた値をpcolormesh, contourでプロットしました。

　コードを以下に示します。

import numpy as np
from scipy import stats
from scipy import interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def main():
    # data (samples)
    xy = np.random.uniform(low=-10, high=10, size=(500, 2))
    norm = stats.multivariate_normal(mean=[2.0, 3.0], cov=[[3, 1],[1,3]])
    z = norm.pdf(xy)
    
    # 3d scatter of samples
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(xy[:,0], xy[:,1], z)
    plt.savefig("scatter_3d.png")

    # mesh x-y
    x = y = np.linspace(-10, 10, 500)
    X, Y = np.meshgrid(x, y)

    # interpolated Z
    i_Z = interpolate.griddata(xy, z, (X, Y))

    # true Z
    t_Z = norm.pdf(np.vstack([X.ravel(), Y.ravel()]).T).reshape(X.shape)

    # plot
    fig, axes = plt.subplots(ncols=2)
    plt.subplots_adjust(wspace=0.4)

    # plot true Z
    im = axes[0].pcolormesh(X, Y, t_Z, cmap="jet")
    plt.colorbar(im, ax=axes[0])
    axes[0].contour(X, Y, t_Z, colors=["black"])
    axes[0].set_title("true")

    # plot interpolation Z
    im = axes[1].pcolormesh(X, Y, i_Z, cmap="jet")
    plt.colorbar(im, ax=axes[1])
    axes[1].contour(X, Y, i_Z, colors=["black"])
    axes[1].set_title("interpolation")

    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　stats.multivariate_normalについてはこちらの記事もご覧くだい。

www.haya-programming.com

　それほど難しいことはやっていないので、詳細はコードを読んでいただければわかると思います。出力される画像を以下に示します。

　まず、参考のために出力したサンプルの三次元散布図です。

scatter_3d.png

　けっこうスカスカな印象を受けます。

　次に、真値と内挿で得られた値をプロットした結果です。

result.png

　完璧とは言えませんが、そこそこ良さそうな結果が得られています。なお、内挿のグラフの外側の白い部分は内挿できなくてnanになっている領域です。griddataのfill_valueオプションなどで挙動を変えられるので、検討してみてください。

まとめ

　このように内挿してプロットすることができます。いまいちなデータしか手持ちにないとき威力を発揮します。

　ただし、元のデータの性質を歪めている側面があるので、ある程度注意が必要になります。それさえ気をつければ十分使えると思います。

【python】sklearnのFeatureAgglomerationを使ってみる

はじめに

　FeatureAgglomerationは階層的クラスタリングを用いた教師なし次元削減のモデルです。特徴量に対して階層的クラスタリングを行い（つまり通常のサンプルに対するクラスタリングと縦横の向きが入れ替わる）、似ている特徴量同士をマージします。マージの方法はデフォルトでは平均のようです。

　使用例をあまり見かけませんが、直感的な次元削減方法なので何かしらの役に立つかもしれないと思って使ってみました。

sklearn.cluster.FeatureAgglomeration — scikit-learn 0.20.1 documentation

使い方

　パラメータは以下の通り。

class sklearn.cluster.FeatureAgglomeration(
    n_clusters=2, affinity=’euclidean’, memory=None, connectivity=None, 
    compute_full_tree=’auto’, linkage=’ward’, pooling_func=<function mean>)

　色々いじれるように見えますが、主要パラメータは2つだけです。

• n_clusters

　PCAでいうところのn_componentsです。変換先の次元数を表します。

• pooling_func

　似ている特徴量をマージする方法。callableが渡せます。何もしなければ平均が使われるので、平均より気の利いた方法を思いつく人以外はそのままで大丈夫です。

　あとは階層的クラスタリングのオプションが色々あります。それはそれで大切なものだと思いますが、今回は無視することにします。

実験

　もう何番煎じかわかりませんが、irisの2次元写像で試します。

import matplotlib.pyplot as plt

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import FeatureAgglomeration

def main():
    iris = load_iris()

    pca = PCA(n_components=4)
    ss = StandardScaler()
    agg = FeatureAgglomeration(n_clusters=2)

    pca_X = pca.fit_transform(iris.data)
    agg_X = agg.fit_transform(
        ss.fit_transform(iris.data))

    print(pca.components_)
    print(agg.labels_)

    fig, axes = plt.subplots(nrows=1, ncols=2)
    axes[0].scatter(pca_X[:,0], pca_X[:,1], c=iris.target)
    axes[0].set_title("PCA")
    axes[1].scatter(agg_X[:,0], agg_X[:,1], c=iris.target)
    axes[1].set_title("FeatureAgglomeration\n{}".format(agg.labels_))
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

　動作原理、目的と用途を考えると、事前にスケーリングしておいた方が恐らく無難です。

　printされた出力。

[[ 0.36138659 -0.08452251  0.85667061  0.3582892 ]
 [ 0.65658877  0.73016143 -0.17337266 -0.07548102]]
[0 1 0 0]

　FeatureAgglomerationは圧倒的に結果の解釈性が良いことがわかります。写像先の0次元目は元の0,2,3次元目の平均で*1、写像先の1次元目は元の1次元目ですね。こういうのはシチュエーション次第ですが、ちょっと嬉しいかもしれません。

　出力される画像。

プロットの結果

　概ねPCAと同等に使えています。うまく言葉で表現はできませんが、FeatureAgglomerationの方はなんとなくギザギザ感？みたいなものがあります。平均するとそうなる、というのがなんとなくわかる気もするし、わからない気もする。

考察

　結果の解釈性が良いのと、まがりなりにすべての特徴量の情報が結果に反映されるので、PCAより使いやすいシチュエーションはあると思います。分類前の次元削減とかで使ったときの性能とかは今回検討していませんが、たぶんそんなに良いということはないはず。

　あとドキュメントをあさっていたら、こんなページがあったので、

Feature agglomeration — scikit-learn 0.20.1 documentation

　真似してPCAでも同じものを出してみたら（コードはほとんど書き換えていないので省略。agglo = の行で代入するモデルをコメントアウトで切り替えて、あとlabels_の出力を外しただけです）、やっぱりFeatureAgglomerationはヘボかった（低次元で元の情報を保持することに関しては性能が低かった）です。

　10次元に落として元の情報をどこまで復元できるかという実験。

PCA

FeatureAgglomeration

　まあ、これは仕方ないか。

まとめ

　とにかく結果の解釈性の良さを活かしたい、とか、なにか特別な理由があって使う分には良いと思います。

*1:厳密にはどれか2つが先に平均されて、更に残りと平均されるはず。つまり3つの比重が違う順番はチェックしていないのでわかりませんが、children_属性をちゃんと読み取ればわかると思います

【python】複数の条件を総なめするときの簡略化

　たとえば、こういうものを書きたいとする。

def f(a, b):
    if a == "0" and b == "0":
        print("a:0, b:0")
    elif a == "0" and b == "1":
        print("a:0, b:1")
    elif a == "1" and b == "0":
        print("a:1, b:0")
    elif a == "1" and b == "1":
        print("a:1, b:1")

　条件式がいかにも冗長。
　（ifをネストすれば良い、printも冗長、フォーマット文字列使えば良い、あるいは辞書に入れて結果の表示を分岐させれば良いというツッコミは受け付けないことにする。実際には処理内容の分岐が入ることを想定する。）

　こういうときは、シーケンスにして値比較するとスマート。文字列だと最も簡単なのはこれ。

def f(a, b):
    c = a + b
    if c == "00":
        print("a:0, b:0")
    elif c == "01":
        print("a:0, b:1")
    elif c == "10":
        print("a:1, b:0")
    elif c == "11":
        print("a:1, b:1")
||<　

　もう少し汎用的な型を想定するなら、listかtupleに入れると良い。

>|python|
def f(a, b):
    c = [a, b]
    if c == ["0", "0"]:
        print("a:0, b:0")
    elif c == ["0", "1"]:
        print("a:0, b:1")
    elif c == ["1", "0"]:
        print("a:1, b:0")
    elif c == ["1", "1"]:
        print("a:1, b:1")

　コードがシンプルになり、可読性も向上した。

ImportError: No module named 'sklearn.cross_validation'の対処

概要

　sklearnで書籍やネットに掲載されているコードを実行した結果、表題のようなエラーが出ることがある。

　一例をあげる。

>>> from sklearn.cross_validation import cross_val_score
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: No module named 'sklearn.cross_validation'

原因

　sklearnのバージョンの不一致が原因。

　sklearn 0.20からcross_validationモジュールは消滅し、代わりにmodel_selectionモジュールを使うことになった。
　（model_selectionモジュール自体は0.19とかでも存在する）

　なお、sklearnのバージョンは次のようなコマンドで確認できる。ただし、環境によっては微妙にコマンドが違う可能性もある（pip3だったりcondaだったり。常識的なことだが、念のために書いておく）。

$ pip show scikit-learn
Name: scikit-learn
Version: 0.20.0  # バージョン
Summary: A set of python modules for machine learning and data mining
Home-page: http://scikit-learn.org
Author: None
Author-email: None
License: new BSD
Location: ***  # インストール先
Requires: numpy, scipy
Required-by: sklearn

対処

　基本的にはcross_validationモジュールに存在した関数・クラス等はmodel_selectionモジュールに移動している。なので、以下のように書き換えれば正常に動作する。

# 古い書き方
from sklearn.cross_validation import cross_val_score

# 新しい書き方
from sklearn.model_selection import cross_val_score

　ただし、一部で使い方の変わっているクラスなどがあるので注意。KFoldクラスが例として挙げられる。変更点はドキュメントを見て確認する必要がある。

class sklearn.cross_validation.KFold(n, n_folds=3, shuffle=False, random_state=None)

sklearn.cross_validation.KFold — scikit-learn 0.19.2 documentation

class sklearn.model_selection.KFold(n_splits=’warn’, shuffle=False, random_state=None)

sklearn.model_selection.KFold — scikit-learn 0.20.1 documentation

　ところで実際問題として、cross_validationモジュールが将来的に消滅することは数年前からの規定事項だったので、cross_validationモジュールを使っている書籍やサイトは情報が古いと言わざるを得ない。なので、場合によっては、その情報源は切り捨てる（使いものにならないとみなす）ことも考慮する。

　どうしてもコードを書き換えないで使いたければ、古いバージョンのsklearnにダウングレードするという手はある。おすすめはできないが、一つの方法ではある。

$ pip uninstall scikit-learn
$ pip install scikit-learn==0.19

まとめ

　機械学習が流行り始めてから5年くらいは経っているので、そろそろ古い情報に惑わされる人も増えてきていると思う。

　どんなことにも言えるけど、情報の鮮度を意識することが重要。

プログラミングのブログにアドセンスを貼る話

はじめに

　一ヶ月ほど前から当ブログはGoogle Adsenseを導入しています*1。

　このブログはいわゆる「技術ブログ」と呼ばれるようなプログラミングのブログですが、プログラミング関連のジャンルでアドセンスを貼った場合の収益性についてはネット上にもあまり情報がありませんでした。なので貼る前は「びっくりするほど儲からなかったらどうしよう」などと思っていたのですが、とりあえずそれに関しては杞憂かな、という水準では稼げそうです。

　ご報告を兼ねて、内実をレポートします。

　なお、規約の絡みでPVや収益などの指標の具体的な数字は公表できません。あしからず。以下に出てくる数字は、有効数字一桁くらいのアバウトなものだと思ってください。基本的にボカして書いています。

このブログのPVと収益性

　まずこのブログの一日のPVですが、先月の数字を平均すると一日2000 PV弱くらいです。記事数はこの記事で183記事目で、1記事あたり10 PV強が得られていることになります。

　アドセンスの収益性の指標としては、ページRPMというものがあります。これは1000 PVで得られる収益のことです。これを1000で割れば収益[円]/閲覧数[PV]になるので、こちらで示します。当ブログでは0.2円/PV前後です。

　なので、アバウトですが一日400円弱ほど儲かるかなぁ、という収益性になっています。一ヶ月では一万円くらいです。これをしょぼいと見るか、小遣い稼ぎとしては優秀と見るかは、人によって分かれると思います。

収益の傾向

　プログラミングのブログというジャンルには、幾つかの特色があるようです。世間で言われているアドセンスとは、ちょっと違う傾向が見えました。

ユーザのほとんどはPCからアクセス

　アクセスの9割はPCからです。プログラミングをしながら調べごとをして当ブログにたどり着く方が多いと思うので、当たり前だと思います。

　このおかげで、PC向けに最適化すればとりあえず問題なく収益があがります。

　SEOとかアフィリエイトのサイトには、「現在はモバイルからのアクセスがほとんどなので、サイドバー広告は効果が少ない」なんて書いてあるところが幾らでもありますが、当ブログではサイドバー広告は立派な主力です。

　モバイル向けに最適化しようとすると、入れるのが面倒な上に本格的に邪魔な記事中広告を主力にせざるを得ないと思うので、これは嬉しい傾向です。サイトのユーザビリティを大きく損なわずにアドセンスで稼げます。
（余談ですが、一番稼いでくれるのはブログの一番上にでかでかと貼ってあるヘッダ広告です。印象は最悪なんだろうけど、全収益の2/5くらいをはじき出しています）

低いクリック率とそこそこのクリック単価

　よく巷のサイトには「ITスキルの高いユーザが多いサイトは広告のクリック率が下がる」などと書いてありますが、その理論で行くとプログラミングのジャンルのユーザなんてITスキル最強みたいなものですからクリック率は低いはずです。このブログも案の定というか、0.1%という「低い」クリック率をはじき出しています（平均で0.3%、本格的なアフィサイトであれば1%くらいは行くそうです）。

　それでも、クリックが「0ではない」というのは励みになる要素で、まあ何千人も来ていれば一定の確率で広告をクリックする人はいるということでしょう。

　「そんなにクリック率が低くて儲かるのかよ」と思う方もいると思いますが、広告単価が高いものが多いっぽく、なんとか補えます。恐らく、単価が高めの転職系、学習系、IT製品系などが表示されていることの結果だと思います。

　アドセンスのクリック単価の平均は30～40円くらいらしいですが、当ブログはその倍は軽く叩き出しています。

　なので、最終的にはクリック率の低さと広告単価の高さが相殺しあって、世間並みにはカネになってるのかな？　という印象です。

準主力のCPM

　アドセンスというと、普通はクリック単価広告の印象が強いと思いますが、実はCPMというインプレッション単価型の広告（広告の表示で収益が発生する/広告主からすると表示される回数にお金を払う）も存在します。

　なんか、軽く調べた感じだと、普通のサイトではCPM型の総収益比は2割未満らしいのですが……何しろ上述の通りのクリック率なので、このブログではCPMも準主力くらいには稼いでくれます。日によって変動はありますが、一日の収益の1/3～1/2くらいは稼いでくれていると思います。インプレッション収益も、クリック単価型と遜色ありません。むしろ日によっては上回ることすらあります。

　CPM型の広告は、PVに比例して着実に収益が伸びていきます。サイト運営者からすれば運次第な要素のあるクリック単価型広告とは違い、十数分おきに確認するたびに（確認してどうするんだ、って話ですが）コツコツ1円ずつ収益が積み上がっていくので、心に優しい広告です。

割に合うの？

　これについて知りたい方も多いと思いますが、結論から言えばたぶん仕事としては割に合わないです。ある程度効率的に運営できたとしても、時給500円くらいだと思います（記事を書く労力、維持管理する労力を適当に時間換算すれば、自ずと結論はでます）。

　ただし、ブログ運営は「趣味」ですから、そう考えると「趣味」をやるだけで、

• 時給500円稼げる
• さやかながら我が国のITの発展に貢献できる（かもしれない）

　なので、皆さんも（書くネタがあれば）プログラミングとか技術系のブログを作ってみてください*2。

結論

　儲かるかというと微妙なのですが、稼げないというほどでもないという程度には稼げます。

　強弁すれば小遣い稼ぎにはちょうど良いと見れなくもありませんが、積極的に人におすすめはしません。でも、このブログを完全放置しても月1万円になるので*3、私は貼り続けます。そんなところです。

*1:邪魔と思われている読者の方も当然おられるかと思いますが、私の小遣いに直結する問題なので、生暖かい目で見守っていただければなぁ、と思います。

*2:でもものすごく流行っちゃって、広告の需給が悪化して単価下がったりすると（私が）困るので、この記事を読んで「プログラミングのブログを始めてみようかな」と思い立つ人は2人か3人いればそれで良いです

*3:ただし1万円をずっと維持できる訳ではないことに注意。プログラミングの記事なんて、数年もすれば検索順位下落・トレンドの移り変わりで読まれなくなることでしょう

【python】sklearnでのカテゴリデータの取り扱いまとめ

概要

　カテゴリデータをone-hot表現として取り扱うという方法は、機械学習などでは一般的に行われます。

　しかし、この辺りの処理に対するsklearnでのサポートが微妙に悪いという問題が長年あり、やれpandasを使えだの、やれサードパーティ製ライブラリで凌げだのといった話題が乱立していました。

　sklearn 0.20からはOneHotEncoderが拡張され、sklearnの枠組みの中で簡単にカテゴリデータ、カテゴリ変数をone-hot特徴量に変換できるようになりました。

　目次

カテゴリデータを取り扱うために使用するクラス

　以下で使うクラスについて説明します。

LabelEncoder

　これはラベルの変換に使います。教師ラベルを変換するのに使えます。教師ラベルはあえて変換しなくても、そのまま受け付けてくれるモデルがsklearnでは多い気がするので、なくても良いような気はします。

>>> from sklearn.preprocessing import LabelEncoder
>>> le = LabelEncoder()
>>> le.fit_transform(["hoge", "fuga", "piyo"])
array([1, 0, 2])

　インデックスの順番は、たぶんソートして昇順でしょう。

sklearn.preprocessing.LabelEncoder — scikit-learn 0.20.1 documentation

OrdinalEncoder

　LabelEncoderに似ていますが、ラベルではなくデータを受け付けるもの、と認識しておけば良いと思います。要するに、慣れ親しんだ[n_samples, n_features]の配列を受け取ります。

>>> from sklearn.preprocessing import OrdinalEncoder
>>> oe = OrdinalEncoder()
>>> oe.fit_transform([["h", "ho", "hoge"], ["f", "fu", "fuga"], ["p", "pi", "piyo"]])
array([[1., 1., 1.],
       [0., 0., 0.],
       [2., 2., 2.]])

sklearn.preprocessing.OrdinalEncoder — scikit-learn 0.20.1 documentation

　これで良いときはこれを使えば良いのですが、one-hot表現がほしいときが多いと思います。その場合は、次節で示すOneHotEncoderが使えます。

OneHotEncoder

　最初に書いた通り、sklearn 0.20からOneHotEncoderがOrdinalEncoder相応のデータを受け付けるようになりました（それ以前はそういう仕様ではなく、LabelEncoderかOrdinalEncoderを通して数値化したデータを入れる必要がありました）。

　なので、今後は事実上これだけ使えば良いということです。

>>> from sklearn.preprocessing import OneHotEncoder
>>> ohe = OneHotEncoder()
>>> ohe.fit_transform([["h", "ho", "hoge"], ["f", "fu", "fuga"], ["p", "pi", "piyo"]])
<3x9 sparse matrix of type '<class 'numpy.float64'>'
	with 9 stored elements in Compressed Sparse Row format>
>>> ohe = OneHotEncoder(sparse=False)
>>> ohe.fit_transform([["h", "ho", "hoge"], ["f", "fu", "fuga"], ["p", "pi", "piyo"]])
array([[0., 1., 0., 0., 1., 0., 0., 1., 0.],
       [1., 0., 0., 1., 0., 0., 1., 0., 0.],
       [0., 0., 1., 0., 0., 1., 0., 0., 1.]])

　注意点として、上でも示していますがデフォルトパラメータでは返り値がsparse matrixです。sparse=False一発でnumpy配列を返すようになります。sklearnのモデルの疎行列対応はモデルによって対応していたり、していなかったりという具合なので、後段のモデルに色々なものを使いたければそちらを指定することをおすすめします。

sklearn.preprocessing.OneHotEncoder — scikit-learn 0.20.1 documentation

まとめ

　OneHotEncoderを使えばシンプルにやりたいことが実現できます。