2017-01-19

Mac: コマンドラインでプリンタを操作する

Mac OS X

ちょっとまとまった量の書類を印刷する機会があって、手作業は大変だからターミナルで作業したいなと思った。今回は、そのとき調べた内容について書いておく。尚、操作するプリンタのドライバ等は既にインストールした状態を前提にしている。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1217

プリンタの一覧を取得する

まずは、操作するプリンタの名称を取得する必要がある。これには lpstat コマンドに -s オプションをつけて実行する。

$ lpstat -s
システムのデフォルトの送信先：Canon_iP2700_series
Canon_iP2700_series のデバイス：usb://Canon/iP2700%20series?serial=XXXXXX

上記の出力で「〜のデバイス」となっている中の「〜」がプリンタの名前になっている。例えば、上記なら「Canon_iP2700_series」がそれに当たる。

プリンタを使って印刷する

実際にプリンタを使って印刷するときは lpr コマンドを使う。印刷に使うプリンタは -P オプションで指定する。つけなければデフォルトの送信先になるかな。

$ lpr -P <プリンタ名> <ファイル名>

例えば、先ほどのプリンタを例にすると、こんな感じ。

$ lpr -P Canon_iP2700_series sample.pdf

コマンドラインで実行できれば、あとは好きに操作を自動化できる。例えば現在の作業ディレクトリにあるファイル (pdf とか) をすべて印刷したいなら、次のようにしたり。

$ for i in $(ls); do lpr -P Canon_iP2700_series $i; done;

プリンタの設定を変更する

プリンタの設定は lpoptions コマンドで行う。使っているプリンタで設定できる項目の一覧を取得したいときは -l オプションをつけて実行しよう。対象のプリンタは -p オプションで指定する。

$ lpoptions -p Canon_iP2700_series -l
ColorModel/Color Model: *RGB16
CNIJProfileID/ProfileID: *1 2 3 4 5 6
Resolution/Output Resolution: 300x300dpi *600x600dpi
PageSize/PageSize: Letter Letter.FullBleed Legal A5 *A4 A4.FullBleed B5 4x6 4x6.FullBleed 4x8 4x8.FullBleed 5x7 5x7.FullBleed 8x10 8x10.FullBleed 89x127mm 89x127mm.FullBleed Postcard Postcard.FullBleed DoublePostcard Env10 EnvDL EnvYou4 98x190mm 55x91mm 55x91mm.FullBleed 101.6x180.62mm 101.6x180.62mm.FullBleed Custom.WIDTHxHEIGHT
CNIJCartridge/BJ Cartridge: 1 2 *4
CNIJMediaType/Media Type: *0 50 52 51 42 22 28 16 56 27 54 7 18 8 36
CNIJMediaSupply/Paper Source: *7
CNIJPrintQuality/Quality: 0 5 *10 15 20
CNIJDitherPattern/Halftoning: 2560 *2595
CNIJGrayScale/Grayscale Printing: *0 1 2
CNIJGamma2/Brightness: 14 *18 22
CNIJAmountOfExtension/Amount of Extension: 0 1 *2 3
CNIJMarginType/MarginType: *0 1
CNIJPaperGapCommand/PaperGapCommand: *10
CNIJDiscTrayForm/DiscTrayForm: *10
CNIJBanner/Banner: *5
CNIJDuplexPrinting/Duplex Printing: *0
CNIJSGColorMode/SGColorMode: *127
CNIJISaveMode/ISaveMode: *30
CNIJMediaGroup/MediaGroup: *99
CNIJDuplexSurface/DuplexSurface: *0
CNIJDataDirection/DataDirection: *0
CNIJSpecialMode/SpecialMode: *0
CNIJManualSelectSupply/ManualSelectSupply: *10
CNIJIntent2/Color Mode: *1 3
CNIJBlackAdjustment/BlackAdjustment: *10
CNIJEmergencyMode/EmergencyMode: *127
CNIJLongLifeMode/LongLifeMode: *30
CNIJCartridgeSelect/CartridgeSelect: *1
CNIJPrintMode2/Print Quality: 1 *2 3 5
CNIJPZVividPosiProcess/Vivid Photo: *0 2048
CNIJGrayScaleCheckBox/Grayscale Printing: *0 1
CNIJPQualitySlider/Quality: 1 2 *3 4 5
CNIJHalfToneRadio/Halftoning: 1 *2
CNIJMediaSupplyMenuItem/MediaSupplyMenuItem: *1
CNIJDuplexPrintArea/Print Area: 0 *1
CNIJProfileType/ProfileType: 0 *1
CNIJPDEEnable/PDEEnable: 0 *1
CNIJRGB2GrayConvert/RGB2GrayConvert: *0 1
CNIJAutoDuplexCheckBox/AutoDuplexCheckBox: *0 1
CNIJPreviewSampleType/Sample Type: *1 3 4
CNIJColorPatternCheckBox/View Color Pattern: *0 1
CNIJColorAdjustMode/ColorAdjustMode: *1
CNIJMediaSense/MediaSense: *10
CNIJStapleSide/Stapling Side: *1 2 3 4
CNIJExtDDIStapleMarginSupport/ExtDDIStapleMarginSupport: *1
CNIJExtDDISupportAppType/ExtDDISupportAppType: *1
CNIJTableSwitchInfo/TableSwitchInfo: *0

上記はプリンタ固有の設定項目なので、使用するプリンタの機種ごとに異なる。

モノクロ印刷の設定をしてみる

設定例として、カラーインクがもったいないのでモノクロ印刷がしたい場合を考えてみよう。そのときは、まず関連しそうな項目を上記の出力から探す。

今回使っているプリンタでは GrayScale という文字列の入った項目が該当しそうだ。例えば CNIJGrayScale/Grayscale Printing という項目には 0, 1, 2 という三つの設定値がある。

$ lpoptions -p Canon_iP2700_series -l | grep -i gray
CNIJGrayScale/Grayscale Printing: *0 1 2
CNIJGrayScaleCheckBox/Grayscale Printing: *0 1
CNIJRGB2GrayConvert/RGB2GrayConvert: *0 1

アスタリスク (*) がついているのが現在の設定なので、今は 0 になっていることが分かる。

試しに、これを 1 に変更してみよう。設定を変更するときは -o オプションで <項目名>=<値> という形で指定する。

$ lpoptions -p Canon_iP2700_series -o CNIJGrayScale=1

これで CNIJGrayScale/Grayscale Printing の項目の値が 1 に変更できた。

$ lpoptions -p Canon_iP2700_series -l | grep CNIJGrayScale
CNIJGrayScale/Grayscale Printing: 0 *1 2
CNIJGrayScaleCheckBox/Grayscale Printing: *0 1

あとは、この状態で印刷すればモノクロで印刷できた。

$ lpr -P Canon_iP2700_series sample.pdf

めでたしめでたし。

2017-01-14

統計: F 分布を使って二つの標本の分散が等しいか調べる

統計

統計の世界には、二つの標本から得られた分散が等しいかそうでないかを確かめるための手法がある。それが、今回紹介する F 分布と、それを用いた F 検定だ。なぜ、そんなものがあるかというと、統計には二つの標本を比べるときに分散が等しいかそうでないかでやり方を変える必要のある手法があるから。例えば二標本間で母平均の差を区間推定もしくは差の有無を検定するときに使われる t 検定がそれに当たる。

F 検定のやり方

検定するといっても、今回は分かりやすさを優先して、ここでは帰無仮説と対立仮説がどうの、という話はしないことにする。代わりに、結果から「差がある」といえるか「差がない」といえるのかを文章で解説していく。

F 分布を使った検定のやり方はとてもシンプル。まずは、次のようにして F 値を計算する。

$F = \frac{s_1^{2}}{s_2^{2}}$

ここで $s_1^{2}$ と $s_2^{2}$ は、それぞれ別々の標本から計算された不偏分散を表している。注意点としては、式の分子の方に分散の大きい値を持ってくるという点だ。

念のため、不偏分散の計算式も示しておく。

$s = \frac{1}{n - 1} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}$

ここで $\bar{x}$ は標本平均を表しているため、こう。

$\bar{x} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

次に、F 分布を参照するために必要な自由度を確認する。自由度は、それぞれの標本のデータ数から 1 を引いた数になる。

$df_1 = n_1 - 1$

$df_2 = n_2 - 1$

あとは、上記の二つの自由度と有意水準を元に F 分布表を参照する。有意水準というのは、その手法を使ったとき、どれくらいの割合で間違った結論を導くかを示す値をいう。有意水準には一般に 0.05 や 0.01 を用いる。これは、例えば 0.05 であれば 5% の確率で間違った結論を導いてしまう、ということを意味している。

分布表(t分布表・X2分布表・標準正規分布表・F分布表・Wilcoxon)

このとき、最初に計算した F 値が F 分布表の値を超えるようであれば、二つの標本から計算した分散には有意な差があるといえる。つまり、分散は等しくないということだ。反対に、もし F 値が F 分布表の値を超えないようなら、両者に有意な差はないということになる。

例

惣菜屋甲と惣菜屋乙はおにぎりを販売している。 A さんは、甲で販売しているおにぎりが、乙で販売しているおにぎりと比べて、重さのばらつきが大きいと感じている。ある日、甲でおにぎりを 11 個、乙でおにぎりを 16 個買い、重さを測ったところ標本平均と不偏分散は次のようになった。

惣菜屋甲: 標本平均 = 118.0 不偏分散 = 9.0

惣菜屋乙: 標本平均 = 120.0 不偏分散 = 4.0

おにぎりの重さの母集団分布はいずれの惣菜屋でも正規分布で近似できるものとする。

(統計検定2級公式問題集 2015年6月分問14より)

上記を有意水準 0.05 で実際に検定してみよう。

まずは両者の不偏分散から F 値を計算する。ポイントは前述した通り、分子に大きい値を持ってくるところ。ここでは、惣菜屋甲のおにぎりの方が分散が大きいので、分子に持ってくる。

$F = \frac{s_1^{2}}{s_2^{2}} = \frac{9.0}{4.0} = 2.25$

次に、F 分布表を参照するために自由度を求めよう。まず、甲の自由度は標本数から 10 になる。

$df_1 = 11 - 1 = 10$

そして、乙の自由度は標本数から 15 になる。

$df_2 = 16 - 1 = 15$

有意水準 0.05 と、各自由度を元に F 分布表を参照する。

分布表(t分布表・X2分布表・標準正規分布表・F分布表・Wilcoxon)

