CUBE SUGAR CONTAINER

技術系のこと書きます。

統計: 共分散と相関係数でデータセットの相関を調べる

まず、二次元の特徴量をもったデータセットがあるときを考えてみよう。

もし、一方の次元の値が高いときに、もう一方も高い傾向があるときは、両者に正の相関があるという。 反対に、一方の次元の値が高いときに、もう一方は低い傾向があるときは、両者に負の相関があるという。

では、それぞれの次元に正または負の相関があるか否かを調べるには、具体的にどうしたら良いのだろうか。

散布図を描いてみる

それにはまず、散布図を描いてみるという選択肢がある。 x 軸と y 軸に、それぞれの次元の値をプロットするやり方だ。

このとき、例えば正の相関があるなら、値は次のように左下から右上にかけてプロットされる。

f:id:momijiame:20160924161630p:plain

これはつまり x 軸の次元の値が高いときに y 軸の次元の値も高くなることを示す。

反対に、負の相関があるなら、値は次のように左上から右下にかけてプロットされる。

f:id:momijiame:20160924161639p:plain

これはつまり x 軸の次元の値が高いときに y 軸の次元の値は低くなることを示す。

共分散・相関係数という統計量

ただ、散布図を描いただけでは、具体的にどれくらいの相関があるのかが分からない。 そう、相関には強弱がある。

例えば、次のようなふたつの散布図がある。 相関の強弱でいえば、前者の方が後者よりも強い。

f:id:momijiame:20160924161630p:plain

上記の相関は、後者と比較すると強い。

f:id:momijiame:20160924161729p:plain

上記の相関は、前者と比較すると弱い。

では、相関の強弱は具体的にどのようにすれば分かるのだろうか。 それには、相関の強弱を表す統計量を計算することになる。 それが今回紹介する共分散や相関係数といったもの。

共分散

共分散というのは、各次元の値から平均値を引いたもの同士をかけ合わせた上で、総和を取ってデータの点数で割ったもの。 これは、一言で表せば各次元の偏差 (平均値を引いた値) の積 (かけ算した値) の平均値 (総和をデータ点数で割った値) を計算している。

数式で表すと、次のようになる。

{ \displaystyle
s_{xy} = \frac{1}{N} \sum_{i = 1}^{N} (x_i - \bar{x}) (y_i - \bar{y})
}

上記で \bar{x}\bar{y} は、各次元の平均値を表している。

この共分散という統計量は、正の相関が強いほど数値がプラスに大きくなる。 そして、反対に負の相関が強いほど数値はマイナスに大きくなる。 また、相関が弱ければ数値はゼロに近づく。

例えば、先ほどの散布図であれば、前者は共分散が 835 であるのに対し、後者は 820 だった。

なぜ、そうなるのか?

先ほど、共分散の性質として、正の相関が強いほど数値がプラスに大きく、負の相関が強いほど数値がマイナスに大きくなる、と説明した。 ここでは、「そういうものなんだ」と納得してしまうよりも、理屈から理解しておきたい。

そこで、まずは補助線と補足を入れた次の散布図を見てほしい。

f:id:momijiame:20160924162131p:plain

補助線は x 軸と y 軸の平均値を示している。

各次元の補助線にもとづいて、それぞれの点が意味するところを考えてみよう。 点が、もし x 軸の補助線よりも右側にあれば偏差はプラスとなり、そして左側にあればマイナスになることが分かる。 同様に、点がもし y 軸の補助線よりも上側にあれば偏差はプラスとなり、そして下側にあればマイナスになることが分かる。

先ほどの例で示した散布図では、補助線を引くと第一象限と第三象限に点が集まっていた。 それにもとづいて、各象限における偏差と偏差の積の関係を図示してみよう。

f:id:momijiame:20160924162313p:plain

第一象限においては各次元の偏差がプラスになるため、積はプラスになる。 同様に第三象限においては各次元の偏差がマイナスになるため、積はプラスになる。

上記のようなパターンでは偏差の積はほとんどがプラスになることから、その平均値もプラスになることが分かる。

同じように、第二象限と第四象限に点が集まっているパターンについても考えてみよう。

f:id:momijiame:20160924162442p:plain

先ほどの例から、既にだいたい分かると思うんだけど、このときは偏差の積がマイナスになる。 そのため、偏差の積の平均値を取ると、その値がマイナスになることが分かる。

では、値がまんべんなく分布しているパターンではどうなるか。

f:id:momijiame:20160924162529p:plain

このようなときは、偏差の積がプラスだったりマイナスだったりとまちまちになる。 そのため、平均値を計算するとゼロに近づくというわけ。

ということで、共分散を見れば相関が正なのか負なのか分かることが理解できた。

相関係数への拡張

じゃあ共分散さえ見ておけばいつでもオッケーかというと、そうはいかない。 なぜかというと、共分散はそれ単独では相関の強弱が分かりづらいし、異なるデータセットで比較ができない。 どういうことかというと、共分散は元々のデータセットの値の大小 (単位) に影響を受けてしまう。

例えば、カブトムシとクジラの体長と体高の共分散について考えてみよう。 それぞれの値は、カブトムシが数 cm オーダーだとしたら、クジラは数千 cm オーダーになる。 それぞれの偏差の積の平均値を考えると、値の大きさが全く異なるだろうことが分かる。

相関の強度を異なるデータセットで比べるには、まずは単位に依存しない無名数に変換しないといけない。 それが次に紹介する相関係数だ。 これは、共分散を、各次元の標準偏差の積で割ったもの。

数式で表すと次のようになる。

r = \frac{s_{xy}}{s_x s_y}

s_{xy} が共分散で、s_xs_y が各次元の標準偏差。

省略せずに書くと、こう。

{ \displaystyle
r = \frac
{
  \frac{1}{N} \sum_{i = 1}^{N} (x_i - \bar{x}) (y_i - \bar{y})}
{
  \sqrt{\frac{1}{N} \sum_{i = 1}^{N} (x_i - \bar{x}) {}^2}
  \sqrt{\frac{1}{N} \sum_{i = 1}^{N} (y_i - \bar{y}) {}^2}
}
}

