CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: Twine を使って PyPI にパッケージをアップロードする

Python のサードパーティ製パッケージは、一般に PyPI (Python Package Index) というリポジトリに登録されている。 このリポジトリはユーザ登録さえすれば誰でも自分で作ったパッケージをアップロードできる。 今回は Twine というツールを使って PyPI にパッケージをアップロードしてみる。

PyPI - the Python Package Index : Python Package Index

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.6
BuildVersion:   16G1114
$ python --version
Python 3.6.3
$ twine --version
twine version 1.9.1 (pkginfo: 1.4.1, requests: 2.18.4, setuptools: 38.2.5,
requests-toolbelt: 0.8.0, tqdm: 4.19.5)

どうして Twine を使うのか

以前は PyPI にパッケージをアップロードするときにはセットアップスクリプト (setup.py) の upload コマンドを使うのが一般的だった。 ただ、このやり方だと古い Python ではプロトコルに HTTP が使われるためユーザ名やパスワードが平文でネットワークに流れてしまう。 また、プロトコルに HTTP over SSL/TLS を使っていても、証明書を検証していないと中間者攻撃のリスクがある。

Issue 12226: use HTTPS by default for uploading packages to pypi - Python tracker

そのため、プロトコルが HTTP over SSL/TLS でちゃんと証明書の検証までやってくれる Twine を使うことがセキュリティ的に望ましい。 また、ファイルをアップロードするときの使い勝手という点でも Twine の方が優れている点がある。

ちなみに Twine は PyPA (Python Packaging Authority) というワーキンググループがメンテナンスしている。 これは setuptoolspipvirtualenv といったパッケージ関連でデファクトなツールのメンテナンスをしているコミュニティ。

Python Packaging Authority — PyPA documentation

インストール

前置きが長くなったけど、ここからは Twine を実際に使っていく。

まずは Twine を pip を使ってインストールする。

$ pip install twine

PyPI にユーザ登録する

続いて PyPI にユーザ登録をする。 以下の Web サイトの「Register」からユーザ登録を完了しよう。

PyPI - the Python Package Index : Python Package Index

あるいは、新しい Web サイトの以下からでも良い。

PyPI - the Python Package Index · Warehouse (PyPI)

同時にテスト用の PyPI にも別途ユーザ登録しておこう。 こちらにも登録しておくとリリース作業の予行演習ができる。

PyPI - the Python Package Index : Python Package Index

テスト用の PyPI についても新しい Web サイトがある。

PyPI - the Python Package Index · Warehouse (TestPyPI)

設定ファイルを用意する (オプション)

必須ではないんだけど、設定ファイルを用意しておくとこれからの作業が捗る。 ホームディレクトリに .pypirc という名前でファイルを作っておくと Twine がデフォルトで読み込んでくれる。

書式はこんな感じ。 <your_username><your_password> の部分は、先ほど自分で登録したユーザ名とパスワードに置き換える。

$ cat << 'EOF' > ~/.pypirc 
[distutils]
index-servers =
  pypi
  pypitest

[pypi]
repository=https://upload.pypi.org/legacy/
username=<your_username>
password=<your_password>

[pypitest]
repository=https://test.pypi.org/legacy/
username=<your_username>
password=<your_password>
EOF

アップロードしたいパッケージを用意する

次は PyPI にアップロードしたいパッケージを用意する。 今回は自作パッケージの diagram-autobuild を例にして説明する。

github.com

パッケージをビルドする

用意ができたら、まずはアップロードしたいパッケージの成果物をビルドする。

例えばソースコード配布物をビルドするなら、次のようにセットアップスクリプトに sdist を指定する。

$ python setup.py sdist

Wheel パッケージなら bdist_wheel を指定する。

$ python setup.py bdist_wheel

するとプロジェクトの dist ディレクトリ以下に成果物ができる。

$ ls dist
diagram-autobuild-0.1.0.tar.gz               diagram_autobuild-0.1.0-py2.py3-none-any.whl

パッケージをアップロードする

あとは twine upload コマンドを使って dist ディレクトリ以下のアップロードしたいファイルを指定するだけ。 アップロード先は --repository オプションで指定する。 まずはテスト用の PyPI にアップロードして内容を確認しよう。