上記で $df_1 = 10, df_2 = 15$ の値を参照すると 2.54 と分かる。

$F_{0.05} = 2.54$

F 分布表を参照した値と、先ほど計算した F 値と比べると、F 値の方が小さいことが分かる。

$F = 2.25 \le F_{0.05} = 2.5437$

上記から、惣菜屋甲と惣菜屋乙が作るおにぎりの分散に有意な差は見られないことが分かった。

まとめ

二つの標本の分散が等しいかを調べるのには F 検定が使える。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

作者: もみじあめ
発売日: 2016/03/12
メディア: Kindle版
この商品を含むブログを見る

2017-01-03

統計: 続・はじめての推定 (母平均編)

統計

今回のエントリは、以下のエントリの続きになっている。

blog.amedama.jp

上記のエントリでは、統計における基本的な推定の考えを解説した。また、その例として考えうる最も単純な推定を扱っている。それというのは、母集団の全ての母数 (パラメータ) が既知な状況での標本の推定だった。

今回は、それの続きとして母集団の母数 (パラメータ) が部分的に未知な状況を扱う。幾つかのパターンを扱う中で一貫しているのは、母平均が未知になっているのでそれを推定するという点だ。

状況ごとに使う手法について

まず、全ての前提として母集団は正規分布ということが分かっているとする。しかし、母集団の平均は分かっておらず、これを推定したい。

このとき、点推定に関しては状況によらず標本平均を使えば良い。これが母平均の点推定になる。

$\bar{x} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

それに対し、区間推定のときは状況によって、使う手法を三つのパターンで変えなければいけない。

母分散が既知のとき
- 母分散と正規分布を使って推定する
母分散が未知で大標本のとき
- 不偏分散と正規分布を使って推定する
母分散が未知で小標本のとき
- 不偏分散と t 分布を使って推定する

上記で、実用的なのは母分散が未知なときの推定だろう。なぜなら、母分散が既知なのに母平均は未知なので知りたいという状況は、そうそうないため。

母分散が既知のとき

母分散が既知のときは、実のところ前回のエントリで使ったのと同じ以下の不等式が使える。

$-z_{\frac{\alpha}{2}} \le \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}} \le z_{\frac{\alpha}{2}}$

上記において $z_{\frac{\alpha}{2}}$ は有意水準にもとづいた値になる。例えば信頼係数 95% なら 1.96 だし 99% なら 2.58 を使う。それ以外では $\bar{x}$ が標本平均で $\sigma$ が母標準偏差、そして $\mu$ が母平均になる。

上記の式で未知なのは母平均 $\mu$ だけだ。母標準偏差 $\sigma$ は母分散の平方根を取れば良いし、標本平均 $\bar{x}$ はいくつか取り出して計算すれば分かる。なので、やることは母平均 $\mu$ について上記の不等式を解くだけで良い。

まずは三辺に $\frac{\sigma}{\sqrt{n}}$ を掛ける。

$-z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}} \le \bar{x} - \mu \le z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}$

次に三辺に $\mu$ を足す。

$\mu - z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}} \le \bar{x} \le \mu + z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}$

左辺に着目する。

$\mu - z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}} \le \bar{x}$

$\mu$ について解く。

$\mu \le \bar{x} + z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}$

同じように右辺に着目する。

$\bar{x} \le \mu + z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}$

$\mu$ について解く。

$\bar{x} - z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}} \le \mu$

二つの式をつなげると、こうなる。

$\bar{x} - z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}} \le \mu \le \bar{x} + z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}$

つまり、母分散が既知なときの区間推定は次の式で求まる。

$\bar{x} \pm z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}$

例題

あるおにぎり製造マシンは作るおにぎりの重さを自由に調整できる。このとき、設定した重さと実際に出来上がる重さには誤差が出るが、その分布は正規分布に従うとする。また、母分散は作るおにぎりの重さに関わらず 9g とカタログスペックから分かっている。この機械を使って試しにおにぎりを 16 個作ってみたところ、その平均は 100g だった。この機械が作るおにぎりの母平均を信頼係数 95% で区間推定せよ。

問題文から得られた値を、先ほどの式に代入する。

$100 \pm 1.96 \frac{\sqrt{9}}{\sqrt{16}}$

上記を解くと、次の区間が得られる。

$[ 98.53, 101.47$ ]

このおにぎり製造マシンの母平均は上記の区間の中に 95% の確率で収まることが分かった。

母分散が未知で大標本のとき

次に、母分散が未知で大標本のときについて。ここでいう大標本とは、標本のデータ数が概ね $n \ge 30$ や $n \ge 100$ のときを指す。

これも、母分散が既知のときと基本的な考え方は変わらない。ただし、母分散が分からないので、代わりに標本を元にした分散を使って推定する。

標本を元にして母分散に近づくように計算した分散を「不偏分散」という。使われる記号は $s^{2}$ だったり $u^{2}$ だったりと色々ある。ここでは $\hat{\sigma}^{2}$ を使うことにする。計算方法は次の通り。

$\hat{\sigma}^{2} = \frac{1}{n - 1} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}$

上記の数式の中で $n$ は標本のデータ数を表している。また、 $\bar{x}$ は標本平均を表している。

$\bar{x} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

あとは、不偏分散を元に計算した不偏標準偏差を母分散が既知なときの式に組み込むだけで良い。

$\bar{x} \pm z_{\frac{\alpha}{2}} \frac{\hat{\sigma}}{\sqrt{n}}$

上記の不偏標準偏差は単に不偏分散の平方根を取ったものだけど、一応示しておく。

$\hat{\sigma} = \sqrt{\frac{1}{n - 1} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}}$

例題

あるおにぎり製造マシンは作るおにぎりの重さを自由に調整できる。このとき、設定した重さと実際に出来上がる重さには誤差が出るが、その分布は正規分布に従うとする。この機械を使って試しにおにぎりを 49 個作ってみたところ、その平均は 100g で不偏分散は 16 だった。この機械が作るおにぎりの母平均を信頼係数 95% で区間推定せよ。

問題文から得られた値を、先ほどの数式に当てはめる。

$100 \pm 1.96 \frac{\sqrt{16}}{\sqrt{49}}$

上記を解くと、次の区間が得られる。

$[ 98.88, 101.12$ ]

このおにぎり製造マシンの母平均は上記の区間の中に 95% の確率で収まることが分かった。

別解

先ほどの例では標本を元に母分散に近づけた「不偏分散」を使うやり方だった。これとは別のやり方として「標本分散」を使うやり方もある。

不偏分散と標本分散の違いは割る数にある。不偏分散では $n - 1$ なのに対して標本分散では $n$ になる。ここでは標本分散を $s^{2}$ で表すことにする。

$s^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}$

標本標準偏差は、標本分散の平方根を取ったもの。

$s = \sqrt{\frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}}$

標本分散 (標準偏差) を使うときは縮み方が変わる。そこで式に組み込むときは $\sqrt{n}$ で割っていたところを $\sqrt{n - 1}$ にする。

$\bar{x} \pm z_{\frac{\alpha}{2}} \frac{s}{\sqrt{n - 1}}$

母分散が未知で小標本のとき

次は母分散が未知で小標本のとき。ここでいう小標本とは、標本のデータ数が概ね $n \le 30$ や $n \le 100$ のときを指す。

この状況も、先ほどの母分散が未知で大標本のときと基本的な考え方は変わらない。母分散が未知なので、不偏分散は引き続き使うことになる。ただし、標本が少ないときは分布が正規分布に従わなくなる。そこで、代わりに登場するのが自由度 f の t 分布になる。

t 分布というのは正規分布を少しつぶしたような形の分布をしている。さらに、自由度というパラメータによっても、その形が異なる。次のグラフは標準正規分布と共に自由度 1 と 2 の t 分布を描いている。 t 分布は自由度が大きくなるほど正規分布に近づいていく。

f:id:momijiame:20170103212624p:plain

ここで自由度というのは、どれだけの値が確定すれば全体の値が確定するかを示した値になっている。例えば標本が一次元の離散値なら、自由度はその標本として取り出したデータ点数から 1 を引いた値になる。

$f = n - 1$

これを数式に組み込むと、こうなる。内容的には、母分散が未知で大標本のときに使った正規分布の値を t 分布の値に入れ替えるだけ。

$\bar{x} \pm t_{\frac{\alpha}{2}}(n - 1) \frac{\hat{\sigma}}{\sqrt{n}}$

上記の $t_{\frac{\alpha}{2}}(n - 1)$ について自分で計算しないときは t 分布表を参照する。次の分布表では自由度 $f$ の t 分布を両側確率 $P$ で書いてある。

付表： t分布表 Student's t distribution — 中川雅央（滋賀大学）

これは、例えば標本の点数が 6 で信頼係数 95% なら $f = 5$ かつ $P = 0.05$ を参照する。つまり、このとき不等式は、次のようになる。

$\bar{x} \pm 2.570582 \frac{\hat{\sigma}}{\sqrt{n}}$

例題

あるおにぎり製造マシンは作るおにぎりの重さを自由に調整できる。このとき、設定した重さと実際に出来上がる重さには誤差が出るが、その分布は正規分布に従うとする。この機械を使って試しにおにぎりを 9 個作ってみたところ、その平均は 100g で不偏分散は 4 だった。この機械が作るおにぎりの母平均を信頼係数 95% で区間推定せよ。

まずは、標本のデータ点数 9 にもとづいて自由度 $f = 9 - 1$ となる。また、信頼係数 95% にもとづいて $P = 0.05$ となる。パラメータをもとに t 分布表を見ると 2.306004 を使えば良いことが分かった。

そして、問題文から得られた値を、先ほどの数式に当てはめる。

$100 \pm 2.306004 \frac{\sqrt{4}}{\sqrt{9}}$

上記を解くと、次の区間が得られる。

$[ 98.462664, 101.537336$ ]

このおにぎり製造マシンの母平均は上記の区間の中に 95% の確率で収まることが分かった。

別解

母分散が未知で大標本のときと同様に不偏分散の代わりに標本分散を使うこともできる。そのときは、こちらも同様に標本標準偏差を $n - 1$ で割る。

$\bar{x} \pm t_{\frac{\alpha}{2}}(n - 1) \frac{s}{\sqrt{n - 1}}$

標本分散の計算方法は先ほども書いたけど念のため。

$s^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}$

標本標準偏差は標本分散の平方根を取るだけ。

$s = \sqrt{\frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}}$

標本平均の計算方法は、こう。

$\bar{x} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

まとめ