分母となる各次元の標準偏差は、二乗して平方根を取ることで符号をすべてプラスにしている。 そのため、共分散と比べるとプラスとマイナスの値が互いに打ち消し合うことがない。 つまり、分母の標準偏差の積は共分散が取りうる最大値となることが分かる。 結果として、相関係数は相関の強弱が -1 から 1 の間で得られることになる。

例えば、共分散で相関の強弱を示したときに使った散布図で、同じように相関係数も計算してみよう。 まず、この散布図では相関係数が 0.9897 として得られた。

f:id:momijiame:20160924161630p:plain

続いて、こちらの場合には相関係数は 0.805 だった。

f:id:momijiame:20160924161729p:plain

相関の強弱が、共分散よりも分かりやすく得られている。

また、次のようにふたつの次元が完全に比例するときは、相関係数が 1.0 になる。

f:id:momijiame:20160924163017p:plain

共分散と相関係数の注意点

ただ、共分散や相関係数でふたつの次元の相関関係が分かるとはいえ、それだけに頼るのは避けた方が良い。

例えば、次の散布図を見てもらいたい。

f:id:momijiame:20160924163634p:plain

上記の散布図には、ふたつの次元に明らかな規則性が見て取れる。

では、上記の共分散と相関係数はどうなるだろうか? なんと、どちらもゼロになるのだ。

共分散と相関係数は、あくまでふたつの次元の間に線形な関係があるか否かしか見ることができない。 先ほどのように、人の目でみればあきらかな法則性があったとして、数値の上ではそれが分からない。 そのため、ふたつの次元の間に関係性を見いだそうとするときは、共分散や相関係数だけを確認して終わることは避ける必要がある。

まとめ

  • ふたつの次元の相関関係の強弱は、共分散や相関係数といった統計量を計算することで分かる
  • ただし、それで分かるのは線形な関係があるか否かだけなので、それだけを確認して終わることは避ける必要がある

おまけ

最初のグラフを描くのに使った Python のソースコードは次の通り。

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

from matplotlib import pyplot as plt

import numpy as np


def main():
    N = 100
    R = 10
    x = np.arange(N) + np.random.rand(N) * R - R // 2
    y = np.arange(N) + np.random.rand(N) * R - R // 2

    xmu, xsigma = x.mean(), x.std()
    ymu, ysigma = y.mean(), y.std()

    covariance = sum([(xi - xmu) * (yi - ymu) for xi, yi in zip(x, y)]) / N
    print('共分散:', covariance)
    correlation_coefficient = covariance / (xsigma * ysigma)
    print('相関係数:', correlation_coefficient)

    plt.scatter(x, y)

    plt.xlabel('x')
    plt.ylabel('y')

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

依存ライブラリとして matplotlib を使っているので、実行する前にインストールする必要がある。

$ pip install matplotlib

統計: 異なるデータセットの標本を標準得点で比較する

例えば、次のような二種類のデータセットをプロットしたヒストグラムがあったとする。 どちらも一般正規分布となっているようだ。

f:id:momijiame:20160902215826p:plain

このデータセットを、あるテストの点数と捉えてみよう。 本来、テストなら上限値と下限値があるはずだけど、そこは簡略化している。

それぞれのデータセット (テストの得点) は次のようなパラメータになっている。 平均値も標準偏差 (値のバラつき) も異なることがわかる。

    • 平均値: 66
    • 標準偏差: 20
    • 平均値: 122
    • 標準偏差: 41

ある人が、このふたつのテストを受けて、それぞれ別々の点数が得られたとする。 そのとき、ふたつの得点がどれだけ優れているのかを比較するには、どうしたら良いだろうか?

例えば A さんは青のテストで 86 点を取って、緑のテストでは 183 点を取ったとする。 どちらも平均値を大きく上回っていることから、良い点数であることは分かる。 しかし、具体的に、どの程度良い点数なのかを知るには、どうしたら良いか。

標準得点で比較する

上記の A さんが取得した得点は、データセット全体 (母集団) から抽出した標本と見なせるはず。 このように、パラメータの異なる分布の標本同士を比べるには、それを標準得点に加工する必要がある。

標準得点というのは、いくつか種類はあるものの、代表的なものに z スコアがある。 z スコアというのは、標本 (x) から平均値 (\mu) を引いて標準偏差 (\sigma) で割ったものをいう。

z = \frac{x - \sigma}{\mu}

この操作を標準化という。 結果として得られるのは、標本が平均値から標準偏差を基準にして何個分離れているかを表したものになる。

標準化の意味

例えば青のデータセットを例にして、z スコアを使った標準化の意味を考えてみよう。

もし、標本が平均値ちょうど (66) なら、最初に平均値を引いた時点で 0 になる。 それを標準偏差 (20) で割っても 0 のままだ。 もし、平均値から標準偏差ひとつ分だけプラスに離れている (66 + 20 = 86) なら、平均値を引くと 20 になる。 それを標準偏差 (20) で割ると 1 だ。

このように、z スコアを使った標準化では、標本が平均値から標準偏差を基準にしていくつ分離れているかが分かる。

z スコアに標準化した標本同士は、比較できる。 なぜなら、どちらも平均値から標準偏差を基準にして、いくつ分離れているかという同じ尺度になっているため。

例を当てはめて比べてみる

例えば、先ほど A さんが取得した得点について考えてみよう。

まず、青のテストで取った 86 点を z スコアに標準化してみる。

z = \frac{x - \sigma}{\mu} = \frac{86 - 66}{20} = 1.0

次に、緑のテストで取った 183 点も z スコアに標準化する。

z = \frac{x - \sigma}{\mu} = \frac{183 - 122}{41} = 1.4878

ふたつの z スコアを比較すると、緑のテストで取った点数の方が高いことが分かった。

