読者です 読者をやめる 読者になる 読者になる

CUBE SUGAR CONTAINER

技術系のこと書きます。

統計: 統計検定2級に合格した

統計

先日、当面の目標にしていた統計検定2級に合格することができた。 今回は、受験に関する諸々について書いてみることにする。

受験のきっかけ

以前から、データ分析や機械学習に興味があった。 そして、それらの書籍を読んだり手法を調べていくうちに、だんだんと統計学に対する興味が湧いてきた。 統計学は、データ分析や機械学習に深い関わりがある。

その後は、初心者向けの統計に関する書籍などを読んで学び始めた。 とはいえ、それだけでは統計が身についているのかがよく分からない状況に陥ってしまう。 そんなとき、統計検定の存在を知った。

ウェブで下調べしたところ、概ね統計検定の 2 級に合格すれば統計の基礎は分かっていると胸を張れるらしい。 現在、統計検定は 1 級 (数理・応用)、準 1 級、2 級、3 級、4 級と五つのレベルに分かれている。 公式では 2 級の試験内容を「大学基礎課程で習得すべきこと」と位置づけている。

www.toukei-kentei.jp

その上で出題範囲を見ると、なかなか広いので体系的な知識が得られそうな感じ。 勉強する過程で例題をたくさん解くことになるので身にもつきやすそうだ。 そこで、これをひとまず当面の目標に勉強を進めてみることにした。

以上が、統計検定2級を受験したきっかけ。 ちなみにスタート地点としては「標準偏差ってどう計算するんだっけ?」というレベルだった。 大学は情報系の学科だったものの、統計については全く触れていない。 ちなみに、今は高校の数学でも統計が範囲に入っているらしい。 統計検定でいうと 3 級に対応するかな。

やってきた勉強

次に、合格までにどんな勉強をしたかについて。 念のため、統計検定を意識する以前の内容についても書いておく。

完全独習 統計学入門

まず、統計学を学び始めるとっかかりとして、次の本を読み始めた。

完全独習 統計学入門

完全独習 統計学入門

上記の本は、統計に関する本を色々と読んだ中でも抜群に分かりやすかった。 本の中に、自身の説明として「これ以上内容を削ると統計ではなくなる」とあるけど、正にという感じ。 最初の一冊としては、ぜひおすすめしたいものになっている。

この本は二部構成になっていて、前半では記述統計を、後半では推測統計を扱っている。 前半の記述統計というのは、データの要約や可視化などを扱うもの。 この本では、基本的な統計量 (平均・分散・標準偏差など) の知識と正規分布について学ぶことができる。 後半の推測統計というのは、標本 (サンプル) を元に母集団を推測するもの。 この本では、ステップバイステップで正規分布・カイ二乗分布・t分布を用いた推定と検定について学ぶことができる。

ちなみに、推定と検定は統計検定2級でも必ず出題される。 その割に分かりやすく解説しているものは少ないので、この本はその意味でも貴重な存在といえる。 基礎的な考え方さえ身につけば、あとはその応用なので別の本も読みやすくなる。

完全独習 ベイズ統計学入門

読んだのはだいぶ後だけど、同じ作者でこの本もおすすめできる。 先ほどの本と何が違うかというと、この本はベイズ推定に焦点を当てた説明になっている。

完全独習 ベイズ統計学入門

完全独習 ベイズ統計学入門

実は統計学はスタンダードな統計 (ネイマン・ピアソン統計とも呼ばれる) とベイズ統計に分かれている。 先ほどの本はスタンダードな統計を扱っていて、こちらはベイズ統計を扱っているという違い。

統計検定2級では確率の問題でベイズ推定も頻繁に出題されるので、この本は強い味方になる。 確率を面積で計算するテクニックを身につけると問題を解くのに大いに役立つだろう。

統計検定 2級 公式問題集

次は、みんな大好き過去問題集。 この問題を解きながら、分からないところを解説を読みながら調べていくのが勉強の基本だった。

日本統計学会公式認定 統計検定 2級 公式問題集[2013〜2015年]

日本統計学会公式認定 統計検定 2級 公式問題集[2013〜2015年]

最初に取り組んだときはあまりに解けなくて絶望したけど、苦手な分野を早めに把握するのが重要。 分からなかったところは、付属している解説を読んだり後述する公式の教科書を参照しながら進めよう。

改訂版 統計検定2級対応 統計学基礎

統計検定2級の教科書的な位置づけの本がこちら。

改訂版 日本統計学会公式認定 統計検定2級対応 統計学基礎

