CUBE SUGAR CONTAINER

技術系のこと書きます。

統計: 最小二乗法を使って回帰直線を求める

今回は、データセットのある次元の値から別の次元の値を予測する方法の一つとして、最小二乗法というやり方で回帰直線を求めてみる。 とはいえ、いきなり最小二乗法や回帰直線といわれても何が何やらという感じなので、最初はその説明からしていく。

回帰

まず、ある次元の値から別の次元の値を予測するというのは具体的にどういうことだろうか。 例えば、賃貸の部屋の広さと家賃の値段に相関関係があるとしよう。 つまり、広い部屋はそれだけ家賃も高いということで、これは直感に反していない。 だとすれば、部屋の広さと家賃の値段の関係が含まれたデータセットがあるとき、部屋の広さから家賃の値段を予測できるはずだ。

統計や機械学習の世界では、このような問題のジャンルに「教師あり学習」の「回帰」という名前がついている。 ただし、ある次元の値から別の次元の値を予測する問題がすべて回帰と呼ばれるわけではない。 なぜかというと、その次元に含まれるデータの種類によって、やることが異なるため。

データの種類には、大きく分けて二つある。 まず、男性や女性、あるいは好きな食べ物といった、いくつかの種類からひとつを選ぶようなものは質的データという。 それに対し、年齢や一日の摂取カロリーなど、値が連続性を持ったものは量的データという。

ある次元から別の次元を予測するとき、予測する対象が質的データなら「分類」で、量的データなら「回帰」になる。 そして、予測の元ネタになるデータのことを説明変数、予測する値のことを応答変数と呼ぶ。 このとき、説明変数がひとつなら単回帰、複数なら重回帰と呼ばれる。

以上から、今回このエントリでやってみるのは単回帰になる。 回帰直線というのは、説明変数と応答変数の関係を次の一次式で表した直線のこと。

y = \alpha + \beta x

このとき、x が説明変数で y は応答変数を表している。 \alpha\beta は切片と傾きに対応していて、これは回帰係数と呼ばれている。

最小二乗法

最小二乗法というのは、教師データに最もフィットする回帰係数を見つける手法のことをいう。

最小二乗法では、次に示す e_i を最小化するようになっている。 ここで \hat{y_i}y の予測値を表している。 予測値にはハットの記号をつける慣例になっているため。 この e_i のことを残差という。 ようするに、これは真の値と予測値との差のこと。

e_i = y_i - \hat{y_i}

より具体的には、この残差を自乗した合計値が最も小さくなる回帰係数を見つけるのが目的になる。

S(\hat{\alpha}, \hat{\beta}) = \sum_{i = 1}^{N} e_i ^2

これを残差平方和という。 自乗しているのは符号をすべてプラスにするため。

これは e_i を展開すると、こうなる。

S(\hat{\alpha}, \hat{\beta}) = \sum_{i = 1}^{N} ( y_i - \hat{ y_i }) ^2 = \sum_{i = 1}^{N} (y_i - (\hat{\alpha} + \hat{\beta} x_i)) ^2

これを最小にする回帰係数を得る手法を最小二乗法という

具体的な手順

先ほどの式は \hat{\alpha}\hat{\beta} で偏微分して、それが 0 に等しいとすると連立一次方程式に直すことができる、らしい。 が、そこに至るまでの過程はひとまず飛ばして、先人の知恵だけ享受することにする。

まず、必要なのは次の回帰直線を表した一次式。

\hat{y} = \hat{\alpha} + \hat{\beta} x

この中で、まずは \hat{\beta} は、次の式で得られる。

\hat{\beta} = \frac{s_{xy}}{s_x ^2}

ここで、s_{xy}xy の共分散を表している。 つまり、こう。

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

共分散や相関関係については、前回のエントリに詳しく書いたので見てもらいたい。

blog.amedama.jp

また、s_x ^2x の分散を表している。 つまり、こう。

s_x ^2 = \frac{1}{N} \sum_{i = 1}^{N} (xi - \bar{x}) ^2

\bar{x}x の平均値を表すので、こうだね。

\bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i

そして、次に \hat{\alpha} は、次のようにして得られる。

\hat{\alpha} = \bar{y} - \hat{\beta}\bar{x}

\bar{x} と同様に、\bar{y}y の平均値を表している。

\bar{y} = \frac{1}{N} \sum_{i = 1}^{N} y_i

数式がずらずら並んでなんのこっちゃという感じになってきたので、まとめると…

\hat{y} = \hat{\alpha} + \hat{\beta} x

という回帰直線の一次式に含まれる回帰係数 \hat{\alpha}\hat{\beta} は、それぞれ