1.0 \lt 1.4878

つまり、緑のテストで取った 183 点の方が、青のテストで取った 86 点よりもすごさでいえば上ということが分かった。

z スコアと偏差値

実は、先ほど計算した z スコアは、一度は耳にしたことがあるはずの、とある指数とも深い関わりがある。 それが、受験などでよく使われる偏差値だ。 これも、標準得点のひとつとなる。

偏差値というのは、実は z スコアを少し加工するだけで作ることができる。 具体的には、z スコアを 10 倍して 50 を足したものが偏差値だ。

z スコアというのは、元々のデータセットをすべてこれに変換すると平均値 (\mu) が 0 で標準偏差 (\sigma) が 1 の分布になる。 それに対し、偏差値は平均値 (\mu) が 50 で標準偏差 (\sigma) が 10 になったものをいう。

平均値を 0 から 50 にするには 50 を足せばいい。 標準偏差を 1 から 10 にするには 10 倍すればいい。 つまり、偏差値は z スコアを 10 倍して 50 を足せば得られるということになる。

Z = 10\frac{x - \sigma}{\mu} + 50

ようするに z スコアが 1.0 というのと、偏差値が 60 というのは本質的に同じものを表している。 まあ、後者の方が数値がある程度大きい分、分かりやすいかもしれない。

一応、最初の A さんの得点を偏差値でも比較しておこう。

まずは青の点数から。

Z = 10\frac{x - \sigma}{\mu} + 50 = 10\frac{86 - 66}{20} + 50 = 60

次に緑の点数。

Z = 10\frac{x - \sigma}{\mu} + 50 = 10\frac{183 - 122}{41} + 50 = 64.878

偏差値は z スコアを元にしているので当たり前のことだけど緑の点数の方がすごいということがわかる。

60 \lt 64.878

つまり、テストの平均値と標準偏差さえ分かれば次の式にもとづいて偏差値は自分で計算できる。

Z = 10\frac{x - \sigma}{\mu} + 50

まとめ

  • 異なるデータセットの標本同士を比較するには標準得点を使う
    • 標準得点にはいくつかの種類がある
  • 最も基本的な標準得点は z スコア
    • z スコアは標本から平均値 (\mu) を引いて標準偏差 (\sigma) で割る
  • z スコアは偏差値の元となる
    • z スコアを 10 倍して 50 を足すと偏差値になる

Mac で USB シリアルケーブルの iBUFFALO BSUSRC0610BS を使う

2021-03-23 追記: macOS Mojave (10.14) 以降、ドライバのインストールは不要になっています。

法人向けのネットワーク機器を使ったりするときなんかは、コンソールを取るのにシリアルケーブルが必要なことがある。 今回は、家で使っている iBUFFALO BSUSRC0610BS を Mac OS X で使う方法についてメモっておく。

商品はこちら。

この製品のドライバはチップメーカーが Web サイトで公開している。

www.ftdichip.com

通常であれば Mac OS X 10.9 and above の x64 (64-bit) を選んでインストールすれば良い。

あるいは Homebrew Cask をインストールしているなら、それ経由で入れることもできる。

$ brew cask install ftdi-vcp-driver

インストールしてケーブルをつなぐと /dev 以下にデバイスが見えるようになる。

$ ls /dev | grep usbserial           
cu.usbserial-FTF6PD08
tty.usbserial-FTF6PD08

あとは、このデバイス経由で通信すれば良い。 後ろに指定している 9600 はボーレート (通信速度) なので、機器に設定されているものを使おう。

$ screen /dev/tty.usbserial-FTF6PD08 9600

機器のログインコンソールが表示されれば成功。

login: 

Mac 用のドライバがない USB シリアルケーブル (のチップ) もたまにあったりするので、使えることが確認できているものは安心できる。

めでたしめでたし。

Python: Fabric をスクリプトに組み込んで使う

Fabric は Python で書かれたデプロイやオペレーションを自動化するためのツール。 Fabric では、タスクと呼ばれるオペレーション内容も Python で書く。

今回は、普段ならコマンドラインツールから使うことが多い Fabric を Python のスクリプトに組み込んで使う方法について書く。 尚、Fabric はリモートのホストに接続して使うことが多いため、そのホストとして Vagrant を使って Ubuntu 16.04 LTS のマシンを用意した。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1004
$ python --version
Python 2.7.10

インストール

まずは pip を使って Fabric をインストールしておく。

$ pip install fabric

Vagrant で仮想マシンを用意する

ここでは Vagrant がすでにインストール済みであることを想定する。 もし、入っていないときは公式サイトのバイナリか Homebrew Cask などを使ってインストールしてほしい。

まずは Vagrant の設定ファイル (Vagrantfile) を用意する。 イメージは Chef 社の提供する bento リポジトリにあるものを使った。

$ cat << 'EOF' > Vagrantfile
# -*- mode: ruby -*-
Vagrant.configure("2") do |config|
  config.vm.box = "bento/ubuntu-16.04"
  config.vm.provider "virtualbox" do |vb|
    vb.cpus = "2"
    vb.memory = "1024"
  end
end
EOF

仮想マシンを起動する。

$ vagrant up

起動できたら vagrant ssh-config コマンドを実行しよう。 ここで Port の項目を確認しておく。 おそらく通常は 2222 が割り振られているはず。 ただ、異なる場合には後述する手順で指定するポート番号を変更する必要がある。

$ vagrant ssh-config
Host default
  HostName 127.0.0.1
  User vagrant
  Port 2222
  UserKnownHostsFile /dev/null
  StrictHostKeyChecking no
  PasswordAuthentication no
  IdentityFile /Users/amedama/Documents/vagrant/ubuntu1604/.vagrant/machines/default/virtualbox/private_key
  IdentitiesOnly yes
  LogLevel FATAL

これは OpenSSH の設定ファイルと同じフォーマットになっている。 上記の設定から、先ほど Vagrant で作った仮想マシンには 127.0.0.1:2222 に SSH することでログインできることがわかる。