今回は母平均の点・区間推定について書いた。点推定については、どういった状況でも標本平均が推定値になる。それに対して区間推定は状況ごとに使う手法を変える必要があることが分かった。

母分散が既知のとき
- 母分散と正規分布を使って推定する
母分散が未知で大標本のとき
- 不偏分散と正規分布を使って推定する
母分散が未知で小標本のとき
- 不偏分散と t 分布を使って推定する

参考文献

完全独習統計学入門

作者: 小島寛之
出版社/メーカー: ダイヤモンド社
発売日: 2006/09/28
メディア: 単行本（ソフトカバー）
購入: 215人クリック: 3,105回
この商品を含むブログ (115件) を見る

2017-01-02

統計: はじめての推定

統計

今回は、統計における重要な手法である「推定」について書いてみることにする。推定は、現実世界の様々な場面で使われている。例えば、選挙で開票作業が始まった直後に当選確実がニュースで流れることがある。一体どうしてそんなことが分かるのか、不思議に思ったことがあるかもしれない。実は、これには正に統計における推定が使われている。

ただ、推定のやり方は多種多様なので一つのエントリで書きつくすことは難しい。それをすると、むしろ分かりにくくなってしまうと思う。そこで、最初は推定の基本的な考え方と共に最もシンプルな推定について書くことにする。続きとなる別のパターンについては、おいおい書いていくと思う。

推定の説明に入る前に、いくつか前提となる知識について書いていく。

記述統計と推測統計

一口に統計学といっても、実は古典的な統計には大きく分けて二つのジャンルがある。それが記述統計と推測統計というもの。

まず、記述統計というのは、生のデータを要約することで分かりやすくすることを主に扱う分野になっている。これには例えば、データをグラフに加工することで見やすくしたり隠れた特徴を見出すことも含まれている。統計というと難しい数式を扱うもの、というイメージがあるかもしれないけど、それだけではないということ。

それに対して推測統計というのは、限られた部分集団から抽出元の情報を推測することを扱う分野になっている。例えば、先ほど例に挙げた選挙速報もこれに当たる。

選挙速報では、まず出口調査で誰に投票したかという情報を集める。この集めた情報は標本 (サンプル) と呼ばれる。標本は、その抽出元の特徴を受け継いだものになっている。ここでいう抽出元というのは有権者全員が誰に投票したか、という情報で母集団と呼ばれる。選挙速報では、出口調査で集めた標本から誰にどれだけ票が入ったかを推定する。

今回扱う推定は、後者の推測統計で扱われる。

どうやって推定するか？

推定のやり方には主に二つある。一つ目は点推定で、二つ目が区間推定というやり方だ。点推定というのは、最もそれらしい値を求めて「ここだ！」と一点賭けするやり方。それに対して区間推定というのは、推定する値はある確率でここからここまでの区間に収まるはずだ、と求めるやり方。

また、推定するものにも色々とある。例えば平均、分散、比率などが主なもの。先ほどの選挙速報の例では比率になる。

さらに、抽出元の情報がどこまで分かっているか、といった点でもやり方が変わってくる。例えば、選挙の例であれば母集団の比率を推定することになる。このとき、母分散は分かっているのか、母分布は分かっているのか、標本のサンプル数はどれだけあるのか。それぞれ状況によって使う手法は変わってくる。

最も単純な推定

今回は、考えうる最も単純な推定を扱ってみる。この場合、前提として母集団についての母数 (パラメータ) があらかじめ一通り分かっている状況になる。その状況で、標本について推定するのが最も単純な推定になると思う。

ここで、母集団というのは抽出元のデータのことをいう。それに対して標本というのは、母集団に含まれる値を取り出したものをいう。

母集団についての母数 (パラメータ) が分かっている状況というのは、次の通り。

平均について: 母平均が分かっている
分散について: 母分散が分かっている
分布について: 母分布が正規分布と分かっている

上記の言葉にあまり馴染みがなくても大丈夫、ここから解説していく。

平均について

平均と一口に言っても実は色々ある。ただし、ここでいう平均は最も一般的な算術平均を指している。

算術平均の定義は次の通り。ここで $\mu$ が平均になる。データ点数が $n$ で $x_i$ が各要素の値になる。ここで、各要素が出現する確率は等しいとする。

$\mu = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

数式は難しいけど、ようするにデータの総和をデータの個数で割ったもの。一例として、次のようなデータがあるとしよう。

$X = 10, 20, 30$

このデータについて平均を計算すると、次のようになる。

$\mu = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i = \frac{10 + 20 + 30}{3} = 20$

分散について

分散というのは、データがどれだけバラついているかを表す統計量の一つ。

分散の定義は次の通り。

$\sigma^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \mu)^{2}$

先ほどのデータについて分散も計算してみよう。

$X = 10, 20, 30$

ここで $\sigma^{2}$ が母分散になる。

$\sigma^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \mu)^{2} = \frac{(10 - 20)^{2} + (20 - 20)^{2} + (30 - 20)^{2}}{3} = 66.\dot{6}\dot{6}$

そして、この分散の平方根を取ると標準偏差という統計量になる。

$\sigma = \sqrt{\frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \mu)^{2}}$

先ほどのデータについて計算すると、こうなる。

$\sigma = \sqrt{\frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \mu)^{2}} = \sqrt{\frac{(10 - 20)^{2} + (20 - 20)^{2} + (30 - 20)^{2}}{3}} = 8.1649658$

この標準偏差という統計量は、後述する正規分布と深いつながりがある。

正規分布について

データの分布の仕方にも色々とある。その中でも、正規分布は自然界にもよくあるし統計でも扱いやすい性質を持っている。

まず、正規分布は平均を中心に釣鐘型の形になっている。

f:id:momijiame:20170102201116p:plain

ちなみに、上記は平均が 0 で標準偏差が 1 となっている。このような正規分布のことを標準正規分布と呼ぶ。

f:id:momijiame:20170102201125p:plain

標準正規分布は、どこまでの値が全体の何パーセントを占めているかが、あらかじめ分かっている。例えば、平均値 0 を中心として標準偏差±一つ分 ( $1 \sigma$ ) までなら全体の約 68% がそこに入る。同じように標準偏差±二つ分 (tex: 2 \sigma]) なら約 95% が入る。

f:id:momijiame:20170102201135p:plain

平均が 0 で標準偏差が 1 の標準正規分布に対して、それ以外の正規分布を一般正規分布と呼ぶ。次の一般正規分布は平均が 50 で標準偏差が 10 で、これは偏差値の分布になっている。

f:id:momijiame:20170102201145p:plain

一般正規分布は標準化という操作をすることで標準正規分布に変換できる。具体的には、一般正規分布の各要素の値から平均値を引いて標準偏差で割る。次の数式は、一般正規分布 $X$ を標準化して標準正規分布 $Z$ を得る過程を示している。

$Z = \frac{X - \mu}{\sigma}$

ここから言えることを考えてみよう。標準正規分布は、あらかじめどういった値がどれくらいの割合で含まれているかが分かっている。そして、一般正規分布は平均値と標準偏差を使って標準化することで標準正規分布に変換できる。つまり、データが正規分布だと知っていて平均値と標準偏差が分かっているとすれば、どういった値がどれくらいの割合で含まれるかが計算できるということ。

テストを例に考えてみる

例えば、あるテストで A さんは 90 点を取ったとしよう。このテストを受けた全員の結果をまとめると正規分布になっていた。そして、平均点が 60 で標準偏差が 15 だったとする。この情報にもとづいて A さんの取った点数を標準化してみることにする。

$\frac{90 - 60}{15} = 2$

結果は 2 となった。次は、この結果がどういった意味を持つのかについて説明する。以下のサイトにある標準正規分布表を見てもらいたい。

標準正規分布表

この標準正規分布表には、標準正規分布において平均からの距離で、それ以上の値が全体のどれだけの割合で含まれるか書かれている。表の見方は、まず縦の行が小数点第一位まで、そして横の列がそれに対応する小数点第二位になっている。表のヘッダ部分に書かれている数字は平均からの距離を表している。そして、表の中身に書かれている数字は、その距離よりも大きな値が全体のどれだけ含まれているかになる。例えば一番左上は 0.0 で平均ちょうどを表している。その中身として書かれている 0.5 という数字は 0.0 よりも大きな値は全体の 50% ある、という意味になる。

表の見方が分かったところで $2 \sigma$ つまり 2.0 のところ、を見てみよう。中身として書かれている値は 0.02275 つまり約 2% になる。つまり A さんの成績は上位 2% に位置していることが分かる。

また、仮に上記の操作を全員分のテストの結果に適用する場合を考えてみよう。すると、元が平均 60 で標準偏差が 15 の一般正規分布が、平均 0 で標準偏差が 1 の標準正規分布になるというわけ。

点推定

だいぶ回り道になってしまったけど、これで前提知識が一通り説明できた。なので、話を推定に戻すことにする。まず、今回はデータ (母集団) が正規分布で、平均値と分散 (標準偏差) が既に分かっている状況を扱うという話だった。

上記のようなシチュエーションで、この母集団から標本を一つ取り出すとしたら、どういった値が出てくるか予測したい。このとき、推定の一つ目の手法である点推定では、最も出やすい値を一つだけ選ぶやり方になる。そして、正規分布において最も出やすい値というのは分布の中央にある平均値を指す。

$\mu = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

つまり、正規分布から点推定をするなら平均値を答えれば良いということ。しかし、理論上最も当たりやすいとはいえ、これが実際に当たるかといえば、そうは思えないはず。

区間推定

そこで登場するのが区間推定という考え方になる。区間推定では、その名の通り「ここからここまで」という区間で推定する。

とはいえ、どの区間を推定で使うかというのが重要になる。当たり前だけど、区間に $[ -\infty, \infty$ ] を指定すれば、その推定は 100% 当たる。しかし、そんな推定は何の意味もない。

そこで、区間推定では信頼係数という概念を用いる。この信頼係数には一般的に 95% や 99% という値を使う。つまり、信頼係数の確率で推定した値がその区間に入りそうだ、ということ。

また、信頼係数は別の表現方法として、次のように有意水準 $\alpha$ を使うこともある。

$100(1 - \alpha) \%$

有意水準というのは、言い換えれば推定した値が求めた区間の外にある確率でもある。もちろん 100% ではない以上、区間内に推定した値が入らない恐れはあるということ。信頼係数 95% というのは、有意水準 5% でその区間推定が間違っている確率でもある。

標本の区間推定

区間推定の考え方の説明が終わったので、次は実際に計算してみよう。信頼係数には 95% (つまり有意水準 5%) を採用する。