$ twine upload \
  --repository pypitest \
  dist/*

ちなみに、先ほどの設定ファイルを作っていない状態でテスト用の PyPI にアップロードしたいときは --repository-url を使えば良い。

$ twine upload \
  --repository-url https://test.pypi.org/legacy/ \
  dist/*

テスト用の PyPI に意図通りパッケージがアップロードできていることを確認しておこう。

diagram-autobuild · Warehouse (TestPyPI)

問題がなければ本番の PyPI にもリリースする。

$ twine upload \
  --repository pypi \
  dist/*

本番の PyPI にパッケージがアップロードできていることを確認しよう。

diagram-autobuild · Warehouse (PyPI)

最後に、ちゃんとアップロードしたパッケージを pip 経由でインストールできることを確認できたら作業はおわり。

$ pip install diagram-autobuild
$ pip list --format=columns | grep diagram-autobuild
diagram-autobuild 0.1.0  

めでたしめでたし。

Python: ERAlchemy を使って ER 図を描く

今回は ERAlchemy という ER 図を描くツールを使ってみる。 このツールは erd という Haskell で書かれた同様のツールにインスパイアされて作られたものらしい。 ただ、機能的にできることは ERAlchemy の方が多いみたいだ。

ERAlchemy が提供する基本的な機能は次の通り。

  • ER フォーマットのテキストファイルから ER 図を生成する
  • SQLAlchemy 経由で既存のデータベースから ER 図を生成する

後者の既存データベースから ER 図を生成するところなんかは、これまでだと MySQL Workbench を使ったりしてた。 ただ、このやり方だと文字通り MySQL でしか使えないのに対して ERAlchemy はそれ以外のデータベースにも対応している。 今回も試しに SQLite3 のデータベースから ER 図を生成してみている。 ただ、この機能が出力する図は、ちょっと直感には反する図になってしまった。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.12.6
BuildVersion:   16G1114
$ python --version
Python 3.6.3

インストール

まずは Homebrew を使って Graphviz を Pango と一緒にインストールしておく。

$ brew reinstall graphviz --with-pango

続いて ERAlchemy 本体をインストールする。

$ brew install eralchemy

ちなみに ERAlchemy 自体は Python で書かれているので pip 経由でインストールしても構わない。

$ pip install eralchemy

ER ファイルから ER 図を生成する

まず、基本となる使い方として ER フォーマットというテキストファイルから ER 図を生成してみる。 ER フォーマットの詳細については erd で詳しく説明されている。 とはいえ、ここでもその概要については説明する。

例えば、次の ER ファイルでは users というテーブルの中に三つのカラムを定義している。 主キーとなる id と、名前と年齢を入れる nameage だ。

$ cat << 'EOF' > example.er 
[users]
*id
name
age
EOF

上記を ER 図にレンダリングしてみよう。 eralchemy-i オプションに上記のファイルを指定したら、出力先を -o で指定する。

$ eralchemy -i example.er -o example.png

するとレンダリングされた画像が次のように生成される。

$ file example.png 
example.png: PNG image data, 89 x 159, 8-bit/color RGBA, non-interlaced

f:id:momijiame:20171230060152p:plain

レンダリング先のフォーマットとしては画像ファイル以外に PDF にも対応している。

$ eralchemy -i example.er -o example.pdf
$ file example.pdf 
example.pdf: PDF document, version 1.3

リレーションを表現する

リレーショナルデータベースの ER 図を描くのだから当然リレーションについても図示できないとまずい。 続いてはリレーションを含む ER 図を描いてみよう。

次の例ではテーブル users とテーブル emails が 1:n 対応しているところを表現している。 テーブルを定義してから、その後ろでそれらの関係性を書いていく感じ。 ちなみに本家の erd では外部キー参照を + で表現するようだけど ERAlchemy ではまだサポートしていないようだ。

$ cat << 'EOF' > example.er 
[users]
*id
name
age

[emails]
*id
address
user_id

users 1--* emails
EOF

上記をレンダリングしてみよう。

$ eralchemy -i example.er -o example.png

すると、次のような ER 図が得られる。

f:id:momijiame:20171230060525p:plain

先の ER ファイルでテーブル間のリレーションを表現するのには -- の前後に対応関係を表す記号を入れていた。 この中で登場するのは 1* だけだったけど、それ以外には ?+ も使うことができる。

意味 記号
0 or 1 ?
exactly 1 1
0 or more *
1 or more +

型情報などのラベルをつける

リレーショナルデータベースの ER 図を描くのであれば、当然それぞれのカラムの型についても情報がほしい。 そういったときは次のように label で情報を付与する。

$ cat << 'EOF' > example.er 
[users]
*id {label: "INTEGER"}
name {label: "TEXT"}
age {label: "INTEGER"}

[emails]
*id {label: "INTEGER"}
address {label: "TEXT"}
user_id {label: "INTEGER"}

users 1--* emails
EOF

同様にレンダリングする。

$ eralchemy -i example.er -o example.png

上記から得られた ER 図は次の通り。

f:id:momijiame:20171230061628p:plain

既存のデータベースから SQLAlchemy 経由で ER 図を描く

ERAlchemy の特徴として、SQLAlchemy 経由で既存のデータベースから ER 図を描く機能が挙げられる。 試しに SQLite3 のデータベースを作って、そこから ER 図を書いてみることにしよう。

まずは、次のように SQLite3 のデータベースを用意する。

$ sqlite3 example.db
SQLite version 3.16.0 2016-11-04 19:09:39
Enter ".help" for usage hints.
sqlite> CREATE TABLE users (
   ...> id INTEGER NOT NULL, 
   ...> name TEXT NOT NULL, 
   ...> PRIMARY KEY (id)
   ...> );
sqlite> CREATE TABLE emails (
   ...> id INTEGER NOT NULL, 
   ...> address TEXT NOT NULL, 
   ...> user_id INTEGER, 
   ...> PRIMARY KEY (id), 
   ...> FOREIGN KEY(user_id) REFERENCES users (id)
   ...> );
sqlite> .exit
$ file example.db 
example.db: SQLite 3.x database, last written using SQLite version 3016000

用意ができたら ERAlchemy を使って ER 図をレンダリングする。 今度は -i オプションに SQLAlchemy でデータベースの接続に使う URI 形式を渡すのがポイントとなる。

$ eralchemy -i sqlite:///example.db -o example.png

レンダリングされた図は以下の通り。 users の方が 0...N になっていて、なんだから emails の方が主体のような図になってしまった。 感覚的には反対になる気がするんだけどなあ。

f:id:momijiame:20171230062450p:plain

まとめ

今回は ERAlchemy を使って ER 図を描いてみた。 既存のデータベースから ER 図を作る機能は、ちょっと「ん?」という結果になってしまった。 とはいえ ER ファイルから図をレンダリングする部分に関してはちゃんと動くみたいなので、上手く活用していきたい。

macOS に Homebrew で PostgreSQL をインストールする

とあるライブラリのために PostgreSQL を使う機会があったので、その作業メモ。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.6
BuildVersion:   16G1114
$ brew --version
Homebrew 1.4.1
Homebrew/homebrew-core (git revision 86b0; last commit 2017-12-29)

Homebrew を使って PostgreSQL をインストールする。

$ brew install postgresql

バージョンは 10.1 が入った。

$ psql --version
psql (PostgreSQL) 10.1

インストールすると brew servicespostgresql が登録される。

$ brew services list
Name       Status  User Plist
postgresql stopped      

start コマンドで PostgreSQL のサービスを起動する。

$ brew services start postgresql
==> Successfully started `postgresql` (label: homebrew.mxcl.postgresql)

インストール直後で作られているデータベースはこんな感じ。

$ psql -l
                                List of databases
   Name    |  Owner  | Encoding |   Collate   |    Ctype    |  Access privileges  
-----------+---------+----------+-------------+-------------+---------------------
 postgres  | amedama | UTF8     | en_US.UTF-8 | en_US.UTF-8 | 
 template0 | amedama | UTF8     | en_US.UTF-8 | en_US.UTF-8 | =c/amedama         +
           |         |          |             |             | amedama=CTc/amedama
 template1 | amedama | UTF8     | en_US.UTF-8 | en_US.UTF-8 | =c/amedama         +
           |         |          |             |             | amedama=CTc/amedama
(3 rows)

psql コマンドを使って接続できる。

$ psql -d postgres
psql (10.1)
Type "help" for help.

postgres=#

めでたしめでたし。

Python: 機械学習で分類問題のモデルを評価する指標について

今回は、機械学習において分類問題のモデルを評価するときに使われる色々な指標について扱う。 一般的な評価指標としては正確度 (Accuracy) が使われることが多いけど、これには問題も多い。 また、それぞれの指標は特徴が異なることから、対象とする問題ごとに重視するものを使い分ける必要がある。

今回扱う代表的な評価指標は次の通り。

  • 正確度 (正解率、Accuracy)
  • 適合率 (精度、陽性反応的中度、Precision)
  • 再現率 (感度、真陽性率、Recall)
  • F-値 (F-score, F-measure)
  • AUC (Area Under the Curve)

上記それぞれの指標について、特徴を解説すると共に Python を使って計算してみる。 データセットには scikit-learn に組み込みの乳がんデータセットを用いた。 今回は「機械学習で」と書いてしまったけど、上記は実際には医療統計の分野でもよく使われる指標だったりする。

尚、今回使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.6
BuildVersion:   16G1114
$ python --version
Python 3.6.3

もくじ

下準備

まずは Python で計算するために必要な準備をする。 しばらく地味な作業が続くため、興味がない人は次のセクションまで飛ばしてもらっても良いかも。

最初に、必要なパッケージをインストールしておこう。

$ pip install numpy scipy scikit-learn matplotlib

続いて Python の REPL を起動する。

$ python

そして、乳がんデータセットをロードしておく。

>>> from sklearn import datasets
>>> breast_cancer = datasets.load_breast_cancer()

この乳がんデータセットには、腫瘍について 30 項目のファクターが特徴量として収められている。

>>> breast_cancer.feature_names
array(['mean radius', 'mean texture', 'mean perimeter', 'mean area',
       'mean smoothness', 'mean compactness', 'mean concavity',
       'mean concave points', 'mean symmetry', 'mean fractal dimension',
       'radius error', 'texture error', 'perimeter error', 'area error',
       'smoothness error', 'compactness error', 'concavity error',
       'concave points error', 'symmetry error', 'fractal dimension error',
       'worst radius', 'worst texture', 'worst perimeter', 'worst area',
       'worst smoothness', 'worst compactness', 'worst concavity',
       'worst concave points', 'worst symmetry', 'worst fractal dimension'],
      dtype='<U23')
>>> len(breast_cancer.feature_names)
30

そして、その腫瘍が悪性 (malignant) か良性 (benign) かラベル付けされている。

>>> breast_cancer.target_names
array(['malignant', 'benign'],
      dtype='<U9')

データセットの説明が済んだところで、特徴量とラベルを取り出しておく。

>>> X, y =  breast_cancer.data, breast_cancer.target

続いては教師データを学習用とテスト用に分割しておく。 学習用のデータで評価指標を測ってしまうと、それは汎化性能を表さないため。 汎化性能というのは、未知のデータに対するモデルの対処能力のことを指す。

>>> from sklearn.model_selection import train_test_split
>>> X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5)

ここらへんは交差検証や交差検定、とかで調べると色々分かる。

今回は分類器のモデルとしてランダムフォレストを用いた。 特にチューニングなどはせず、そのまま使う。

>>> from sklearn.ensemble import RandomForestClassifier
>>> clf = RandomForestClassifier()

ランダムフォレストを学習用のデータを使って学習する。

>>> clf.fit(X_train, y_train)
RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=None, max_features='auto', max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
            oob_score=False, random_state=None, verbose=0,
            warm_start=False)

続いてテスト用データを使って推定する。

>>> y_pred = clf.predict(X_test)

ここまでで評価指標を紹介するための準備が整った。

評価指標を計算する

続いて、先ほどのセクションで得られた y_predy_test を使って評価指標を計算していく。 補足しておくと y_pred はモデルが推定した値で y_test は真の値になっている。

まずは、評価指標を説明する上で必要な定義について確認しておく。 説明の中では、モデルの判定結果を 4 種類に分けて扱う。 具体的には、次のようなマトリックスで表すことができる。 このマトリックスは混同行列 (Confusion Matrix) と呼ばれる。 横軸は事実を、縦軸はモデルや検査の推定を表している。

悪性 良性
陽性 TP FP
陰性 FN TN

セルに入った英字は、それぞれ以下のような意味がある。

  • TP: True Positive
    • 正しく陽性と判定できた場合
  • TN: True Negative
    • 正しく陰性と判定できた場合
  • FP: False Positive
    • 本来は陰性なところを、誤って陽性と判定してしまった場合
  • FN: False Negative
    • 本来は陽性なところを、誤って陰性と判定してしまった場合

scikit-learn で混同行列を計算するときは confusion_matrix() 関数が使える。

>>> from sklearn.metrics import confusion_matrix
>>> confusion_matrix(y_test, y_pred)
array([[118,   1],
       [ 45, 121]])

正確度 (正解率、Accuracy)

まずは評価指標として最も一般的な正確度 (Accuracy) から紹介していく。 これは、ようするに推定した値と真の値が一致した割合になる。

正確度が使われる場面では FPFN の重要度について特に考慮しなくても良い場合が考えられる。 ちなみに、今回のケースでは FPFN の重要度は、あきらかに異なる。 なぜなら、本当は悪性の患者を良性と判断してしまう (FN) のは患者にとって文字通り致命的なため。 このような場面では、正確度だけを見てモデルの良し悪しを判断することはできない。

また、正確度だけ見ているとデータセットの性質によってはモデルを正しく評価するのが難しいケースもある。 これは、具体的には陽性と陰性の出現する割合が極端に異なる場合が挙げられる。 例えば、医療統計の分野では実際に病気の人は健康な人に比べるとずっと少ないことがよくある。 そのような場合、数の多い健康な人をとりあえず健康だと判断しているだけでも正確度は高く出てしまう。 これではモデルの良し悪しを正しく評価することは難しい。

正確度の定義は次の通り。

 Accuracy = \frac{TP + TN}{TP + FP + FN + TN}

定義通りに計算しても良いけど scikit-learn に専用の関数が用意されているので使う。

>>> from sklearn.metrics import accuracy_score
>>> accuracy_score(y_test, y_pred)
0.93333333333333335

このように、評価指標としてよく使われる正確度だけど、実際の運用では適用が難しいケースも多いことが分かった。

適合率 (精度、陽性反応的中度、Precision)

続いて紹介する評価指標は適合率 (Precision) というもの。 これは、モデルが陽性と判断したものの中に、どれだけ本当に陽性なものが含まれていたかを示す指標になる。 いわば正確性を見るための指標といえる。 今回のケースであれば、モデルが悪性と判断したものの中にどれだけ本当に悪性のものが含まれていたか。 後述する再現率とはトレードオフの関係にある。

この評価指標は、あきらかに陽性と分かりやすいものだけを見つけたいときに重視するものと考えられる。 この評価指標を重視すると、分かりにくいものは基本的に陰性と判断することになる。 当然ながら、その中には実際には陽性だったものがたくさん含まれているはず。 つまり、適合率を重視するときは FN が発生することが許容できるケースといえる。 それでもなお FP があっては困るときに適合率を使うことができる。

今回のケースでいうと、この評価指標を重視すべきでないことはあきらかだ。 FN が発生した場合、それは患者の生命に関わる。

定義は次の通り。

 Precision = \frac{TP}{TP + FP}

定義通りに計算しても良いけど、先ほどと同じように専用の関数を使ってしまうのが楽ちん。

>>> from sklearn.metrics import precision_score
>>> precision_score(y_test, y_pred)
0.95402298850574707

再現率 (感度、真陽性率、Recall)

続いて紹介するのは再現率 (Recall) という評価指標になる。 これは、実際に陽性だったもののうちモデルが陽性と判断したものの割合を指す。 いわば網羅性を見るための指標といえる。 今回のケースであれば本当に悪性だったものの中でモデルが悪性と判断できたものの割合になる。 先述した適合率とはトレードオフの関係にある。

この評価指標は、怪しいものはとりあえず全て見つけ出したいときに重視すべきものと考えられる。 この評価指標を重視すると、分かりにくいものは基本的に陽性と判断することになる。 当然ながら、その中には実際には陰性だったものがたくさん含まれているはず。 つまり、再現率を重視するときは FP が発生することが許容できるケースといえる。 それでもなお FN があっては困るときに再現率を使うことができる。

今回のケースでいうと、この評価指標を重視するのが適当とかんがえられる。 FP は、再検査などを受けて実際には陰性と分かったとき「よかったですね」で済む。 もちろん無いに越したことはないものの FN が発生することに比べれば、ずっと重要度は低い。

定義は次の通り。

 Recall = \frac{TP}{TP + FN}

これについても scikit-learn に計算用の関数が用意されている。

>>> from sklearn.metrics import recall_score
>>> recall_score(y_test, y_pred)
0.93785310734463279

適合率と再現率のトレードオフについて

前述した通り、適合率と再現率はトレードオフの関係にある。 ここでは、それについて単純化した例を使って説明したい。

まずは次の図を見てほしい。 これは、ある因子が取る値によって疾病の有無の確率が変化することを示したグラフになっている。 例えば血液検査して得られた特定の項目が高いと、ある疾病を持っている確率が高い、みたいな感じ。 もちろん、これは説明するために用意したものなので実際にこのような関係になっているものが具体的にあるわけではない。

f:id:momijiame:20171218003752p:plain

ここで、要因は一つしかないので横軸の何処かに識別境界を置いて疾病の有無を判断することになる。 ようするに、この閾値より上であれば疾病ありと判断、下であれば疾病なしと判断するというわけ。

f:id:momijiame:20171218003809p:plain

例えば識別境界をど真ん中に置いてみたときはどうなるだろう。 誤って疾病あり (FP) と判断される人と誤って疾病なし (FN) と判断される人が同じように出ることになる。

f:id:momijiame:20171218003819p:plain

続いては、識別境界を右にずらしてみよう。 このときは、誤って疾病あり (FP) と判断される人は減るため適合率は上がることになる。 しかし、反対に誤って疾病なし (FN) と判断される人は増えることから再現率は下がってしまう。

f:id:momijiame:20171218003837p:plain

同じように、識別境界を今度は左にずらしてみよう。 すると、誤って疾病なし (FN) と判断される人は減るため再現率は上がる。 反対に、誤って疾病あり (FP) と判断される人は増えるので適合率は下がってしまうことが分かる。

f:id:momijiame:20171218003845p:plain

このように、単純化した例を使うと両者がトレードオフの関係にあることが分かりやすい。

適合率と再現率の調整について

続いては、実際に scikit-learn を使って適合率と再現率を調整してみよう。

元々の predict() メソッドでは閾値を 0.5 にしてある。

>>> clf.predict(X_test)
array([1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0,
       0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0,
       1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1,
       1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1,
       1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0,
       1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0,
       1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0,
       1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1,
       0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1,
       1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0,
       1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1,
       1, 0, 1, 1, 1, 1, 0, 1, 1])

一体何が 0.5 なのかというと、そのクラスに分類される閾値のこと。 predict_proba() メソッドを使うと、そのクラスに分類される確率が得られる。 例えば 4 番目の項目はクラス 0 に分類される確率が 0.6 でクラス 1 に分類される確率が 0.4 だった。 閾値が 0.5 にあるとき、この項目はクラス 0 に分類されることになる。

>>> clf.predict_proba(X_test)
array([[ 0. ,  1. ],
       [ 0.1,  0.9],
       [ 1. ,  0. ],
       [ 0.6,  0.4],
...(snip)...
       [ 0.9,  0.1],
       [ 0. ,  1. ],
       [ 0. ,  1. ]])

次のようにすると、自分で閾値を変更した分類結果が得られる。

>>> (clf.predict_proba(X_test)[:, 0] < 0.5).astype(int)
array([1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0,
       0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0,
       1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1,
       1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1,
       1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0,
       1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0,
       1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0,
       1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1,
       0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1,
       1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0,
       1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1,
       1, 0, 1, 1, 1, 1, 0, 1, 1])

実際に閾値をずらしながら評価指標が変化することを確認してみよう。 まずは閾値を変えることで再現率を上げてみる。

>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.5).astype(int)
>>> recall_score(y_test, y_pred)
0.96385542168674698
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.6).astype(int)
>>> recall_score(y_test, y_pred)
0.98795180722891562
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.7).astype(int)
>>> recall_score(y_test, y_pred)
0.99397590361445787
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.8).astype(int)
>>> recall_score(y_test, y_pred)
1.0

閾値を反対に変化させれば適合率を上げることができる。

>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.5).astype(int)
>>> precision_score(y_test, y_pred)
0.97560975609756095
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.4).astype(int)
>>> precision_score(y_test, y_pred)
0.97530864197530864
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.3).astype(int)
>>> precision_score(y_test, y_pred)
0.98709677419354835
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.2).astype(int)
>>> precision_score(y_test, y_pred)
0.99305555555555558
>>> y_pred = (clf.predict_proba(X_test)[:, 0] < 0.1).astype(int)
>>> precision_score(y_test, y_pred)
0.99180327868852458

また、次のように precision_recall_curve() 関数を使うと両者の関係をグラフ化できる。

>>> from sklearn.metrics import precision_recall_curve
>>> precision, recall, thresholds = precision_recall_curve(y_test, clf.predict_proba(X_test)[:, 1])
>>> plt.step(recall, precision, color='b', alpha=0.2,
...          where='post')
[<matplotlib.lines.Line2D object at 0x1189cac18>]
>>> plt.fill_between(recall, precision, step='post', alpha=0.2,
...                  color='b')
<matplotlib.collections.PolyCollection object at 0x1189cadd8>
>>>
>>> plt.xlabel('Recall')
Text(0.5,0,'Recall')
>>> plt.ylabel('Precision')
Text(0,0.5,'Precision')
>>> plt.ylim([0.0, 1.05])
(0.0, 1.05)
>>> plt.xlim([0.0, 1.0])
(0.0, 1.0)
>>> plt.show()

今回のモデルであれば次のようなグラフが得られた。

f:id:momijiame:20171218003905p:plain

F-値 (F-score, F-measure)

続いて紹介する評価指標は、F-値 (F-score, F-measure) というもの。 これは、先述した適合率と再現率の調和平均を取ったものになっている。 両者のバランスが極端に悪くないものを作りたいときは、この評価指標を使うのが良いようだ。

正確度では陽性と陰性の出現する割合が極端に異なると、正しくモデルを評価するのが難しかった。 それに対しF-値はそのような場合であっても評価しやすい。 そのため、常に正確度ではなくF-値を使って評価することを推奨する人もいるようだ。

定義は次の通り。

 F-measure = 2 \frac{Precision \cdot Recall}{Precision + Recall}

これも scikit-learn に計算用の関数がある。

>>> from sklearn.metrics import f1_score
>>> f1_score(y_test, y_pred)
0.94586894586894577

AUC (Area Under the Curve)

続いて紹介するのは AUC (Area under the curve) という評価指標になる。 ただ、この AUC を理解するには、その前に ROC (Receiver Operating Characteristic) 曲線というものを理解しなきゃいけない。 ROC 曲線は、日本語だと受信者動作特性曲線とも言ったりする。

ROC 曲線というのは、縦軸に再現率を、横軸に特異度をプロットしたグラフのことをいう。 特異度という新しい言葉が登場したけど、これは実際に陰性だったもののうちモデルが陰性と判断したものの割合を指す。 特異度は、真陰性率とも呼ばれるほか英語では Specificity となる。

特異度の定義は次の通り。

 Specificity = \frac{TN}{TN + FP}

上記、再現率と特異度の関係を二次元でグラフ化したものが ROC 曲線と呼ばれる。

実際に今回のモデルを使って ROC 曲線を描いてみよう。 scikit-learn で ROC 曲線を得るには roc_curve() 関数が使える。

>>> from sklearn.metrics import roc_curve
>>> fpr, tpr, _ = roc_curve(y_test, clf.predict_proba(X_test)[:, 1])
>>> plt.step(fpr, tpr, color='b', alpha=0.2, where='post')
[<matplotlib.lines.Line2D object at 0x131100cc0>]
>>> plt.fill_between(fpr, tpr, step='post', alpha=0.2, color='b')
<matplotlib.collections.PolyCollection object at 0x12b7b22b0>
>>> plt.plot([0, 1], [0, 1], color='navy', linestyle='--')
[<matplotlib.lines.Line2D object at 0x12b7b2160>]
>>> plt.ylabel('True Positive Rate')
Text(0,0.5,'True Positive Rate')
>>> plt.xlabel('False Positive Rate')
Text(0.5,0,'False Positive Rate')
>>> plt.ylim([0.0, 1.0])
(0.0, 1.0)
>>> plt.xlim([0.0, 1.0])
(0.0, 1.0)
>>> plt.show()

上記で得られたグラフが次の通り。

f:id:momijiame:20171218003922p:plain

ROC 曲線では、上記の青い部分が多いほど (白い部分が少ないほど) 優れたモデルということを表している。 そして、その優れたモデルかどうかを表す評価指標が実は AUC ということになる。 AUC というのは、ようするに青い部分の面積を 0.5 ~ 1 の間で表現したものだ。

どうして 0.5 以上になるかというと、それは ROC 曲線の特性とも関わってくる。 ROC 曲線に斜めの破線が描かれていたことに気づいただろうか。 実は、二値分類を完全にランダムに実行すると、ROC 曲線はこの斜めの破線上に存在して AUC も 0.5 になる。 もし、ROC 曲線が斜めの破線より下にあって AUC が 0.5 未満ということは、つまりランダムに分類するよりも性能が悪いことになる。 その場合だと、予測内容を反転してしまった方がマシということだ。 つまり AUC は 0.5 以上でないとおかしい。

scikit-learn で AUC を計算するには auc() 関数が使える。 先ほど計算した再現率と特異度を渡してやろう。

>>> from sklearn.metrics import auc
>>> auc(fpr, tpr)
0.99141945935000508

まとめ

今回は、機械学習において分類問題のモデルを評価するときに使われる色々な指標を紹介した。 評価指標はそれぞれ特性が異なるため、問題に応じて適材適所で使い分ける必要がある。

Python: 多様体学習 (Manifold Learning) を用いた次元縮約

今回は多様体学習を使ってデータの次元を縮約する方法について。 これはデータの前処理として、主に二つの目的で使われる。 一つ目は、次元を縮約することで二次元や三次元の形でデータを可視化できるようにするため。 もう一つは、次元を縮約した結果を教師データとして用いることでモデルの認識精度を上げられる場合があるため。

データの次元を縮約する手法としては主成分分析 (PCA) が有名だけど、これは線形な変換になっている。 ただ、実際に取り扱うデータは必ずしもそれぞれの次元が線形な関係になっているとは限らない。 そこで、非線形な変換をするのが多様体学習ということらしい。

今回使った環境は次の通り。

$ sw_vers      
ProductName:    Mac OS X
ProductVersion: 10.12.6
BuildVersion:   16G1114
$ python --version
Python 3.6.3

下準備

まずは、今回使う Python のライブラリをインストールしておく。

$ pip install numpy scipy scikit-learn matplotlib

扱うデータセットとしては scikit-learn の digits データセットにした。 これは 8 x 8 のピクセルで表現された手書きの数値データになっている。 8 x 8 ピクセルを扱うので 64 次元になっていて、それが 1797 点ある。

>>> from sklearn import datasets
>>> dataset = datasets.load_digits()
>>> dataset.data.shape
(1797, 64)

上記のデータセットが具体的にどういったものなのかを可視化しておく。 データセットからランダムに 25 点を取り出して画像として表示してみよう。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
from matplotlib import cm
import numpy as np
from numpy import random
from sklearn import datasets


def main():
    dataset = datasets.load_digits()
    X = dataset.data
    y = dataset.target

    # データの中から 25 点を無作為に選び出す
    sample_indexes = random.choice(np.arange(len(X)), 25, replace=False)

    # 選んだデータとラベルを matplotlib で表示する
    samples = np.array(list(zip(X, y)))[sample_indexes]
    for index, (data, label) in enumerate(samples):
        # 画像データを 5x5 の格子状に配置する
        plt.subplot(5, 5, index + 1)
        # 軸に関する表示はいらない
        plt.axis('off')
        # データを 8x8 のグレースケール画像として表示する
        plt.imshow(data.reshape(8, 8), cmap=cm.gray_r, interpolation='nearest')
        # 画像データのタイトルに正解ラベルを表示する
        plt.title(label, color='red')
    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を適当な名前をつけて保存したら実行する。

$ python digits.py

すると、次のようなグラフが得られる。 MNIST (28 x 28 ピクセル) に比べると、だいぶ荒いことが分かる。

f:id:momijiame:20171209032808p:plain

とはいえデータセットのダウンロードが生じないので取り回しが良い。

多様体学習を使ってデータセットの次元を縮約する

続いては、上記で確認したデータセットを実際に多様体学習アルゴリズムを使って次元縮約してみる。 元々のデータセットは 64 次元なので、個々を画像として表示はできるものの全体を散布図のように図示することはできない。 そこで、多様体学習を用いて 2 次元に縮約することで散布図として図示できるようにしてしまおう、ということ。

scikit-learn に組み込まれている多様体学習アルゴリズムの一覧は次のページで確認できる。

2.2. Manifold learning — scikit-learn 0.19.1 documentation

次のサンプルコードでは、上記からいくつか主要なアルゴリズムを使っている。 digits データセットの 64 次元を 2 次元に縮約した結果を散布図として図示した。 一応、比較対象として主成分分析も入れている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
from sklearn import datasets

from sklearn.decomposition import PCA
from sklearn.manifold import MDS
from sklearn.manifold import LocallyLinearEmbedding
from sklearn.manifold import SpectralEmbedding
from sklearn.manifold import Isomap
from sklearn.manifold import TSNE

def main():
    dataset = datasets.load_digits()

    X = dataset.data
    y = dataset.target

    plt.figure(figsize=(12, 8))

    # 主な多様体学習アルゴリズム (と主成分分析)
    manifolders = {
        'PCA': PCA(),
        'MDS': MDS(),
        'Isomap': Isomap(),
        'LLE': LocallyLinearEmbedding(),
        'Laplacian Eigenmaps': SpectralEmbedding(),
        't-SNE': TSNE(),
    }
    for i, (name, manifolder) in enumerate(manifolders.items()):
        plt.subplot(2, 3, i + 1)

        # 多様体学習アルゴリズムを使って教師データを 2 次元に縮約する
        X_transformed = manifolder.fit_transform(X)

        # 縮約した結果を二次元の散布図にプロットする
        for label in np.unique(y):
            plt.title(name)
            plt.scatter(X_transformed[y == label, 0], X_transformed[y == label, 1])

    plt.show()

if __name__ == '__main__':
    main()

上記を適当な名前をつけて保存したら実行する。

$ python manifoldlearning.py

すると、次のようなグラフが得られる。 グラフの各色は、それぞれの数字 (0 ~ 9) を表している。 それぞれのクラスタがキレイにまとまった上で分かれているほど、上手く縮約できているということだと思う。 この中だと t-SNE が頭一つ抜けてるなという感じ。

f:id:momijiame:20171209035139p:plain

次元を縮約した結果を教師データとして用いる

続いては多様体学習を使う目的の二つ目、認識精度の向上について見ていく。 素の教師データ、主成分分析の結果、t-SNE の結果それぞれをランダムフォレストの教師データとして渡してみよう。 汎化性能は、どのように変化するだろうか。

まずは素の教師データから。 digits データセットをそのままランダムフォレストに渡している。 その際の精度をK-分割交差検証で確かめる。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
from sklearn import datasets

from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold

from sklearn.ensemble import RandomForestClassifier

def main():
    dataset = datasets.load_digits()

    X = dataset.data
    y = dataset.target

    scores = np.array([], dtype=np.bool)

    # K-分割交差検証で汎化性能を調べる (分割数は 10)
    cross_validator = KFold(n_splits=10)
    for train, test in cross_validator.split(X):
        X_train, X_test = X[train], X[test]
        y_train, y_test = y[train], y[test]

        # ランダムフォレストで素の教師データを学習する
        cls = RandomForestClassifier()
        cls.fit(X_train, y_train)

        # テストデータに対して分類する
        y_pred = cls.predict(X_test)
        scores = np.hstack((scores, y_test == y_pred))

    # テストデータに対する精度から汎化性能を求める
    accuracy_score = sum(scores) / len(scores)
    print(accuracy_score)


if __name__ == '__main__':
    main()

実行結果は次の通り。 これはランダムフォレストの使う木構造の作られ方にもよるので、出力される数値は毎回微妙に異なる。 今回については約 93% となった。

$ python randomforest.py 
0.929326655537

続いては主成分分析を使って次元を縮約した結果をランダムフォレストに渡すパターン。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
from sklearn import datasets
from sklearn.decomposition import PCA
from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold
from sklearn.ensemble import RandomForestClassifier


def main():
    dataset = datasets.load_digits()

    X = dataset.data
    y = dataset.target

    # 次元縮約に主成分分析を使う
    manifolder = PCA()

    scores = np.array([], dtype=np.bool)

    cross_validator = KFold(n_splits=10)
    for train, test in cross_validator.split(X):
        # 教師データを主成分分析を使って次元縮約する
        X_transformed = manifolder.fit_transform(X)

        X_train, X_test = X_transformed[train], X_transformed[test]
        y_train, y_test = y[train], y[test]

        # 縮約した教師データを学習する
        cls = RandomForestClassifier()
        cls.fit(X_train, y_train)

        y_pred = cls.predict(X_test)
        scores = np.hstack((scores, y_test == y_pred))

    accuracy_score = sum(scores) / len(scores)
    print(accuracy_score)


if __name__ == '__main__':
    main()

実行結果は次の通り。 残念ながら精度は約 87% に低下してしまった。

$ python randomforestpca.py 
0.870339454647

これは、元々のデータセットの性質が非線形なため上手く主成分を取り出すことができなかったということだろう。 まあ、可視化した段階でもそれぞれのクラスタがごちゃっと固まってたしね。

続いては t-SNE を使った場合。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
from sklearn import datasets
from sklearn.manifold import TSNE
from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold
from sklearn.ensemble import RandomForestClassifier


def main():
    dataset = datasets.load_digits()

    X = dataset.data
    y = dataset.target

    # 次元縮約に t-SNE を使う
    manifolder = TSNE()

    scores = np.array([], dtype=np.bool)

    cross_validator = KFold(n_splits=10)
    for train, test in cross_validator.split(X):
        # 教師データを t-SNE で次元縮約する
        X_transformed = manifolder.fit_transform(X)

        X_train, X_test = X_transformed[train], X_transformed[test]
        y_train, y_test = y[train], y[test]

        # 縮約した教師データを学習する
        cls = RandomForestClassifier()
        cls.fit(X_train, y_train)

        y_pred = cls.predict(X_test)
        scores = np.hstack((scores, y_test == y_pred))

    accuracy_score = sum(scores) / len(scores)
    print(accuracy_score)


if __name__ == '__main__':
    main()

実行結果は次の通り。 今度は精度が約 97% に向上した。

$ python randomforesttsne.py 
0.973288814691

素のデータセットをそのまま渡す場合に比べると 4% も認識精度が良くなっている。

t-SNE の高速化

ところで、実際に scikit-learn の t-SNE を使うコードを実行してみると、かなり遅いことに気づくはず。 どうやら、元々 t-SNE は他の多様体学習アルゴリズムに比べると計算量が多いようだ。 ただ、正直このままだと実際のデータセットに使うのは厳しいなあと感じた。 そこで、もうちょっと高速に動作する実装がないか調べたところ Multicore-TSNE という実装があった。

次はこれを試してみよう。 まずはパッケージをインストールする。

$ pip install MulticoreTSNE

次のサンプルコードでは t-SNE の実装を scikit-learn から Multicore-TSNE に切り替えている。 やっていることは scikit-learn 版と変わらない。 コード上の変更点もインポート文を変えたのとジョブ数を指定しているくらい。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import multiprocessing

from matplotlib import pyplot as plt
import numpy as np
from sklearn import datasets
from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold
from sklearn.ensemble import RandomForestClassifier
from MulticoreTSNE import MulticoreTSNE as TSNE

def main():

    dataset = datasets.load_digits()

    X = dataset.data
    y = dataset.target

    manifolder = TSNE(n_jobs=multiprocessing.cpu_count())

    scores = np.array([], dtype=np.bool)

    cross_validator = KFold(n_splits=10)
    for train, test in cross_validator.split(X):
        X_transformed = manifolder.fit_transform(X)

        X_train, X_test = X_transformed[train], X_transformed[test]
        y_train, y_test = y[train], y[test]

        cls = RandomForestClassifier()
        cls.fit(X_train, y_train)

        y_pred = cls.predict(X_test)
        scores = np.hstack((scores, y_test == y_pred))

    accuracy_score = sum(scores) / len(scores)
    print(accuracy_score)


if __name__ == '__main__':
    main()

上記の実行時間は次の通り。 今回使った環境では、だいたい 1 分ほどで終わった。

$ time python multicoretsne.py
0.9816360601
python multicoretsne.py  59.58s user 1.51s system 95% cpu 1:03.96 total

比較として scikit-learn 版の実行時間も示しておく。 こちらは、なんと 6 分ほどかかっている。

$ time python randomforesttsne.py 
0.979410127991
python randomforesttsne.py  382.42s user 34.13s system 96% cpu 7:11.93 total

64 次元 1797 点で 1 分かあ、と思うところはあるものの 6 倍速いという結果は頼もしい限り。

まとめ

今回は多様体学習アルゴリズムを使ってデータセットの次元を縮約してみた。 多様体学習を使って次元を縮約することで、データセットを可視化できたりモデルの認識精度を向上できる場合がある。

Python: NumPy で正方行列を三角行列に加工する

今回は NumPy で正方行列を扱うとき、上三角行列とか下三角行列を取り出す方法について。 三角行列というのは、正方行列において対角要素より上の成分が全て 0 だったり、下の成分が全て 0 だったりする行列のこと。 ちなみに、最初この呼び方を知らなくて「行列 斜め 上」とかでたくさんぐぐった。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.12.6
BuildVersion:   16G1036
$ python --version
Python 3.6.3

インストール

ひとまず、何はともあれ NumPy はインストールしておく。

$ pip install numpy

サンプル用の正方行列を用意する

まずは、加工する前の正方行列を用意する。 今回は 1 ~ 25 の数字を入れた配列を 5 x 5 の正方行列にして使おう。

>>> import numpy as np
>>> array = np.arange(1, 5 * 5 + 1).reshape(5, 5)

中身はこんな感じ。

>>> array
array([[ 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]])

上三角行列 (対角要素あり)

まず、一番簡単なのは対角要素のある上三角行列かな。 これは単純に numpy.triu() 関数に正方行列を渡すだけで得られる。

>>> np.triu(array)
array([[ 1,  2,  3,  4,  5],
       [ 0,  7,  8,  9, 10],
       [ 0,  0, 13, 14, 15],
       [ 0,  0,  0, 19, 20],
       [ 0,  0,  0,  0, 25]])

簡単だね。

上三角行列 (対角要素なし)

続いて対角要素のない上三角行列を。 これは numpy.triu() 関数の k オプションに 1 を渡してやれば良い。

>>> np.triu(array, k=1)
array([[ 0,  2,  3,  4,  5],
       [ 0,  0,  8,  9, 10],
       [ 0,  0,  0, 14, 15],
       [ 0,  0,  0,  0, 20],
       [ 0,  0,  0,  0,  0]])

ちなみに k オプションの値は増やしたり減らすことで削る領域を調整できる。

>>> np.triu(array, k=2)
array([[ 0,  0,  3,  4,  5],
       [ 0,  0,  0,  9, 10],
       [ 0,  0,  0,  0, 15],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0]])
>>> np.triu(array, k=-1)
array([[ 1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10],
       [ 0, 12, 13, 14, 15],
       [ 0,  0, 18, 19, 20],
       [ 0,  0,  0, 24, 25]])

下三角行列 (対角要素あり)

(2017/12/06 追記)

下三角行列の効率的な加工方法を教えて頂きました。 転置行列を使うことで、短く書けるみたいです。 本文も、このやり方を使ったものに変更しました。 大変勉強になりました。この場を借りてお礼を申し上げます。


続いては下三角行列を見ていく。 下三角行列には、転置行列を使うことで効率的に加工できるようだ。

まず NumPy では行列の T アトリビュートから転置行列が得られる。 転置行列では、上三角と下三角が反転した状態になる。

>>> array.T
array([[ 1,  6, 11, 16, 21],
       [ 2,  7, 12, 17, 22],
       [ 3,  8, 13, 18, 23],
       [ 4,  9, 14, 19, 24],
       [ 5, 10, 15, 20, 25]])

そこで NumPy.triu() 関数を使うと上三角だけが残る。 転置する前においては下三角だった要素たちだ。

>>> np.triu(array.T)
array([[ 1,  6, 11, 16, 21],
       [ 0,  7, 12, 17, 22],
       [ 0,  0, 13, 18, 23],
       [ 0,  0,  0, 19, 24],
       [ 0,  0,  0,  0, 25]])

あとは、この状態から、さらに転置行列を取得することで元々の下三角行列が得られる。

>>> np.triu(array.T).T
array([[ 1,  0,  0,  0,  0],
       [ 6,  7,  0,  0,  0],
       [11, 12, 13,  0,  0],
       [16, 17, 18, 19,  0],
       [21, 22, 23, 24, 25]])

下三角行列 (対角要素なし)

対角要素なしの下三角行列の取り出しは、これまでの考え方の組み合わせでいける。

先ほどの考え方はそのままに、上三角行列を取り出すタイミングで k オプションを指定すれば良いだけ。

>>> np.triu(array.T, k=1).T
array([[ 0,  0,  0,  0,  0],
       [ 6,  0,  0,  0,  0],
       [11, 12,  0,  0,  0],
       [16, 17, 18,  0,  0],
       [21, 22, 23, 24,  0]])

めでたしめでたし。

電球形蛍光灯を買ったけどすぐに LED 電球を買い直した話

今回は、珍しくコンピュータではなく照明器具の話題で。 実際にやった (やらかした) のは今から 2 年ほど前のこと。 照明の節電と長寿命化のために電球形蛍光灯を選んでみたら失敗した、という話。

今となっては電球形蛍光灯と LED 電球の価格差はだいぶ縮まってきている。 なので、これから電球形蛍光灯をあえて選ぼうという人も減ってきているはず。 とはいえ、同じ失敗をする人がいるとも限らないので、今さらながら書いてみることにした。

照明を見直したきっかけ

ことの始まりは、今住んでいる賃貸マンションが初期状態で廊下などの照明ソケットに白熱電球を入れていることだった。 ただ、白熱電球は寿命が短いので、だいたい半年から一年くらいするとフィラメントが切れてしまう。 最初こそ同じタイプの白熱電球を買って付け替えていたんだけど、何度も交換するうちにうんざりしてきた。

最初に購入してたのは、こういうやつ。 廊下なんかで使われるソケットは基本的には E26 口金というやつを選べば良いようだ。 以下の商品であれば、定格電力は 54W で定格寿命は 1,500 時間となっている。 Amazon の購入履歴を見ると当時の価格は ¥369 だった。

パナソニック レフ電球(屋内用) E26口金 100V60形 散光形(ビーム角=60°)

パナソニック レフ電球(屋内用) E26口金 100V60形 散光形(ビーム角=60°)

代替製品の検討

続いて、節電と長寿命化を狙って電球形蛍光灯か LED 電球への交換を考えた。 両者を比較すると、消費電力の面でも定格寿命の面でも LED 電球の方が勝っている。 ただ、当時は電球形蛍光灯なら LED 電球の 1/2 ~ 1/3 くらいの値段で買うことができた。 今だと LED 電球の普及が進んだことで、両者の価格差はもっと縮まっている。 とはいえ、同じ会社で同じシリーズならやはり電球形蛍光灯の方が安く買える。

例えば以下の商品であれば、定格電力が 11W で定格寿命が 10,000 時間となっている。 白熱電球に比べると定格電力で約 1/5 だし、定格寿命は約 6 倍になる。

パナソニック 電球形蛍光灯 パルックボール 電球60W形相当 口金直径26mm 電球色相当 EFD15EL11E

パナソニック 電球形蛍光灯 パルックボール 電球60W形相当 口金直径26mm 電球色相当 EFD15EL11E

実際に購入したのは上記の旧モデルだけど、当時は ¥400 で購入できた。 製品の値段が白熱電球と大差ないことを考えると、しめしめこれは良さそうだぞとなった。 しかし、実際に購入して交換してみると、思った以上にデメリットがある。

具体的には、照明が点灯するまでに大きなディレイ (遅れ) があるため。 一般的な蛍光灯を使ったことがある人なら、スイッチを入れてから「・・・チチッ、パッ」という感覚が分かると思う。 それと、完全に明るくなるまでつけっぱなしでしばらく (数分) かかる。 考えてみれば当たり前のことなんだけど、電球形蛍光灯はその特性を完全に引き継いでしまっている。

実際に使うまでは、そんなコンマ数秒の遅れがあっても気にならないだろうと思っていた。 しかし、元々が白熱電球を使っていた場所にそれを交換すると、これがもう気になって仕方がない。 スイッチを入れてから、ほんの少しの間を置かないと明るくならないだけで、こんなにストレスになるとは思わなかった。

LED 電球の買い直し

結局、それから我慢しきれずに LED 電球を買い直すことになった。 例えば、以下の商品なら定格電力が 4.9W で定格寿命が 40,000 時間に及ぶ。 電球形蛍光灯と比べても電力で約半分、寿命は 4 倍になる。 白熱電球と比べれば電力は約 1/11 で、寿命に至っては約 26 倍の計算になる。 何より、使い勝手が白熱電球とほとんど変わらず、点灯までのディレイがないし光量が変化したりもしない。

実際に購入したのは上記の旧モデルだけど、当時は ¥752 だったようだ。 電球形蛍光灯と比較したとき、金額的には大差なくても定格電力の差を考えるとペイするのは相当先のことだった。 それもあって最初、電球形蛍光灯を選んだように思う。 もちろん、白熱電球との比較であれば製品の寿命的にすぐペイできるんだけど。

ところで上記の商品は白熱電球 40 形相当、となっている。 つまり、白熱電球であれば 40W タイプに相当することになる。 あれ、じゃあ最初に出てきた白熱電球よりも暗いの?というと、そんなこともなかった。 どうやら LED 電球は前者二つに比べると明るく感じやすいのか、一つ大きさを落としたくらいでちょうど良いらしい。 実際に購入して、今使っている製品も 40 形相当なんだけど元の白熱電球に比べて暗いという感じは全くない。 また、サイズが小さい方が価格が安いというメリットもある。

ということで、実際に電球形蛍光灯を使っていたのは一ヶ月くらいだったようだ。 安物買いの銭失いとは、正にこのこと。 ちなみに、当たり前だけど交換した LED 電球は 2 年経った今も元気に我が家を照らしている。

ところで「賃貸だから長寿命の照明器具を買ってもなあ」と購入に消極的な考えを持つ必要はないと思う。 というのも、元々ついている電球を、取り外して保管しておけば済む話なので。 別の家に引っ越すときは、保管していたものに付けなおして、自分で購入したものは引越し先に持って行けば良い。 ソケットの規格は日本に住んでいる限りは変わらないだろうから、また使うことができるはず。