さて、これで Fabric を使う対象となる仮想マシンの準備ができた。

Fabric の設定ファイルを用意する

Fabric はデフォルトで fabfile という名前のモジュール (またはパッケージ) を、タスクが記述された処理対象として扱う。 なので fabfile.py という名前でファイルを用意しよう。 Python ではスクリプトファイル (*.py) とモジュールが 1:1 で対応している。 例えば fabfile.py であれば Python インタプリタはそれを fabfile というモジュールとして扱うということ。

それでは、今回の主役となる fabfile.py を用意する。 これには task() 関数として Fabric のタスクが定義されている。 この task() 関数はリモートのホスト上で、通常ユーザ権限を使って uptime コマンドを実行することを示している。 また、main() 関数では task 関数を引数にして execute() 関数を実行することもわかる。

$ cat << 'EOF' > fabfile.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from fabric.api import run
from fabric.api import execute

def task():
    run('uptime')


def main():
    execute(task, hosts=['vagrant@127.0.0.1:2222'])


if __name__ == '__main__':
    main()
EOF

普段の使い方

Fabric がインストールされていると fab コマンドが使えるようになる。 この fab コマンド経由でタスクを実行するのが Fabric の普段の使い方だ。

先ほど用意した fabfile.py も fab コマンド経由で実行してみよう。 -H オプションでホスト名、--port オプションでポート名、-u オプションでユーザ名を指定する。 最後の引数は実行するタスク名だ。 仮想マシンのログインパスワードを聞かれるので「vagrant」答えよう。

$ fab -H localhost --port 2222 -u vagrant task
[localhost] Executing task 'task'
[localhost] run: uptime
[localhost] Login password for 'vagrant':
[localhost] out:  10:46:20 up 6 min,  1 user,  load average: 0.00, 0.05, 0.02
[localhost] out:


Done.
Disconnecting from localhost:2222... done.

上手く実行できた。

Python スクリプトとして実行する

先ほどは fab コマンドからタスクを実行した。 今度は Python スクリプトからタスクを実行してみよう。

これは、ただ単に python コマンドで fabfile.py ファイルを実行するだけだ。 先ほどと同じようにログインパスワードを聞かれるので「vagrant」と答える。

$ python fabfile.py
[vagrant@127.0.0.1:2222] Executing task 'task'
[vagrant@127.0.0.1:2222] run: uptime
[vagrant@127.0.0.1:2222] Login password for 'vagrant':
[vagrant@127.0.0.1:2222] out:  10:48:08 up 8 min,  1 user,  load average: 0.00, 0.03, 0.01
[vagrant@127.0.0.1:2222] out:

ちゃんと実行できたことがわかる。

実行した fabfile.py の内容をもう一度掲載しておく。 Python スクリプトとして実行したときは、この中の main() 関数が実行される。 そして、その中では execute() 関数にタスクとログイン名、ホスト名、ポート番号を指定していた。 つまり、組み込みで実行したいときはこの execute() 関数を使えば良いということが分かる。

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

from fabric.api import run
from fabric.api import execute

def task():
    run('uptime')


def main():
    execute(task, hosts=['vagrant@127.0.0.1:2222'])


if __name__ == '__main__':
    main()

環境を設定する

先ほどの例ではパスワードを手動で入力しなければいけなかった。 今度は、そこも自動化してみよう。

次の fabfile.py ではパスワードの手動入力が不要になっている。 ホスト名やユーザ名、パスワードといった Fabric の環境設定は fabric.api.env を指定すれば良い。

$ cat << 'EOF' > fabfile.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from fabric.api import run
from fabric.api import execute
from fabric.api import env

def task():
    run('uptime')


def main():
    env.hosts = ['127.0.0.1:2222']
    env.user = 'vagrant'
    env.password = 'vagrant'
    execute(task)


if __name__ == '__main__':
    main()
EOF

実行してみよう。

$ python fabfile.py
[127.0.0.1:2222] Executing task 'task'
[127.0.0.1:2222] run: uptime
[127.0.0.1:2222] out:  10:48:44 up 8 min,  1 user,  load average: 0.00, 0.03, 0.00
[127.0.0.1:2222] out:

ちゃんと実行できて、今度はパスワードを聞かれることもない。

まとめ

Fabric を Python スクリプトに組み込んで使うときは fabric.api.execute() 関数を使おう。

追記

このやり方にはちょっとした注意点があるため、次の記事で補足している。

blog.amedama.jp

統計: 変動係数で値のバラつきを比べる

まず初めに、次のようなヒストグラムがあったとする。

f:id:momijiame:20160902215826p:plain

このヒストグラムには、青色と緑色のふたつのグループが含まれている。 それぞれのグループは、平均値や度数が異なるようだ。 果たして、それぞれのグループはどちらの方が値のバラつきが大きいのだろうか?

標準偏差だけでは比較できない

通常、データセットの値のバラつきは分散や標準偏差、四分位数といった統計量で表される。 しかし、これらの統計量は、平均値や単位などが異なると単純に比較することはできない。

例えば、あるカブトムシの大きさの標準偏差が 1 cm で、あるクジラの大きさの標準偏差が 1 m だとしよう。 クジラの方が標準偏差にして 100 倍の大きさがある。 しかし、だからといってクジラの方がバラつきも大きいとは限らない。 元々、その生物がだいたいどれくらいの大きさなのかが分からなければ判断がつかない。

先ほどの例であれば、カブトムシがだいたい 10 cm で標準偏差が 1 cm なのと、クジラがだいたい 20 m で標準偏差が 1 m だとしたら? なんとなくカブトムシの方が値のバラつきが大きそうだ、というのが感覚的にも分かる。

変動係数を使う

こういったときは、両者を比較するのに変動係数という統計量を使う。 なんだか大仰な名前がついてるけど、これは単に標準偏差を平均値で割ったもの。 先ほど、感覚的にバラつきの大小を思い浮かべたときは、暗にこの値を比較していたはず。