区間推定では、なるべく狭い区間で信頼係数を満たす確率を達成したい。そのため、正規分布の区間推定では、平均値を中心として左右に同じ距離を指定する。ちょうど、先に出したグラフに信頼係数 95% の区間推定の範囲があった。ただし、グラフでは $\pm 2 \sigma$ と書いていたけど、実際には $\pm 1.96 \sigma$ を使う

f:id:momijiame:20170102201135p:plain

この 1.96 という数字が何処から出てきたかというと、もう一度標準正規分布表を見ると分かる。

標準正規分布表

ここで 1.96 のところを見ると 0.024998 と、ほぼ 2.5% になっていることが分かる。正規分布では平均値を中心に大きい値と小さな値が均等に出るので、有意水準 5% を左右に均等に割り振る。そのため、この値になる。もし、信頼係数が 99% なら、有意水準 1% を左右に均等に割り振って 0.5% になる。そのときは代わりに 2.58 を使う。

つまり、標準正規分布で信頼係数 95% の区間推定をするなら平均値 0 を中心に $\pm 1.96$ を指定すれば良いことが分かる。これを不等式で表現してみよう。この式を満たす $x$ が、標準正規分布から標本をランダムに取り出したとき 95% で当たる値になる。

$-1.96 \le x \le 1.96$

余談だけど、これを信頼係数を決めずに書くときは、次のような式で表現する。

$-z_{\frac{\alpha}{2}} \le x \le z_{\frac{\alpha}{2}}$

ただし、上記は標準正規分布の場合、ということに注意が必要となる。巷に存在する正規分布は一般正規分布だからだ。

少し前に説明した内容を思い出すと、一般正規分布は標準正規分布に変換できた。それは標準化という操作で、平均値を引いて標準偏差で割る、というものだった。これを、先ほどの不等式に組み込んでみよう。

$-1.96 \le \frac{x - \mu}{\sigma} \le 1.96$

これで標準正規分布の区間推定を一般正規分布に拡張できた。

この式を $x$ について解いてみよう。まずは三辺に $\sigma$ を掛ける。

$-1.96 \sigma \le x - \mu \le 1.96 \sigma$

次に三辺に $\mu$ を足す。

$\mu - 1.96 \sigma \le x \le \mu + 1.96 \sigma$

つまり、母数 (パラメータ) が既知な一般正規分布から標本をランダムに一つ取り出したときの信頼係数 95% の区間推定は、この値になる。

$\mu \pm 1.96 \sigma$

例題

例えば、あるプロ野球チームの平均身長が 182cm で、標準偏差が 9cm だったとする。この分布が正規分布に従うとして、一人の選手をランダムに抽出することを考える。抽出した選手の身長を信頼係数 95% で区間推定しよう。

問題文から得られた、分布の母数 (パラメータ) を先ほどの式に代入する。

$182 \pm 1.96 \times 9$

上記を解くと、次のような区間が得られる。

$[ 164.36, 199.64$ ]

問題文にある母集団から一人の選手をランダムに抽出すると 95% の確率で上記の範囲に収まることが分かった。

標本平均の区間推定

先ほどまでの内容は、母集団から一つの標本を取り出すときのことを考えていた。次は母集団からいくつかの標本を取り出して平均を取ったときのことを考えてみよう。このいくつかの標本を元にした平均を標本平均と呼ぶ。記号では $\bar{x}$ で表されることが多い。

考え方は一つの標本の区間推定と変わらない。ただし、全く同じではない。なぜなら、標本平均の分布は、母分布に比べると標準偏差が縮むため。

ここでいう標本平均の分布というのは、母集団から何度も何度も無限に標本平均を計算してプロットしたもの。このとき、平均値は母分布と変わらないものの標準偏差は変わってしまう。具体的には母分布の標準偏差が $\sigma$ のとき標本平均の分布では $\frac{\sigma}{\sqrt{n}}$ となる。つまり、標本平均を計算するのに使ったデータ点数に応じて標準偏差が縮むということ。

これを標準化するときの公式にも反映する必要がある。そのため、先ほどの公式は次のように変わる。推定するものが標本平均 $\bar{x}$ になったのと、標準偏差が $\frac{\sigma}{\sqrt{n}}$ と小さくなる。

$-1.96 \le \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}} \le 1.96$

上記を標本平均 $\bar{x}$ について解くと、こうなる。

$\mu -1.96 \frac{\sigma}{\sqrt{n}} \le \bar{x} \le \mu + 1.96 \frac{\sigma}{\sqrt{n}}$

つまり、標本平均の信頼係数 95% の区間推定はこう。

$\mu \pm 1.96 \frac{\sigma}{\sqrt{n}}$

例題

例えば、あるプロ野球チームの平均身長が 182cm で、標準偏差が 9cm だったとする。この分布が正規分布に従うとして、四人の選手をランダムに抽出することを考える。抽出した選手の身長を信頼係数 95% で区間推定しよう。

問題文から得られた、分布の母数 (パラメータ) を先ほどの式に代入する。

$182 \pm 1.96 \frac{9}{\sqrt{4}}$

これを計算すると、こう。

$[ 173.18, 190.82$ ]

標本平均の区間推定は、先ほどの一つの標本の区間推定よりも幅が狭まっている。これは標本平均に使ったデータ点数だけ標準偏差が縮んでいるため。

まとめ

ということで今回は母集団の母数 (パラメータ) が既知なときの標本について推定をしてみた。ここまでに扱ってきた内容は次の通り。

統計学には二つのジャンルがある
- 記述統計
- 推測統計
推定には二つのやり方がある
- 点推定
- 区間推定
推定する対象や既知な内容によって手法が枝分かれする
- 推定する対象: 平均、分散、比率など
- 既知になりうるもの: 分布、平均、分散など
点推定の手順
- 既知な平均値を使う
区間推定の手順
- 信頼係数を決める
- 一般正規分布を標準正規分布に標準化する
- 不等式にする
- 標本について解く
標本平均の性質
- 標準偏差が $\frac{\sigma}{\sqrt{n}}$ に縮む

参考文献

完全独習統計学入門

作者: 小島寛之
出版社/メーカー: ダイヤモンド社
発売日: 2006/09/28
メディア: 単行本（ソフトカバー）
購入: 215人クリック: 3,105回
この商品を含むブログ (115件) を見る

2016-12-23

統計: Python と R で重回帰分析してみる

統計 scikit-learn Python NumPy Mac OS X

今回は R と Python の両方を使って重回帰分析をしてみる。モチベーションとしては、できるだけ手に慣れた Python を使って分析をしていきたいという気持ちがある。ただ、計算結果が意図通りのものになっているのかを R の結果と見比べて確かめておきたい。

また、分析にはボストンデータセットを用いる。このデータセットはボストンの各地区ごとの不動産の平均価格と、それに付随する情報が入っている。

今回使った環境は次の通り。 R と Python は、あらかじめインストール済みであることを想定している。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1212
$ python --version
Python 3.5.2
$ R --version
R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under the terms of the
GNU General Public License versions 2 or 3.
For more information about these matters see
http://www.gnu.org/licenses/.

もし、まだ R をインストールしていないのであれば、こちらの記事が参考になるかも。

blog.amedama.jp

R

まずは R のやり方から。

最初に R を起動する。

$ R

続いてボストンデータセットの入っているライブラリ MASS を読み込む。

> library(MASS)

データセットの中身は次のようになっている。

> head(Boston)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
  medv
1 24.0
2 21.6
3 34.7
4 33.4
5 36.2
6 28.7

各変数は次のような意味を持っている。ただ、今回はボストンデータセットについて詳しく見ることが目的ではないので詳しくは立ち入らない。

1. CRIM: per capita crime rate by town
2. ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
3. INDUS: proportion of non-retail business acres per town
4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
5. NOX: nitric oxides concentration (parts per 10 million)
6. RM: average number of rooms per dwelling
7. AGE: proportion of owner-occupied units built prior to 1940
8. DIS: weighted distances to five Boston employment centres
9. RAD: index of accessibility to radial highways
10. TAX: full-value property-tax rate per $10,000
11. PTRATIO: pupil-teacher ratio by town
12. B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
13. LSTAT: % lower status of the population
14. MEDV: Median value of owner-occupied homes in $1000's
(https://archive.ics.uci.edu/ml/datasets/Housing より)

今回はシンプルに線形回帰モデルを使うので lm() 関数を使う。 medv は目的変数である「その地区の不動産の平均価格」になっている。次の処理では、目的関数をその他の全ての変数を使って重回帰している。

> lmfit = lm(medv ~ ., data=Boston)

重回帰した結果の詳細は summary() 関数で取得できる。

> summary(lmfit)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max
-15.595  -2.730  -0.518   1.777  26.199

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 **
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288
chas         2.687e+00  8.616e-01   3.118 0.001925 **
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 **
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

自由度調整済み決定係数が 0.7338 ということで、あまり上手く説明できていなさそうな感じ。

Python

続いて Python でのやり方について。今回は statsmodels というライブラリを使うことにする。ちなみに、線形回帰モデルは別のライブラリにも色々な実装がある。ただ、フィッティングしたときに得られるパラメータが statsmodels は細かく取れて R に近いので、今回はこちらを使ってみた。

まずは必要なライブラリをインストールする。今回は回帰に使わない scikit-learn を入れているのはボストンデータセットを読み込むため。

$ pip install statsmodels scikit-learn

最初に Python の REPL を起動しよう。

$ python

起動したら scikit-learn からボストンデータセットを読み込む。

>>> from sklearn import datasets
>>> dataset = datasets.load_boston()

このデータセットには次のようなメンバがある。

>>> dir(dataset)
['DESCR', 'data', 'feature_names', 'target']

target メンバが目的変数となる各地区ごとの不動産物件の平均価格となっている。

>>> dataset.target[:10]
array([ 24. ,  21.6,  34.7,  33.4,  36.2,  28.7,  22.9,  27.1,  16.5,  18.9])

そして説明変数が data に入っている。ボストンを 506 地区に分割して 13 の変数を与えたデータになっている。

>>> dataset.data.shape
(506, 13)

変数名は feature_names から取得できる。

>>> dataset.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT'],
      dtype='<U7')

データセットを読み込めたところで statsmodels をインポートしよう。

>>> from statsmodels import api as sm

線形回帰モデルを使うときは OLS クラスを使う。第一引数に目的変数を、第二変数に説明変数を渡す。

>>> model = sm.OLS(dataset.target, dataset.data)