改訂版 日本統計学会公式認定 統計検定2級対応 統計学基礎

  • 作者: 田中豊,中西寛子,姫野哲人,酒折文武,山本義郎,日本統計学会
  • 出版社/メーカー: 東京図書
  • 発売日: 2015/12/10
  • メディア: 単行本
  • この商品を含むブログを見る

ただ、これを頭から読んで勉強していくのはおすすめしない。 なぜなら、ストーリー性が皆無なので順番には読みづらいから。 個人的には辞書として使うのがおすすめ。 分からなかったところをピンポイントに読んで理解の助けにするのが良い。

統計検定3級対応 データの分析

最初に紹介した「完全独習 統計学入門」は素晴らしい本なんだけど、ちょっと残念なところもある。 それは、本のコンセプトからして仕方ないんだけど、推定や検定を説明するまでが最短距離すぎて、色々と省いているところ。 例えば色々な種類のグラフとか、箱ひげ図と五数要約みたいな一般的な話が出てこない。 そういった点では、この本で前提知識をまとめてインプットしておくと良いと思う。

日本統計学会公式認定 統計検定3級対応  データの分析

日本統計学会公式認定 統計検定3級対応 データの分析

  • 作者: 藤井良宜,竹内光悦,後藤智弘,日本統計学会,竹村彰通,岩崎学,美添泰人
  • 出版社/メーカー: 東京図書
  • 発売日: 2012/07/07
  • メディア: 単行本(ソフトカバー)
  • 購入: 1人 クリック: 65回
  • この商品を含むブログ (6件) を見る

この本は割りとサクサク読めるので頭から教科書的に使っても大丈夫そうだった。

高校数学の美しい物語

こちらは書籍ではなくウェブサイトだけど、統計に関する記事が豊富にある。 特定のトピックについてよく分からないなーと思ってぐぐると、こちらがヒットする確率が高い。

mathtrain.jp

役に立つ薬の情報~専門薬学 (統計)

こちらも書籍ではなくウェブサイトだけど、ぐぐるとヒットする確率が高い。 特にカイ二乗検定やt検定の説明の分かりやすさはピカイチだと思う。

統計学

ブログ執筆

あとは、理解があやふやだったり重要だと思ったところはこのブログで積極的に記事にした。 他の人に説明することは自分が理解していないとできないことなので、解説することは理解につながる。

受験について

統計検定は年に二回、夏と冬に定期的に開催されている。 ただし、2級と3級に関しては CBT 方式の試験も随時受けることができる。

www.toukei-kentei.jp

今回合格したのは、上記の CBT 方式の試験だった。 公式サイトにも記述があるけど、受けた後の感想としては問題の難易度は定期試験のそれと変わらない。 だいたい 35 問くらいが出題されて、そのうち 60% (21 問) を正解すれば合格となる。

CBT 方式の試験では、あらかじめテストセンターに予約を入れて、その時間に現地に赴いて受験する。 今回は有楽町にあるオデッセイ公式のテストセンターで受けた。

otc.odyssey-com.co.jp

持参するものは受験票と電卓と身分証だけで、紙とペンは会場で用意してもらえる。 あとはディスプレイに表示される問題を解いて、ひたすら答えをマウスとキーボードで入力していく。

実際に受けた感想としては (他の人はどうか分からないけど) 個人的には試験時間が全然足りなかった。 選択肢から、あらかじめありえないものは枝刈りするなどして時間を節約していくのが重要そう。

受験を終えて

勉強を通して、統計について最低限知っておくべき知識は身につけることができたと思う。 現実世界でも使われている統計の色々な手法について知ることができたのは面白かった。 例えば、選挙の当選確実をどうやって出しているか (母比率の推定) とか、治験の効果の調べ方 (一標本t検定) とか。

とはいえ、あくまで最低限なのは確かで、より専門的な知識をこれから身につけていかないとなという感じ。 これからは、より実務的に使われているデータ分析や機械学習の手法についての勉強にシフトしていこうと思う。

いじょう。

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

めでたしめでたし。

統計: 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.050.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 検定が使える。

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

統計

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

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 分布を使って推定する

参考文献

完全独習 統計学入門

完全独習 統計学入門

統計: はじめての推定

統計

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

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

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

記述統計と推測統計

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

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

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

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

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

どうやって推定するか?

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

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

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

最も単純な推定

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

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

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

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

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

平均について

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

算術平均の定義は次の通り。 ここで  \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}} に縮む

参考文献

完全独習 統計学入門

完全独習 統計学入門

統計: 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) で線形モデルを使った重回帰分析をしてみた。

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

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.9176.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

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

まとめ

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