変動係数の定義としては、こう

C.V. = \frac{\sigma}{\mu}

\sigma は母集団の標準偏差で、\mu は母集団の平均値を表している。

あるいは、対象が標本であれば標準偏差を s で、平均値を \bar{x} で表す。

C.V. = \frac{s}{\bar{x}}

統計の世界では、記号を使い分けることでその意味を伝える。 母集団が、鍋の中に入っているたくさんのスープだとしよう。 標本は、そのスープを味見するためにスプーンですくった一杯を表している。

変動係数で比べてみる

さて、話を変動係数に戻そう。 変動係数の定義通りに、標準偏差を平均値で割ると結果はどうなるだろうか。 これは、無名数といって単位がなくなる。 異なる単位のものを比べるには、この操作が必要になる。

先ほどのカブトムシとクジラの例で計算してみよう。 カブトムシは平均 10 cm で標準偏差が 1 cm だとすると変動係数は 1 / 10 = 0.1 になる。 クジラは平均 20m で標準偏差が 1 m なので変動係数は 1 / 20 = 0.05 だ。 変動係数は値が大きい方がバラつきが大きいことを表している。 感覚的にカブトムシの方がバラつきが大きそうだということを、数値の上でも確かめることができた。

最初に紹介したヒストグラムは?

最初に紹介したヒストグラムの母数 (パラメータ) は、実は次のようになっていた

    • 平均値: 66
    • 標準偏差: 20
    • 度数: 370000
    • 平均値: 122
    • 標準偏差: 41
    • 度数: 520000

ただし、変動係数の定義からすると度数は必要ない。

両方の変動係数を計算すると、まず青は次のようになる

C.V. = \frac{\sigma}{\mu} = \frac{20}{66} \fallingdotseq 0.303

続いて緑

C.V. = \frac{\sigma}{\mu} = \frac{41}{122} \fallingdotseq 0.336

青と緑を比較すると、緑の方が大きい

0.303 \lt 0.336

ということで、緑の方がバラついていることが分かった。

まとめ

異なる平均値や単位をもったデータセットのバラつきを比べるときには変動係数を使おう。

おまけ (その一)

最初のヒストグラムを書くのに使った Python のプログラムは次の通り。

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

import numpy as np

from matplotlib import pyplot as plt


def main():
    x1 = np.random.normal(66, 20, 370000)
    x2 = np.random.normal(122, 41, 520000)

    plt.hist(x1, 100, alpha=0.5)
    plt.hist(x2, 100, alpha=0.5)

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

実行には matplotlib が必要なので pip でインストールする。

$ pip install matplotlib

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G31
$ python --version
Python 3.5.2

おまけ (その二)

ヒストグラムを生成するのに使った母数 (パラメータ) は大学入試センターの資料を使った。 具体的には、青が数学I/数学Aを、緑が英語の内容を、キリの良い数字で使っている。

www.dnc.ac.jp

ただし、実際のデータには最高値があるが、今回のヒストグラムにはない。

Python: line_profiler でボトルネックを調べる

前回は Python の標準ライブラリに用意されているプロファイラの profile/cProfile モジュールについて書いた。

blog.amedama.jp

今回は、同じ決定論的プロファイリングを採用したプロファイラの中でも、サードパーティ製の line_profiler を使ってみることにする。 line_profiler は profile/cProfile モジュールに比べるとコード単位でプロファイリングが取れるところが魅力的なパッケージ。 また、インターフェースについても profile/cProfile に近いものとなっている。

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

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G31
$ python --version
Python 3.5.2

インストール

まずは line_profiler をインストールする。 サードパーティ製のパッケージなので pip を使う。

$ pip install line_profiler

題材とするプログラム

ひとまず、前回も使った FizzBuzz を表示するソースコードをプロファイリングの対象にしよう。

次のようなソースコードを用意する。

$ cat << 'EOF' > fizzbuzz.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


def fizzbuzz(n):
    if n % 3 == 0 and n % 5 == 0:
        return 'FizzBuzz'

    if n % 3 == 0:
        return 'Fizz'

    if n % 5 == 0:
        return 'Buzz'

    return str(n)


def main():
    for i in range(1, 100):
        message = fizzbuzz(i)
        print(message)


if __name__ == '__main__':
    main()
EOF

このプログラムは実行すると 1 から 99 までの FizzBuzz を表示する。

$ python fizzbuzz.py
1
2
Fizz
...(省略)...
97
98
Fizz

プロファイリングしてみる

さて、それでは用意したソースコードを早速プロファイリングしてみよう。

ここでは Python の REPL を使って動作を確認していく。 先ほどのソースコードを用意したのと同じディレクトリで REPL を起動しよう。

$ python

REPL が起動したら、まずはプロファイリング対象の fizzbuzz モジュールをインポートしておく。 補足しておくと Python ではソースコード (*.py) が、そのままモジュールに対応している。 例えば先ほどの fizzbuzz.py であれば fizzbuzz モジュールとして使えるわけだ。

>>> import fizzbuzz

続いて line_prpfiler モジュールをインポートして LineProfiler クラスのインスタンスを生成しよう。 このインスタンスを使ってプロファイリングをする。

>>> import line_profiler
>>> pr = line_profiler.LineProfiler()

プロファイリング対象の関数を LineProfiler#add_function() メソッドを使って登録しよう。 ここでは、前述の fizzbuzz モジュールの中にある fizzbuzz 関数を登録した。

>>> pr.add_function(fizzbuzz.fizzbuzz)

そして、続いて LineProfiler#runcall() メソッドでプロファイリングを実行する。 このメソッドで実行した処理の中で、先ほど LineProfiler#add_function() メソッドで登録した関数が呼び出されることを期待する。 もし呼びだされた場合には、それが処理される過程が line_profiler によって逐一記録されるという寸法だ。