準備ができたら fit() メソッドを使って分析する。

>>> result = model.fit()

分析結果からは summary() メソッドを使って R の summary() 関数に近い内容が取得できる。

>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.959
Model:                            OLS   Adj. R-squared:                  0.958
Method:                 Least Squares   F-statistic:                     891.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):               0.00
Time:                        22:59:04   Log-Likelihood:                -1523.8
No. Observations:                 506   AIC:                             3074.
Df Residuals:                     493   BIC:                             3129.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1            -0.0916      0.034     -2.675      0.008        -0.159    -0.024
x2             0.0487      0.014      3.379      0.001         0.020     0.077
x3            -0.0038      0.064     -0.059      0.953        -0.130     0.123
x4             2.8564      0.904      3.160      0.002         1.080     4.633
x5            -2.8808      3.359     -0.858      0.392        -9.481     3.720
x6             5.9252      0.309     19.168      0.000         5.318     6.533
x7            -0.0072      0.014     -0.523      0.601        -0.034     0.020
x8            -0.9680      0.196     -4.947      0.000        -1.352    -0.584
x9             0.1704      0.067      2.554      0.011         0.039     0.302
x10           -0.0094      0.004     -2.393      0.017        -0.017    -0.002
x11           -0.3924      0.110     -3.571      0.000        -0.608    -0.177
x12            0.0150      0.003      5.561      0.000         0.010     0.020
x13           -0.4170      0.051     -8.214      0.000        -0.517    -0.317
==============================================================================
Omnibus:                      204.050   Durbin-Watson:                   0.999
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1372.527
Skew:                           1.609   Prob(JB):                    9.11e-299
Kurtosis:                      10.399   Cond. No.                     8.50e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.5e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

めでたしめでたし、と行きたいところだけどちょっと待ってほしい。先ほど取得した R の結果と内容が異なっている。

実は R の分析では説明変数にバイアス項という変数が自動で挿入されていた。これは、不動産物件のベースとなる価格と捉えることができる。

statsmodels を使っているときにバイアス項を説明変数に挿入するときは statsmodels.add_constant() 関数を使う。バイアス項を挿入した説明変数を X という名前で保存しておこう。ついでに目的変数も Y という変数で保存しておく。

>>> X = sm.add_constant(dataset.data)
>>> Y = dataset.target

バイアス項を追加した説明変数を用いて、先ほどと同じように回帰してみよう。

>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          6.95e-135
Time:                        23:00:34   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         36.4911      5.104      7.149      0.000        26.462    46.520
x1            -0.1072      0.033     -3.276      0.001        -0.171    -0.043
x2             0.0464      0.014      3.380      0.001         0.019     0.073
x3             0.0209      0.061      0.339      0.735        -0.100     0.142
x4             2.6886      0.862      3.120      0.002         0.996     4.381
x5           -17.7958      3.821     -4.658      0.000       -25.302   -10.289
x6             3.8048      0.418      9.102      0.000         2.983     4.626
x7             0.0008      0.013      0.057      0.955        -0.025     0.027
x8            -1.4758      0.199     -7.398      0.000        -1.868    -1.084
x9             0.3057      0.066      4.608      0.000         0.175     0.436
x10           -0.0123      0.004     -3.278      0.001        -0.020    -0.005
x11           -0.9535      0.131     -7.287      0.000        -1.211    -0.696
x12            0.0094      0.003      3.500      0.001         0.004     0.015
x13           -0.5255      0.051    -10.366      0.000        -0.625    -0.426
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                     1.51e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.51e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

今度は R で得られた内容と一致している！

AIC

重回帰分析では、モデリングに用いる説明変数の選択基準として AIC という数値を小さくするやり方があるので、その計算方法について。

R

R では AIC() 関数が用意されているので、それに回帰結果を渡す。

> AIC(lmfit)
[1] 3027.609

Python

Python では statsmodels の回帰結果のインスタンスに aic というメンバがあって、そこに格納されている。ちなみに AIC はサマリーの中にも表示されていた。

>>> result.aic
3025.6767200074596

標準化

標準化というのはデータセットを平均が 0 で標準偏差が 1 になるように加工すること。説明変数が標準化されていると、それぞれの説明変数が目的変数に対してどれくらい影響を与えているかを見たりするのに分かりやすい。

R

R には標準化のための関数として scale() が用意されている。

> X <- scale(Boston[, -14])
> Y <- Boston[, 14]

説明変数を標準化して回帰してみよう。

> lmfit <- lm(Y ~ X)
> summary(lmfit)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max
-15.595  -2.730  -0.518   1.777  26.199

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.53281    0.21095 106.814  < 2e-16 ***
Xcrim       -0.92906    0.28269  -3.287 0.001087 **
Xzn          1.08264    0.32016   3.382 0.000778 ***
Xindus       0.14104    0.42188   0.334 0.738288
Xchas        0.68241    0.21884   3.118 0.001925 **
Xnox        -2.05875    0.44262  -4.651 4.25e-06 ***
Xrm          2.67688    0.29364   9.116  < 2e-16 ***
Xage         0.01949    0.37184   0.052 0.958229
Xdis        -3.10712    0.41999  -7.398 6.01e-13 ***
Xrad         2.66485    0.57770   4.613 5.07e-06 ***
Xtax        -2.07884    0.63379  -3.280 0.001112 **
Xptratio    -2.06265    0.28323  -7.283 1.31e-12 ***
Xblack       0.85011    0.24521   3.467 0.000573 ***
Xlstat      -3.74733    0.36216 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

説明変数が標準化されたことで、どの変数がどれくらい影響を与えているのかが分かりやすくなった。例えば lstat などは影響が大きそうだ。

Python

Python の場合は scikit-learn の preprocessing() 関数を使うと簡単に標準化できる。

>>> X = dataset.data
>>> from sklearn import preprocessing
>>> X_scaled = preprocessing.scale(X)

標準化した上でバイアス項を追加して回帰してみよう。