\hat{\beta} = \frac{s_{xy}}{s_x ^2}

\hat{\alpha} = \bar{y} - \hat{\beta}\bar{x}

という式で得られる、ということ。

実際に計算してみよう

では、上記を Python のプログラムで実際に計算して上手くいくか試してみよう。

サンプルプログラムでは matplotlib を使ってグラフを描いているので、まずはインストールしておく。

$ pip install matplotlib

次が肝心のサンプルプログラム。 中で最小二乗法を使って回帰直線を得ている。

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

from __future__ import division

from matplotlib import pyplot as plt

import numpy as np


def main():
    # データ点数は 100 にする
    N = 100
    # 0 ~ 10 (-5 ~ +5) の範囲で、ランダムに値を散らばらせる
    R = 10

    x = np.arange(N)
    # y = 50 + x の直線に (-5 ~ +5) でランダム性をもたせる
    y = np.arange(N) + np.random.rand(N) * R - R // 2 + 50

    # 平均値と標準偏差
    xmu, xsigma2 = x.mean(), x.var()
    # 平均値
    ymu = y.mean()

    # 共分散
    sxy = np.sum(x * y) / N - xmu * ymu

    # 回帰係数ベータ
    beta = sxy / xsigma2
    # 回帰係数アルファ
    alpha = ymu - beta * xmu

    # 得られた回帰直線の一次式
    print('y = {0} + {1}x'.format(alpha, beta))

    # グラフに回帰直線をひく
    lsm = np.array([alpha + beta * xi for xi in x])
    plt.plot(x, lsm, color='r')

    # 元データもプロットする
    plt.scatter(x, y)

    # グラフを表示する
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

これを lsm.py という名前をつけて実行してみる。

$ python lsm.py 
y = 49.52693410299069 + 0.9984255719901004x

データにはランダム性を持たせているため、微妙に異なるものの、ほぼ y = 50 + x の一次式が得られている!

また、可視化したグラフは次のようになった。

f:id:momijiame:20161008224310p:plain

データに沿って、ちゃんと回帰直線がひけていることが分かる。

まとめ

  • データセットのある次元から別の次元を予測するのに、回帰という問題のジャンルがある
  • 予測するデータの種類によって問題の呼び方は異なる
    • 質的データなら「分類」になる
    • 量的データなら「回帰」になる
  • 予測の元ネタを説明変数、予測する値のことを応答変数という
  • 説明変数がひとつの次元なら単回帰、ふたつ以上なら重回帰という
  • 回帰のやり方のひとつに最小二乗法という手法がある
    • 最小二乗法では残差平方和を最小化する
    • 予測値と真の値との差を残差という
    • 残差の自乗をすべて足したものが残差平方和
  • 説明変数と応答変数の関係を一次式で表した直線を回帰直線という
    • 回帰直線のパラメータとなる \hat{\alpha}\hat{\beta} を回帰係数という

おまけ

回帰係数の \hat{\beta} は式変形すると別の求め方もできる。

これには、回帰係数 r_{xy}xy の標準偏差を使う。

\hat{\beta} = r_{xy} \frac{s_y}{s_x}

回帰係数は、次のように求める。

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

ここで、s_{xy} は共分散なので、次のように求める。

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

s_xs_y は標準偏差なので、こう。

s_x = \sqrt{\frac{1}{N} \sum_{i = 1}^{N} (x_i - \bar{x}) ^2 }

s_y = \sqrt{\frac{1}{N} \sum_{i = 1}^{N} (y_i - \bar{y}) ^2 }

\bar{x}\bar{y} は平均値なので、こう。

\bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i

\bar{y} = \frac{1}{N} \sum_{i = 1}^{N} y_i

これを、先ほどのプログラムに適用してみる。

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

from __future__ import division

from matplotlib import pyplot as plt

import numpy as np


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

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

    sxy = np.sum(x * y) / N - xmu * ymu

    rxy = sxy / (xsigma * ysigma)

    beta = rxy * (ysigma / xsigma)
    alpha = ymu - beta * xmu

    print('y = {0} + {1}x'.format(alpha, beta))

    lsm = np.array([alpha + beta * xi for xi in x])
    plt.plot(x, lsm, color='r')

    plt.scatter(x, y)

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

    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

実行結果は次の通り。

$ python lsm2.py 
y = 50.359765922390736 + 0.9931534165033539x

グラフもこのように。 f:id:momijiame:20161008225924p:plain

どちらでも、大丈夫だね。

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

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

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

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

散布図を描いてみる

それにはまず、散布図を描いてみるという選択肢がある。 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 はコード一行単位でプロファイリング結果が見えるので分かりやすい