>>> pr.runcall(fizzbuzz.main)
1
2
Fizz
...(省略)...
97
98
Fizz

ちなみに、上記の LineProfiler#runcall() メソッドは、代わりに次のようにしても構わない。 こちらのやり方では LineProfiler#enable() メソッドと LineProfiler#disable() メソッドによってプロファイリングのオン・オフを切り替えている。 LineProfiler#enable() メソッドと LineProfiler#disable() メソッドの間に実行された処理がプロファイリングの対象ということになる。

>>> pr.enable()
>>> fizzbuzz.main()
1
2
Fizz
...(省略)...
97
98
Fizz
>>> pr.disable()

プロファイリングが終わったら LineProfiler#print_stats() メソッドで結果を表示しよう。 すると、次のようにプロファイリング結果が得られる。

>>> pr.print_stats()
Timer unit: 1e-06 s

Total time: 0.000292 s
File: /Users/amedama/Documents/temporary/fizzbuzz.py
Function: fizzbuzz at line 5

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     5                                           def fizzbuzz(n):
     6        99           95      1.0     32.5      if n % 3 == 0 and n % 5 == 0:
     7         6            0      0.0      0.0          return 'FizzBuzz'
     8
     9        93           57      0.6     19.5      if n % 3 == 0:
    10        27           13      0.5      4.5          return 'Fizz'
    11
    12        66           46      0.7     15.8      if n % 5 == 0:
    13        13           29      2.2      9.9          return 'Buzz'
    14
    15        53           52      1.0     17.8      return str(n)

上記の表のカラムは次のような意味がある。

  • Line

    • ソースコードの行番号
  • Hits

    • その行が実行された回数
  • Time

    • その行の実行に費やした時間の合計
    • 単位は先頭行の Timer unit で示されている
  • Per Hit

    • その行を 1 回実行するのに費やした時間
  • % Time

    • その行の実行に費やした時間の全体に対する割合
  • Line Contents

    • 行番号に対応するソースコード

とても分かりやすい。

コマンドラインから調べる

line_profiler には、もうひとつコマンドラインからプロファイリングを実行するインターフェースがある。 むしろ、公式の説明を見るとこちらを推している雰囲気がある。 ただ、個人的には先ほどのやり方が便利なので、こちらをあまり使うことはないかな…という感じ。

コマンドラインからの処理では、まずプロファイリングを実行したい関数に一工夫が必要になる。 具体的には、プロファイリングを実行したい関数やメソッドを @profile デコレータで修飾する。

以下は fizzbuzz() 関数に @profile デコレータをつけたサンプルコード。

$ cat << 'EOF' > fizzbuzz.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


@profile
def fizzbuzz(n):
    if n % 3 == 0 and n % 5 == 0:
        return 'FizzBuzz'

    if n % 3 == 0:
        return 'Fizz'

    if n % 5 == 0:
        return 'Buzz'

    return str(n)


def main():
    for i in range(1, 100):
        message = fizzbuzz(i)
        print(message)


if __name__ == '__main__':
    main()
EOF

ちなみに、するどい人は上記が @profile デコレータがインポートされていないのでエラーになると気づくかもしれない。

その通りで、実はこのソースコードは通常の Python インタプリタでは実行できなくなる。

$ python fizzbuzz.py 
Traceback (most recent call last):
  File "fizzbuzz.py", line 5, in <module>
    @profile
NameError: name 'profile' is not defined

ただ、line_profiler を使ってコマンドラインでプロファイリングするときは専用のコマンドを使うことになる。 具体的には kernprof というコマンドを使う。 このコマンドではソースコードの中の @profile デコレータがちゃんと動作するように細工して実行するので問題ない。

早速 kernprof コマンドで fizzbuzz.py をプロファイリングしてみよう。

$ kernprof -l fizzbuzz.py
1
2
Fizz
...(省略)...
97
98
Fizz
Wrote profile results to fizzbuzz.py.lprof

実行すると、プロファイリング対象のファイル名に .lprof という拡張子のついたファイルができる。

そして、そのできたファイルを line_profiler モジュールをコマンドライン経由で読み込ませる。 すると、解析結果が表示される。

$ python -m line_profiler fizzbuzz.py.lprof
Timer unit: 1e-06 s

Total time: 0.000275 s
File: fizzbuzz.py
Function: fizzbuzz at line 5

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     5                                           @profile
     6                                           def fizzbuzz(n):
     7        99          116      1.2     42.2      if n % 3 == 0 and n % 5 == 0:
     8         6            3      0.5      1.1          return 'FizzBuzz'
     9
    10        93           53      0.6     19.3      if n % 3 == 0:
    11        27           15      0.6      5.5          return 'Fizz'
    12
    13        66           33      0.5     12.0      if n % 5 == 0:
    14        13            5      0.4      1.8          return 'Buzz'
    15
    16        53           50      0.9     18.2      return str(n)

ただ、上記のやり方だとソースコードに手を加えないといけないから面倒くさい。 それに、モジュールが実行可能になっていないとプロファイリングできないところも汎用性が低いと感じた。

まとめ

  • 標準ライブラリには profile/cProfile というプロファイラがある
  • それ以外にもサードパーティ製で line_profiler というパッケージがある
  • line_profiler はコード一行単位でプロファイリング結果が見えるので分かりやすい

Python: profile/cProfile モジュールでボトルネックを調べる

プログラムを作ってみたは良いけれど、想定していたよりもパフォーマンスが出ないという事態はよくある。 そんなときはパフォーマンス・チューニングをするわけだけど、それにはまず測定が必要になる。 具体的に何を測定するかというと、プログラムの中のどこがボトルネックになっているかという点だ。 この測定の作業は一般にプロファイリングと呼ばれる。

プロファイリングにはふたつの種類がある。 ひとつ目は今回紹介する profile/cProfile モジュールが採用している決定論的プロファイリング。 そして、もうひとつが統計的プロファイリングだ。