>>> X_scaled_with_constant = sm.add_constant(X_scaled)
>>> model = sm.OLS(Y, X_scaled_with_constant)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          6.95e-135
Time:                        23:50:36   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         22.5328      0.211    106.807      0.000        22.118    22.947
x1            -0.9204      0.281     -3.276      0.001        -1.472    -0.368
x2             1.0810      0.320      3.380      0.001         0.453     1.709
x3             0.1430      0.421      0.339      0.735        -0.685     0.971
x4             0.6822      0.219      3.120      0.002         0.253     1.112
x5            -2.0601      0.442     -4.658      0.000        -2.929    -1.191
x6             2.6706      0.293      9.102      0.000         2.094     3.247
x7             0.0211      0.371      0.057      0.955        -0.709     0.751
x8            -3.1044      0.420     -7.398      0.000        -3.929    -2.280
x9             2.6588      0.577      4.608      0.000         1.525     3.792
x10           -2.0759      0.633     -3.278      0.001        -3.320    -0.832
x11           -2.0622      0.283     -7.287      0.000        -2.618    -1.506
x12            0.8566      0.245      3.500      0.001         0.376     1.338
x13           -3.7487      0.362    -10.366      0.000        -4.459    -3.038
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                         9.82
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
"""

結果が R の内容と一致していることが分かる。

ちなみに標準化の操作自体は単純なので、自前でやることもできる。数式に直すと、こんな感じ。各要素から平均 ( $\mu$ ) を引いて標準偏差 ( $\sigma$ ) で割る。

$Z = \frac{X - \mu}{\sigma}$

コードにするとこんな感じ。

>>> X_scaled = (X - X.mean(axis=0)) / X.std(axis=0)

交互作用モデル

交互作用モデルというのは、複数の変数をまとめて説明変数として扱うやり方。ある説明変数と別の説明変数の相乗作用があるときに使うと良い。

R

R では lm() 関数を実行するときに説明変数を ^2 するだけで交互作用モデルを使える。

> lmfit = lm(medv ~ .^2, data=Boston)
> summary(lmfit)

Call:
lm(formula = medv ~ .^2, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max
-7.9374 -1.5344 -0.1068  1.2973 17.8500

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)   -1.579e+02  6.800e+01  -2.323 0.020683 *
crim          -1.707e+01  6.554e+00  -2.605 0.009526 **
zn            -7.529e-02  4.580e-01  -0.164 0.869508
indus         -2.819e+00  1.696e+00  -1.663 0.097111 .
chas           4.451e+01  1.952e+01   2.280 0.023123 *
nox            2.006e+01  7.516e+01   0.267 0.789717
rm             2.527e+01  5.699e+00   4.435 1.18e-05 ***
age            1.263e+00  2.728e-01   4.630 4.90e-06 ***
dis           -1.698e+00  4.604e+00  -0.369 0.712395
rad            1.861e+00  2.464e+00   0.755 0.450532
tax            3.670e-02  1.440e-01   0.255 0.798978
ptratio        2.725e+00  2.850e+00   0.956 0.339567
black          9.942e-02  7.468e-02   1.331 0.183833
lstat          1.656e+00  8.533e-01   1.940 0.053032 .
crim:zn        4.144e-01  1.804e-01   2.297 0.022128 *
crim:indus    -4.693e-02  4.480e-01  -0.105 0.916621
crim:chas      2.428e+00  5.710e-01   4.251 2.63e-05 ***
crim:nox      -1.108e+00  9.285e-01  -1.193 0.233425
crim:rm        2.163e-01  4.907e-02   4.409 1.33e-05 ***
crim:age      -3.083e-03  3.781e-03  -0.815 0.415315
crim:dis      -1.903e-01  1.060e-01  -1.795 0.073307 .
crim:rad      -6.584e-01  5.815e-01  -1.132 0.258198
crim:tax       3.479e-02  4.287e-02   0.812 0.417453
crim:ptratio   4.915e-01  3.328e-01   1.477 0.140476
crim:black    -4.612e-04  1.793e-04  -2.572 0.010451 *
crim:lstat     2.964e-02  6.544e-03   4.530 7.72e-06 ***
...(省略)...
rad:tax        3.131e-05  1.446e-03   0.022 0.982729
rad:ptratio   -4.379e-02  8.392e-02  -0.522 0.602121
rad:black     -4.362e-04  2.518e-03  -0.173 0.862561
rad:lstat     -2.529e-02  1.816e-02  -1.392 0.164530
tax:ptratio    7.854e-03  2.504e-03   3.137 0.001830 **
tax:black     -4.785e-07  1.999e-04  -0.002 0.998091
tax:lstat     -1.403e-03  1.208e-03  -1.162 0.245940
ptratio:black  1.203e-03  3.361e-03   0.358 0.720508
ptratio:lstat  3.901e-03  2.985e-02   0.131 0.896068
black:lstat   -6.118e-04  4.157e-04  -1.472 0.141837
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.852 on 414 degrees of freedom
Multiple R-squared:  0.9212, Adjusted R-squared:  0.9039
F-statistic: 53.18 on 91 and 414 DF,  p-value: < 2.2e-16

結果から説明変数の数が一気に増えてることが分かる。これは 13 種類の説明変数から 2 種類を組み合わせで取り出して、それらを交互作用モデルの説明変数として追加しているため。

また、自由度調整済み決定係数が 0.9039 となっており、ただ 13 種類の変数を使うときよりも上手くモデリングできていることが分かる。まあ、この値だけを見ていても過適合の恐れがあるけど。

Python

statsmodels には回帰するときのやり方を R のような書き方 (R-style らしい) で指定できる。このやり方を使うとさくっと相互作用モデルも使えそうな感じ。

Fitting models using R-style formulas — statsmodels 0.7.0 documentation

ただまあ、これを使うよりも自前で用意する方が楽しそう。なので、今回は自分で計算してみることにする。

Python で組み合わせを計算するときは itertools モジュールを使うと良い。この中にある combinations() 関数を使えば簡単に組み合わせが取れる。例えば説明変数の組み合わせを計算したいなら、こう。

>>> import itertools
>>> list(itertools.combinations(dataset.feature_names, 2))
[('CRIM', 'ZN'), ('CRIM', 'INDUS'), ('CRIM', 'CHAS'), ('CRIM', 'NOX'), ('CRIM', 'RM'), ('CRIM', 'AGE'), ('CRIM', 'DIS'), ('CRIM', 'RAD'), ('CRIM', 'TAX'), ('CRIM', 'PTRATIO'), ('CRIM', 'B'), ('CRIM', 'LSTAT'), ('ZN', 'INDUS'), ('ZN', 'CHAS'), ('ZN', 'NOX'), ('ZN', 'RM'), ('ZN', 'AGE'), ('ZN', 'DIS'), ('ZN', 'RAD'), ('ZN', 'TAX'), ('ZN', 'PTRATIO'), ('ZN', 'B'), ('ZN', 'LSTAT'), ('INDUS', 'CHAS'), ('INDUS', 'NOX'), ('INDUS', 'RM'), ('INDUS', 'AGE'), ('INDUS', 'DIS'), ('INDUS', 'RAD'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('INDUS', 'B'), ('INDUS', 'LSTAT'), ('CHAS', 'NOX'), ('CHAS', 'RM'), ('CHAS', 'AGE'), ('CHAS', 'DIS'), ('CHAS', 'RAD'), ('CHAS', 'TAX'), ('CHAS', 'PTRATIO'), ('CHAS', 'B'), ('CHAS', 'LSTAT'), ('NOX', 'RM'), ('NOX', 'AGE'), ('NOX', 'DIS'), ('NOX', 'RAD'), ('NOX', 'TAX'), ('NOX', 'PTRATIO'), ('NOX', 'B'), ('NOX', 'LSTAT'), ('RM', 'AGE'), ('RM', 'DIS'), ('RM', 'RAD'), ('RM', 'TAX'), ('RM', 'PTRATIO'), ('RM', 'B'), ('RM', 'LSTAT'), ('AGE', 'DIS'), ('AGE', 'RAD'), ('AGE', 'TAX'), ('AGE', 'PTRATIO'), ('AGE', 'B'), ('AGE', 'LSTAT'), ('DIS', 'RAD'), ('DIS', 'TAX'), ('DIS', 'PTRATIO'), ('DIS', 'B'), ('DIS', 'LSTAT'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('RAD', 'B'), ('RAD', 'LSTAT'), ('TAX', 'PTRATIO'), ('TAX', 'B'), ('TAX', 'LSTAT'), ('PTRATIO', 'B'), ('PTRATIO', 'LSTAT'), ('B', 'LSTAT')]

組み合わせは ${}_{13} C_2 = \frac{13!}{2!(13 - 2)!} = 78$ 通り。

>>> len(list(itertools.combinations(dataset.feature_names, 2)))
78

それぞれの説明変数を取り出すための添字にするとこんな感じ。

>>> list(itertools.combinations(range(dataset.data.shape[1]), 2))
[(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10), (0, 11), (0, 12), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10), (1, 11), (1, 12), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10), (2, 11), (2, 12), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10), (3, 11), (3, 12), (4, 5), (4, 6), (4, 7), (4, 8), (4, 9), (4, 10), (4, 11), (4, 12), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10), (5, 11), (5, 12), (6, 7), (6, 8), (6, 9), (6, 10), (6, 11), (6, 12), (7, 8), (7, 9), (7, 10), (7, 11), (7, 12), (8, 9), (8, 10), (8, 11), (8, 12), (9, 10), (9, 11), (9, 12), (10, 11), (10, 12), (11, 12)]

交互作用モデルで使うそれぞれの説明変数は、説明変数同士を乗算することで計算できる。

>>> b = np.array([X[:, 0] * X[:, 1]])

計算した説明変数を、既存の説明変数にカラムとして追加するには次のようにすれば良い。

>>> a = X
>>> np.concatenate((a, b.T), axis=1)
array([[  6.32000000e-03,   1.80000000e+01,   2.31000000e+00, ...,
          3.96900000e+02,   4.98000000e+00,   1.13760000e-01],
       [  2.73100000e-02,   0.00000000e+00,   7.07000000e+00, ...,
          3.96900000e+02,   9.14000000e+00,   0.00000000e+00],
       [  2.72900000e-02,   0.00000000e+00,   7.07000000e+00, ...,
          3.92830000e+02,   4.03000000e+00,   0.00000000e+00],
       ...,
       [  6.07600000e-02,   0.00000000e+00,   1.19300000e+01, ...,
          3.96900000e+02,   5.64000000e+00,   0.00000000e+00],
       [  1.09590000e-01,   0.00000000e+00,   1.19300000e+01, ...,
          3.93450000e+02,   6.48000000e+00,   0.00000000e+00],
       [  4.74100000e-02,   0.00000000e+00,   1.19300000e+01, ...,
          3.96900000e+02,   7.88000000e+00,   0.00000000e+00]])

これで説明変数が一つ増える。

>>> np.concatenate((a, b.T), axis=1).shape
(506, 14)

上記にもとづいて相互作用モデルで使う全ての変数を揃える。

>>> X = dataset.data
>>> for i, j in itertools.combinations(range(X.shape[1]), 2):
...     interaction_column = np.array([X[:, i] * X[:, j]])
...     X = np.concatenate((X, interaction_column.T), axis=1)
...
>>> X.shape
(506, 91)

そして、バイアス項を追加する。

>>> X = sm.add_constant(X)
>>> X.shape
(506, 92)

上記の説明変数を使って回帰する。

>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.921
Model:                            OLS   Adj. R-squared:                  0.904
Method:                 Least Squares   F-statistic:                     53.21
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          5.85e-181
Time:                        19:16:55   Log-Likelihood:                -1197.3
No. Observations:                 506   AIC:                             2579.
Df Residuals:                     414   BIC:                             2967.
Df Model:                          91
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const       -159.4473     67.980     -2.346      0.019      -293.076   -25.819
x1           -16.9621      6.549     -2.590      0.010       -29.835    -4.089
x2            -0.0769      0.458     -0.168      0.867        -0.977     0.823
x3            -2.8163      1.695     -1.662      0.097        -6.148     0.516
x4            44.6397     19.519      2.287      0.023         6.271    83.008
x5            22.1561     75.150      0.295      0.768      -125.567   169.879
x6            25.4191      5.700      4.459      0.000        14.214    36.624
x7             1.2569      0.273      4.612      0.000         0.721     1.793
x8            -1.6398      4.604     -0.356      0.722       -10.690     7.410
x9             1.8056      2.462      0.733      0.464        -3.034     6.646
x10            0.0374      0.144      0.260      0.795        -0.246     0.320
x11            2.7520      2.849      0.966      0.335        -2.849     8.353
x12            0.1008      0.074      1.354      0.177        -0.046     0.247
x13            1.6587      0.853      1.944      0.053        -0.018     3.336
x14            0.4124      0.180      2.287      0.023         0.058     0.767
x15           -0.0467      0.448     -0.104      0.917        -0.927     0.834
x16            2.4377      0.571      4.271      0.000         1.316     3.560
x17           -1.2206      0.920     -1.327      0.185        -3.029     0.588
x18            0.2117      0.049      4.345      0.000         0.116     0.307
x19           -0.0031      0.004     -0.828      0.408        -0.011     0.004
x20           -0.1956      0.106     -1.848      0.065        -0.404     0.012
x21           -0.6599      0.581     -1.135      0.257        -1.803     0.483
x22            0.0349      0.043      0.815      0.416        -0.049     0.119
x23            0.4895      0.333      1.471      0.142        -0.164     1.143
x24           -0.0004      0.000     -2.513      0.012        -0.001 -9.71e-05
x25            0.0295      0.007      4.512      0.000         0.017     0.042
...(省略)...
x82         3.151e-05      0.001      0.022      0.983        -0.003     0.003
x83           -0.0439      0.084     -0.523      0.601        -0.209     0.121
x84           -0.0004      0.003     -0.165      0.869        -0.005     0.005
x85           -0.0251      0.018     -1.384      0.167        -0.061     0.011
x86            0.0078      0.003      3.135      0.002         0.003     0.013
x87        -2.303e-06      0.000     -0.012      0.991        -0.000     0.000
x88           -0.0014      0.001     -1.167      0.244        -0.004     0.001
x89            0.0012      0.003      0.349      0.727        -0.005     0.008
x90            0.0038      0.030      0.128      0.898        -0.055     0.062
x91           -0.0006      0.000     -1.510      0.132        -0.001     0.000
==============================================================================
Omnibus:                      121.344   Durbin-Watson:                   1.524
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              628.079
Skew:                           0.942   Prob(JB):                    4.11e-137
Kurtosis:                       8.123   Cond. No.                     1.18e+08
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.18e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

上記を見ると R で相互作用モデルを使って計算したときの内容と一致していることが分かる。

まとめ

今回は R と Python (statsmodels) で線形モデルを使った重回帰分析をしてみた。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

作者: もみじあめ
発売日: 2016/03/12
メディア: Kindle版
この商品を含むブログを見る

2016-12-23

統計: 二つのグループの平均と分散を合成する

Mac OS X Python 統計

例えば、あるグループ A と B が別々にテストを受けたとする。それぞれのグループの人数と平均点、そして分散は分かっているとしよう。このとき、グループ A と B を合わせた全体での平均や分散は計算することができるだろうか？

結論から言うと、これはできる。今回は、その手順について書いてみることにする。また、同時に Python を使ってサンプルデータを生成したり計算の裏付けをしてみよう。

使った環境は次の通り。

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

題材とするデータを生成する

平均や分散の計算を楽にするために、まずは NumPy をインストールしておこう。

$ pip install numpy

そして Python の REPL を起動する。

$ python

起動したら NumPy をインポートしよう。

>>> import numpy as np

そして、まずはグループ A を生成する。パラメータは適当だけど、点数は 40 から 90 の範囲で 10 人分作っておく。

>>> A = np.random.randint(40, 90, 10)
>>> A
array([57, 67, 74, 79, 82, 53, 80, 46, 74, 47])

上記の平均と分散はそれぞれ 65.9 と 176.09 になった。

>>> A.mean()
65.900000000000006
>>> A.var()
176.09

同じようにグループ B も生成する。こちらは平均が 72.466... で分散が 60.9155... になった。

>>> B = np.random.randint(55, 85, 15)
>>> B
array([83, 77, 77, 59, 67, 57, 77, 71, 68, 77, 72, 76, 84, 64, 78])
>>> B.mean()
72.466666666666669
>>> B.var()
60.915555555555549

そして、上記を合わせたグループについても計算しておこう。便宜的に、ここではこれをグループ C と呼ぶことにする。

>>> C = np.r_[A, B]
>>> C
array([57, 67, 74, 79, 82, 53, 80, 46, 74, 47, 83, 77, 77, 59, 67, 57, 77,
       71, 68, 77, 72, 76, 84, 64, 78])

グループ C の平均と分散は次のようになった。今回は、この値をグループ A と B の平均と分散とデータ点数からだけで求めるのが最終的な目的となる。

>>> C.mean()
69.840000000000003
>>> C.var()
117.3344

平均を合成する

平均の合成はすごく簡単。まず、それぞれのグループの平均値に人数をかけると各グループの点数の合計が得られる。あとは、それをまとめてから全体の人数で割れば合成した平均が得られる。

数式にすると、こんな感じ。

$\bar{x_C} = \frac { \bar{x_A} \times n_A + \bar{x_B} + n_B }{ n_A + n_B }$

ここで $\bar{x_A}$ と $\bar{x_B}$ は各グループの平均値を表している。そして $n_A$ と $n_B$ は各グループの人数を表している。

実際に Python を使って計算してみよう。

>>> C_average = (A.mean() * 10 + B.mean() * 15) / (10 + 15)
>>> C_average
69.840000000000003

上記の値は、実際にグループをまとめて計算した平均値と一致している。

>>> C.mean()
69.840000000000003

分散を合成する

分散の合成については、もうちょっと厄介な計算が必要になる。また、分散の計算方法についても詳しく知っておかないといけない。

分散について

まず、最も基本的な分散の計算式は次のようになる。分散というのは、各要素から平均値を引いたもの (偏差) を二乗して足し合わせて、それをデータ点数で割ることで得られる。

$s^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}$

上記で $\bar{x}$ は $x$ の平均値となっているので、計算式で書くとこう。

$\bar{x} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i$

また、分散の公式は次のように式変形できる。この式は、ようするに各要素を二乗して足し合わせてデータ点数で割ったものから、平均値の二乗を引いて計算することができる、ということ。

$s^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i^{2} - \bar{x}^{2}$

また、上記の「二乗したものを足し合わせてデータ点数で割る」という処理に着目してみよう。これは言い換えると二乗したものの平均を計算していることになる。バー記号は平均を表現する。

$\frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i^{2} = \bar{x^{2}}$

つまり、分散は次のような式でも表現できることが分かった。ようするに、分散は各要素を二乗した平均値 ( $\bar{x^{2}}$ )から各要素の平均値の二乗 ( $\bar{x}^2$ ) を引くことで計算できる。

$s^{2} = \bar{x^{2}} - \bar{x}^2$

分散を合成する手順

さて、ここまでの説明で分散の計算方法は分かった。改めて、求めなければいけないものの式を書いてみよう。

$s_C^{2} = \bar{x_C^{2}} - \bar{x_C}^2$

上記の中で $\bar{x_C}$ については、既に平均の合成をしたところで計算が済んでいる。問題は $\bar{x_C^{2}}$ の方で、これはまだ未知の値なので計算しなきゃいけない。

しかし $\bar{x_C^{2}}$ は、どうすれば求まるだろうか？実は、これも先ほどの分散の公式から導くことができる。

$s^{2} = \bar{x^{2}} - \bar{x}^2$

上記の中で、各グループ毎の $s^{2}$ と $\bar{x}$ は既に分かっている。つまり、式を移行して $\bar{x^{2}}$ を求めるようにできるということ。

$\bar{x^{2}} = s^{2} + \bar{x}^2$

次のように、各グループごとに二乗した平均値を求めることができる。

$\bar{x_A^{2}} = s_A^{2} + \bar{x_A}^2$

$\bar{x_B^{2}} = s_B^{2} + \bar{x_B}^2$

各グループごとの二乗した平均値を求めることができれば、あとは簡単で最初にやった平均の合成をすれば良い。

$\bar{x_C^{2}} = \frac{ \bar{x_A^{2}} \times n_A + \bar{x_B^{2}} \times n_B }{ n_A + n_B }$

上記が計算できれば合成した分散を計算するのに必要な要素が揃うので、あとは定義通りにやるだけ。

$s_C^{2} = \bar{x_C^{2}} - \bar{x_C}^2$

確かめてみる

それでは、上記の式を Python で実行して確かめてみよう。

まずはグループ A の二乗した平均値から求める。

>>> A_square_average = A.var() + A.mean() ** 2
>>> A_square_average
4518.9000000000005

上記は、次の式に対応している。

$\bar{x_A^{2}} = s_A^{2} + \bar{x_A}^2$

そして、同じようにグループ B の二乗した平均値を求める。

>>> B_square_average = B.var() + B.mean() ** 2
>>> B_square_average
5312.333333333333

対応する式は次の通り。

$\bar{x_B^{2}} = s_B^{2} + \bar{x_B}^2$

そして、グループをまとめたときの二乗した平均値を求める。

>>> C_square_average = (A_square_average * 10 + B_square_average * 15) / (10 + 15)
>>> C_square_average
4994.96

対応する式はこちら。

$\bar{x_C^{2}} = \frac{ \bar{x_A^{2}} \times n_A + \bar{x_B^{2}} \times n_B }{ n_A + n_B }$

上記で必要な材料は揃ったので、あとは定義通りに計算するだけ。

>>> C_var = C_square_average - C_average ** 2
>>> C_var
117.33439999999973

実際に計算した値と比較してみよう。

>>> C.var()
117.3344
>>> C_var
117.33439999999973

微妙に誤差が出てるけど、バッチリ合ってるね！

まとめ

二つのグループの平均と分散とデータ点数だけが分かっている状況であってもグループ全体の平均と分散は計算できる

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

作者: もみじあめ
発売日: 2016/03/12
メディア: Kindle版
この商品を含むブログを見る

2016-12-17

Python: 多変数の関数から勾配法で最小値を探す

Mac OS X NumPy Python matplotlib 機械学習

以前、このブログで一変数の関数から勾配法で最小値を探す記事を書いた。

blog.amedama.jp

このときは題材として $f(x) = x^{2}$ という一変数の関数を扱った。今回は、これを多変数の関数に拡張してみることにする。ちなみに、この多変数というのは機械学習における多次元と同じ意味を持っている。つまり、これができると多次元のデータセットで損失関数の出力から最小値を探すことができることを意味している。

使った環境は次の通り。

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

下準備

まずは、前回と同じように可視化に使う matplotlib と、数値計算用のライブラリである numpy をインストールしておこう。

$ pip install matplotlib numpy

多変数の関数を解析的に偏微分する

今回題材として扱う関数は $f(x_0, x_1) = x_0^{2} + x_1^{2}$ にする。

勾配法では最初に関数を微分する必要があった。これは、関数が多変数になったときも変わらない。ただし、多変数になると名前だけは偏微分と呼ばれるようになる。また、変数が複数あるので、それぞれの関数について微分することになる。

試しに、題材とする関数を $x_0$ について解析的に偏微分してみよう。このとき、微分する対象ではない変数 $x_1$ は固定値と捉える。そのため、結果は $2x_0$ となる。

$\frac{ \partial f }{ \partial x_0 } = 2x_0$

次のサンプルコードでは $x_0 = 3.0, x_1 = 4.0$ のときの $\frac{ \partial f }{ \partial x_0 }$ を計算している。とはいえ、上記の内容から結果は $2 \times 3 = 6$ とすぐに分かる。

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


def f(x0, x1):
    """偏微分したい関数: f(x) = x_0^2 + x_1^2"""
    return x0 ** 2 + x1 ** 2


def df_x0(x0, x1):
    """関数 f の解析的な x0 についての偏微分: f'(x) = 2x_0"""
    return 2 * x0  # x_1 は固定値と捉えて微分することで消滅しており登場しない