決定論的プロファイリングでは、プログラムが実行されていく過程を逐一記録していく。 そのため、時間のかかっている処理を正確に把握することができる。 ただし、測定自体がプログラムに与える影響 (プローブ・エフェクト) も大きくなる。

それに対し統計的プロファイリングでは、プログラムが実行されていく過程をランダムに記録していく。 これは、統計的にいえばプログラムが実行されていく過程の全体を母集団と捉えて、そこから標本を抽出している。 そして、抽出した標本から母集団を推定する形で、呼び出しの多い処理や時間のかかっている処理を推定することになる。 この場合、結果はあくまで推定 (一定の割合で誤りを許容すること) になるものの、測定自体が与える影響は比較的小さくなる。

profile/cProfile モジュール

前置きが長くなったけど、今回は Python の標準ライブラリに用意されているプロファイリング用のモジュールを紹介する。 具体的には profile と cProfile というふたつの実装で、どちらも前述した通り決定論的プロファイリングを採用している。 ふたつある理由は profile が Pure Python の実装で、cProfile が C 拡張モジュールを使った実装になっている。 それぞれの特徴は、profile に比べると cProfile は測定にかかるオーバーヘッドが少ない代わりに、ユーザの拡張性が低くなっている。 使い分けについては、一般的なユースケースでは cProfile を常に使っていれば良さそうだ。

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

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G31
$ python --version
Python 3.5.1

まずは使ってみる

何はともあれ profile/cProfile を使ってみることにしよう。 まずは Python の REPL を使って動作を確認していくことにする。

$ python

起動したら cProfile をインポートする。

>>> import cProfile

もちろん profile を使っても構わない。

>>> import profile

動作確認には time モジュールを使うことにしよう。 このモジュールを使うとスレッドを一定時間停止させることができる。

>>> import time

最もシンプルな使い方としてはモジュールの run() 関数に呼び出したい Python のコードを渡す形になる。 こうすると、渡した Python のコードが cProfile 経由で実行されて、プロファイリング結果が得られる。

>>> cProfile.run('time.sleep(1)')
         4 function calls in 1.001 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    1.001    1.001 <string>:1(<module>)
        1    0.000    0.000    1.001    1.001 {built-in method builtins.exec}
        1    1.001    1.001    1.001    1.001 {built-in method time.sleep}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

結果として得られる表の読み方は次の通り。

  • ncalls

    • 呼び出し回数
  • tottime

    • その関数の実行にかかった時間の合計
    • (別の関数呼び出しにかかった時間を除く)
  • percall

    • かかった時間を呼び出し回数で割ったもの
    • (つまり 1 回の呼び出しにかかった平均時間)
  • cumtime

    • その関数の実行にかかった時間の合計
    • (別の関数呼び出しにかかった時間を含む)
  • filename

    • lineno(function): 関数のファイル名、行番号、関数名

再帰的な処理が含まれる場合

上記の ncalls (呼び出し回数) は、再帰呼び出しがあるときとないときで表示が異なる。 例えば、次のような再帰を使ってフィボナッチ数を計算する関数を用意する。

>>> def fib(n):
...     if n == 0:
...         return 0
...     if n == 1:
...         return 1
...     return fib(n - 2) + fib(n - 1)
...

これをプロファイリングにかけると 177/1 といったように分数のような表示が得られる。

>>> cProfile.run('fib(10)')
         180 function calls (4 primitive calls) in 0.000 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    177/1    0.000    0.000    0.000    0.000 <stdin>:1(fib)
        1    0.000    0.000    0.000    0.000 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

これは、元々は 1 回の呼び出しから関数が再帰的に 176 回呼びだされたことを意味している。 右側が再帰によるものでない呼び出し回数、左側は再帰を含む呼び出し回数になっているため。

より実践的な使い方

先ほどの例では Python のコードを文字列で渡すことでプロファイリングした。 とはいえ、より柔軟に一連の Python のコードをプロファイリングしたいという場合も多いはず。 次は、そんなときはどうしたら良いかについて。

まずは Profile クラスのインスタンスを用意する。

>>> pr = cProfile.Profile()

そして、プロファイリングを開始するために enable() メソッドを呼びだそう。

>>> pr.enable()

その上で、プロファイリングしたいコードを実行する。

>>> fib(10)
55

そして、プロファイリングを終了するために disable() メソッドを呼び出す。

>>> pr.disable()

ちなみに、これは代わりに runcall() メソッドを使っても構わない。

>>> pr.runcall(fib, 10)
55

プロファイリングが終わったら print_stats() メソッドで結果を出力する。

>>> pr.print_stats()
         179 function calls (3 primitive calls) in 0.000 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.000    0.000 <stdin>:1(<module>)
    177/1    0.000    0.000    0.000    0.000 <stdin>:1(fib)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

結果を詳しく見る

これまでの例では、プロファイリングの対象となるコードがシンプルだったので結果もシンプルだった。 とはいえ、実際にプロファイリングしたいようなコードは色々な処理が呼び出しあうような複雑なものになるはず。 結果も自ずと重厚長大なものになるはずだ。 そんなときは、プロファイリングの結果を pstats モジュールで詳しく調べることができる。

まずは pstats モジュールをインポートする。

>>> import pstats

そして、先ほどの Profile クラスのインスタンスを渡して Stats クラスをインスタンス化しよう。

>>> stats = pstats.Stats(pr)

Stats クラスのインスタンスでは、例えば特定のカラムで結果をソートできる。 ここでは ncalls (呼び出し回数) で降順ソートしてみよう。

>>> stats.sort_stats('ncalls')
<pstats.Stats object at 0x102518278>

>>> stats.print_stats()
         178 function calls (2 primitive calls) in 0.000 seconds

   Ordered by: call count

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    177/1    0.000    0.000    0.000    0.000 <stdin>:1(fib)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


<pstats.Stats object at 0x102518278>

これで、呼び出し回数の多い fib() 関数が上にくる。

プロファイリング結果をファイルに残す

プロファイリングの結果はデータとして残せるようになっていると色々と便利なはず。

そんなときは Profile クラスの dump_stats() というメソッドが使える。 このメソッドは、ファイル名を指定してプロファイリング結果をファイルに保存できる。

>>> pr.dump_stats('fib.profile')

上記であれば Python の REPL を起動したディレクトリに fib.profile という名前のファイルができる。

$ ls | grep ^fib
fib.profile

上記のファイルは pstats モジュールから読み込んで使うことができる。 Stats クラスをインスタンス化するときにファイル名を文字列で指定すれば結果を読み込むことができる。

>>> import pstats
>>> stats = pstats.Stats('fib.profile')
>>> stats.print_stats()
Tue Aug 30 19:30:07 2016    fib.profile

         178 function calls (2 primitive calls) in 0.000 seconds

   Random listing order was used

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    177/1    0.000    0.000    0.000    0.000 <stdin>:1(fib)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


<pstats.Stats object at 0x10ec9acf8>

コマンドラインでプロファイリングする

次は、そんなに使う場面はないかもしれないけどコマンドラインでプロファイリングする方法について。

例えば、次のような FizzBuzz を実行するソースコードを用意しておこう。

$ cat << 'EOF' > fizzbuzz.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


def fizzbuzz(n):
    if n % 3 == 0 and n % 5 == 0:
        return 'FizzBuzz'

    if n % 3 == 0:
        return 'Fizz'

    if n % 5 == 0:
        return 'Buzz'

    return str(n)


def main():
    for i in range(1, 100):
        message = fizzbuzz(i)
        print(message)


if __name__ == '__main__':
    main()
EOF

このソースコードは実行すると 1 から 99 までの FizzBuzz を出力する。

$ python fizzbuzz.py 
1
2
Fizz
...(省略)...
97
98
Fizz

このとき -m cProfile というオプションをつけて cProfile モジュールが実行されるようにしてみよう。 これで、自動的にプロファイリングが実行される。

$ python -m cProfile fizzbuzz.py
1
2
Fizz
...(省略)...
97
98
Fizz
         202 function calls in 0.001 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.001    0.001 fizzbuzz.py:18(main)
        1    0.000    0.000    0.001    0.001 fizzbuzz.py:5(<module>)
       99    0.000    0.000    0.000    0.000 fizzbuzz.py:5(fizzbuzz)
        1    0.000    0.000    0.001    0.001 {built-in method builtins.exec}
       99    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

この使い方では -s オプションを使って結果のソートができる。

$ python -m cProfile -s ncalls fizzbuzz.py
1
2
Fizz
...(省略)...
97
98
Fizz
         202 function calls in 0.001 seconds

   Ordered by: call count

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       99    0.000    0.000    0.000    0.000 {built-in method builtins.print}
       99    0.000    0.000    0.000    0.000 fizzbuzz.py:5(fizzbuzz)
        1    0.000    0.000    0.001    0.001 fizzbuzz.py:5(<module>)
        1    0.000    0.000    0.001    0.001 fizzbuzz.py:18(main)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        1    0.000    0.000    0.001    0.001 {built-in method builtins.exec}

同じように -o オプションを使えばプロファイリング結果をファイルに保存できる。

$ python -m cProfile -o fizzbuzz.profile fizzbuzz.py
$ ls | grep .profile$
fizzbuzz.profile

呼び出し関係を調べる

profile/cProfile は処理の呼び出し回数や時間などの他に、呼び出し関係も調べることができる。 例えば、先ほどの FizzBuzz のプロファイリング結果を例に使ってみよう。

まずは改めて Python の REPL を起動する。

$ python

そしてプロファイリング結果をファイル経由で読み込む。

>>> import pstats
>>> p = pstats.Stats('fizzbuzz.profile')

そして print_callers() メソッドで、呼び出し元が分かる。 左が呼び出された関数で、右が呼び出し元の関数になっている。

>>> p.print_callers()
   Random listing order was used

Function                                          was called by...
                                                      ncalls  tottime  cumtime
{built-in method builtins.print}                  <-      99    0.000    0.000  fizzbuzz.py:18(main)
fizzbuzz.py:5(<module>)                           <-       1    0.000    0.001  {built-in method builtins.exec}
{built-in method builtins.exec}                   <-
{method 'disable' of '_lsprof.Profiler' objects}  <-
fizzbuzz.py:18(main)                              <-       1    0.000    0.001  fizzbuzz.py:5(<module>)
fizzbuzz.py:5(fizzbuzz)                           <-      99    0.000    0.000  fizzbuzz.py:18(main)


<pstats.Stats object at 0x100d24a58>

その逆が print_callees() メソッドで得られる。 今度は左が呼び出し元で、右が呼び出し先になっている。

>>> p.print_callees()
   Random listing order was used

Function                                          called...
                                                      ncalls  tottime  cumtime
{built-in method builtins.print}                  ->
fizzbuzz.py:5(<module>)                           ->       1    0.000    0.001  fizzbuzz.py:18(main)
{built-in method builtins.exec}                   ->       1    0.000    0.001  fizzbuzz.py:5(<module>)
{method 'disable' of '_lsprof.Profiler' objects}  ->
fizzbuzz.py:18(main)                              ->      99    0.000    0.000  fizzbuzz.py:5(fizzbuzz)
                                                          99    0.000    0.000  {built-in method builtins.print}
fizzbuzz.py:5(fizzbuzz)                           ->


<pstats.Stats object at 0x100d24a58>

まとめ

  • パフォーマンスチューニングには、まずプロファイリングが必要になる
  • プロファイリングには決定論的プロファイリングと統計的プロファイリングのふたつがある
  • Python には profile/cProfile という決定論的プロファイリングができるモジュールがある
  • このモジュールを使うと特定の処理における関数の呼び出し回数や、かかった時間を調べることができる

参考

27.4. Python プロファイラ — Python 3.6.5 ドキュメント