def main():
    grad = df_x0(3.0, 4.0)  # ぶっちゃけ x_1 は影響を与えないので何を使っても良い
    print('解析的な x0 についての偏微分: {0}'.format(grad))


if __name__ == '__main__':
    main()

実行結果は次の通り。まあ、当たり前だね。

解析的な x0 についての偏微分: 6.0

多変数の関数から数値解法で偏微分の近似値を得る

先ほどの例では人が数式を見て解析的に偏微分したけど、これだと機械学習アルゴリズムには適用が難しい。近似値で良いから、先ほどの偏微分をコンピュータにやらせたい。そこで、今度は前回のブログでも扱った中央差分を元にした数値微分を適用してみることにしよう。

次のサンプルコードを見てほしい。ここでは、先ほどと同じように $x_0 = 3.0, x_1 = 4.0$ のときの $\frac{ \partial f }{ \partial x_0 }$ を計算している。ポイントとしては、偏微分というのは微分に使わない変数を固定値として扱うところだ。なので $x_1$ を 4.0 に固定した関数 f_x0 を定義している。そして、その関数を前回のブログ記事にも登場した数値微分の関数 numerical_diff で微分している。

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


def f(x0, x1):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x0 ** 2 + x1 ** 2


def f_x0(x0):
    """関数 f を x0 について偏微分するために x1 を固定値にした関数"""
    return f(x0, 4.0)


def numerical_diff(f, x):
    """中央差分を元に数値微分する関数"""
    h = 1e-4
    return (f(x + h) - f(x - h)) / (2 * h)


def main():
    grad = numerical_diff(f_x0, 3.0)
    print('数値微分を使った x0 についての偏微分: {0}'.format(grad))


if __name__ == '__main__':
    main()

上記を実行してみよう。

数値微分を使った x0 についての偏微分: 6.00000000000378

おお！解析的な偏微分で得られた値とほとんど同じ値 (近似値) が得られている！

NumPy を使って汎用的にしてみる

先ほどのサンプルコードでは偏微分をするために、微分に関係しない変数を固定値にした関数を新たに用意していた。しかし、これを毎回やるのは大変なので、今度は数値微分の関数 numerical_diffを偏微分用に拡張してみよう。

次のサンプルコードでは、偏微分する関数 f の引数を変更している。具体的には、これまで別々の引数としていた $x_0$ と $x_1$ を NumPy の配列で渡すことにしている。また、偏微分に拡張した関数 numerical_diff には引数 i を追加している。これは、引数 x のうち何番目の要素で偏微分するかを表すインデックスになっている。例えば $i = 0$ なら $x_0$ について偏微分することを意味している。

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

import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def main():
    args = [3.0, 4.0]  # 偏微分に使う引数
    grad = numerical_diff(f, args, 0)  # 0 番目の要素 (x0) について偏微分する
    print('数値微分を使った x0 についての偏微分: {0}'.format(grad))


if __name__ == '__main__':
    main()

ポイントとしては偏微分する要素にだけ丸め誤差で無視されない程度に小さな値を加えて数値微分している点だ。

上記を実行してみよう。

数値微分を使った x0 についての偏微分: 6.00000000000378

先ほどと同じように上手く偏微分できていることが分かる。

偏微分を元に勾配を計算する

多変数の関数において、各変数について偏微分して得られた値をベクトルにまとめたものを勾配と呼ぶ。次は、この勾配を計算してみよう。とはいえ、ある変数について偏微分するやり方は先ほどの例で分かったので、あとはそれを繰り返すだけ。

次のサンプルコードでは numerical_gradient という関数で購買を計算している。ポイントは特に無くて、各変数について淡々と偏微分するだけ。

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

import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def numerical_gradient(f, x):
    """勾配を計算する関数

    勾配というのは、全ての変数について偏微分した結果をベクトルとしてまとめたものを言う。
    """
    # 勾配を入れるベクトルをゼロで初期化する
    grad = np.zeros_like(x)

    for i, _ in enumerate(x):
        # i 番目の変数で偏微分する
        grad[i] = numerical_diff(f, x, i)

    # 計算した勾配を返す
    return grad


def main():
    grad = numerical_gradient(f, np.array([3.0, 4.0]))
    print('勾配: {0}'.format(grad))


if __name__ == '__main__':
    main()

実行すると $x_0 = 3.0, x_1 = 4.0$ のときの勾配が得られる。

勾配: [ 6.  8.]

題材となる関数を可視化してみる

さて、先ほどのサンプルコードで勾配を計算することができた。とはいえ、何のために勾配を計算したのかよく分からないかもしれないので、ここで補足しておく。先ほど計算した勾配には、それぞれの変数について偏微分したときの傾きが入っている。これはつまり、各次元ごとのどちらに進めば値が小さくなるかという情報だ。

上記の説明だけでは分かりづらいかもしれないので、より原点に立ち返った説明もしてみよう。例えば $x_0, x_1, y$ という三つの軸を持った三次元の空間をイメージしてほしい。通常であれば変数の名前は $x, y, z$ になるけど、ここではあえて今回使った変数の名前にしている。まず $x_0$ と $x_1$ が関数の入力で $y$ が出力だとしよう。縦と横が $x_0$ と $x_1$ に対応していて高さが $y$ ということになる。今回の問題は高さ $y$ が最も低くなる位置を $x_0$ と $x_1$ の平面の中から探したいということだ。

試しに、今回題材とする関数 $f(x) = x_0^{2} + x_1^{2}$ を上記のように三次元で捉えて可視化してみることにしよう。

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

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


def main():
    x0 = np.arange(-2, 2, 0.1)
    x1 = np.arange(-2, 2, 0.1)
    X0, X1 = np.meshgrid(x0, x1)
    Y = X0 ** 2 + X1 ** 2

    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_wireframe(X0, X1, Y)

    plt.show()


if __name__ == '__main__':
    main()

上記を実行すると、次のようなグラフが得られる。 f:id:momijiame:20161217194456p:plain 上記より、最小値は $x0 = 0, x1 = 0$ ということが分かる。

偏微分を元にした勾配の計算というのは、上記のグラフでいえば各点における傾きを $x_0, x_1$ という、それぞれの軸ごとに計算するという意味になる。

勾配法への組み込み

ここで、前回も扱った勾配法の数式を多変数に拡張して示す。

$x_{ij + 1} = x_{ij} - \eta f_{x_i}'(x)$

上記で現在の位置を $x_{ij}$ としたとき、より最小値に近づいたものが $x_{ij + 1}$ になる。 $\eta$ は学習率と呼ばれるもので、どれくらいの勢いで最小値に近づいていくかを表している。そして $f_{x_i}'(x)$ は $f$ を $x_i$ について偏微分した関数になる。

上記は、扱う変数が増えて多変数 (多次元) になったとしてもアルゴリズムの基本は全く変わらない。単に、各次元ごとにそれぞれで勾配法を実行するだけだ。次のサンプルコードでは初期位置を $x_0 = 1, x_0 = 2$ として、学習率 0.01 ステップ数 50 回で勾配法を実行している。

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

from matplotlib import pyplot as plt
import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def numerical_gradient(f, x):
    """勾配を計算する関数

    勾配というのは、全ての変数について偏微分した結果をベクトルとしてまとめたものを言う。
    """
    # 勾配を入れるベクトルをゼロで初期化する
    grad = np.zeros_like(x)

    for i, _ in enumerate(x):
        # i 番目の変数で偏微分する
        grad[i] = numerical_diff(f, x, i)

    # 計算した勾配を返す
    return grad


def gradient_descent(f, initial_position, learning_rate=0.1, steps=50):
    """勾配法で最小値を求める関数

    :param function f: 最小値を見つけたい関数
    :param numpy.ndarray initial_position: 関数の初期位置
    :param float learning_rate: 学習率
    :param int steps: 学習回数
    """
    # 現在地を示すカーソル
    x = initial_position

    # 学習を繰り返す
    for _ in range(steps):
        # 現在地の勾配 (どちらにどれだけ進むべきか) を得る
        grad = numerical_gradient(f, x)
        # 勾配を元にして現在地を移動する
        x -= learning_rate * grad

    # 最終的な位置を返す
    return x


def main():
    # 勾配法を使って関数 f() の最小値を探す (初期位置は 1, 2)
    min_pos = gradient_descent(f, [1, 2])
    print('勾配法が見つけた最小値: {0}'.format(min_pos))


if __name__ == '__main__':
    main()

上記を実行した結果は次の通り。

勾配法が見つけた最小値: [  1.78405962e-05   3.56811923e-05]

最小値である $x_0 = 0, x_1 = 0$ に、どちらの次元についても近づけていることが分かる。

最小値に収束していく様子を可視化してみる

先ほどの結果から、最小値に収束させることができることは分かった。次は、どのように収束していっているのかをグラフで可視化してみよう。

次のサンプルコードでは各ステップごとの $x_0$ と $x_1$ の位置を折れ線グラフで示している。

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

from matplotlib import pyplot as plt
import numpy as np


def f(x):
    """偏微分したい関数: f(x) = x0^2 + x1^2"""
    return x[0] ** 2 + x[1] ** 2


def numerical_diff(f, x, i):
    """中央差分を元に数値微分する関数 (偏微分)

    :param function f: 偏微分する関数
    :param numpy.ndarray x: 偏微分する引数
    :param int i: 偏微分する変数のインデックス
    """
    # 丸め誤差で無視されない程度に小さな値を用意する
    h = 1e-4
    # 偏微分する変数のインデックスにだけ上記の値を入れる
    h_vec = np.zeros_like(x)
    h_vec[i] = h
    # 数値微分を使って偏微分する
    return (f(x + h_vec) - f(x - h_vec)) / (2 * h)


def numerical_gradient(f, x):
    """勾配を計算する関数

    勾配というのは、全ての変数について偏微分した結果をベクトルとしてまとめたものを言う。
    """
    # 勾配を入れるベクトルをゼロで初期化する
    grad = np.zeros_like(x)

    for i, _ in enumerate(x):
        # i 番目の変数で偏微分する
        grad[i] = numerical_diff(f, x, i)

    # 計算した勾配を返す
    return grad


def gradient_descent(f, initial_position, learning_rate=0.1):
    """勾配法で最小値を求める関数

    :param function f: 最小値を見つけたい関数
    :param numpy.ndarray initial_position: 関数の初期位置
    :param float learning_rate: 学習率
    """
    # 現在地を示すカーソル
    x = initial_position

    # 学習を繰り返す
    while True:
        # 現在地の勾配 (どちらにどれだけ進むべきか) を得る
        grad = numerical_gradient(f, x)
        # 勾配を元にして現在地を移動する
        x -= learning_rate * grad
        # 現在の位置を返す
        yield np.copy(x)


def main():
    g = gradient_descent(f, [1, 2])

    x = range(50)
    y = np.array([next(g) for _ in x])

    plt.plot(x, y[:, 0], label='f_x0(x)')
    plt.plot(x, y[:, 1], label='f_x1(x)')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 f:id:momijiame:20161217195558p:plain それぞれの軸で最小値に収束していく様子が見て取れる。

まとめ

勾配法は多変数の関数でも使える
ただし傾きを得るときの計算が偏微分になる
偏微分は解析的に計算する必要はなく数値解法での近似値が使える
各変数について偏微分した結果のベクトルを勾配という
勾配を元に各次元ごとに最小値を探すことができる

参考文献

ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装

作者: 斎藤康毅
出版社/メーカー: オライリージャパン
発売日: 2016/09/24
メディア: 単行本（ソフトカバー）
この商品を含むブログ (6件) を見る

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

作者: もみじあめ
発売日: 2016/03/12
メディア: Kindle版
この商品を含むブログを見る