タグ「統計学」が付けられているもの

偏相関係数の式の導出をやってみた

偏相関係数とは、XとZ、YとZに相関がある時の、Zの影響を抜いたXとYの相関係数である。XとYの相関係数の絶対値が大きくても、それぞれが間接的な要因Zと因果関係がある場合、それによって相関が強くなっている可能性がある。

そのような場合に、Zの影響を除いたXとYの相関係数を知りたくなる。
相関係数が高いが実はXとYに因果関係が無く、それぞれ間接的な要因Zと因果関係がある疑似相関（見せかけの相関関係）の例としては、アイスクリームの売上と水難事故の発生数の相関（気温が間接要因）や、警察官の人数と犯罪の発生数の相関（人口が間接要因）などが知られているようで、筆者も何度も聞いた覚え、見た覚えがある。

XとY、XとZ、YとZの相関係数をそれぞれrxy, rxz, ryzと書くと、Zの影響を除いたXとYの偏相関係数rxy⋅zは次の式で表される、と色々な所に書いてある。
r_{xy\cdot z}=\frac{r_{xy}-r_{xz}r_{yz}}{\sqrt{\vphantom{r_{yz}^2}(1-r_{xz}^2)}\sqrt{(1-r_{yz}^2)}}

何故こうなるのだろうと思って、筆者が30年前に通っていた時の大学の教科書を読み返したら、「簡単な計算により（上の式）を得る。」とだけ書かれていて、導出過程が載っていなかった。
しかし、その上に考え方は書いてあったので、そこから計算してみた。
いくつか前提知識の復習が必要で、やり始めると面倒くさそうで気が削がれるが、一踏ん張りすれば、確かに簡単だった。

考え方としては、X,YをそれぞれZの一次関数として回帰すると、それらのZの項がZの影響の部分なので、それぞれのZの項の部分を抜いたもの同士の相関係数を偏相関係数とするということらしい。定数項を抜いても相関係数は変わらないので、X,Yそれぞれの残差同士の相関係数ということになる。

以下、Xの標本をxi 、Xの標本標準偏差をsx 、XとYの標本共分散をsxy 、XとYの標本相関係数をrxyのように表記する。
XをZで一次回帰すると
x=\bar{x}+\frac{s_{xz}}{s_z^2}(z-\bar{z})
となることから、X,YをZで回帰した時の残差X',Y'は次のようになる。
\matrix{ x'_i=x_i-\bar{x}-\frac{s_{xz}}{s_z^2}(z_i-\bar{z}) \cr y'_i=y_i-\bar{y}-\frac{s_{yz}}{s_z^2}(z_i-\bar{z}) }
求めるものは、相関係数の定義より
r_{x'y'}=\frac{s_{x'y'}}{s_{x'}s_{y'}}
である。このsx'y', sx', sy'をそれぞれ計算していく。
共分散の定義より、
s_{x'y'}=\frac{1}{n}\sum_i(x'_i-\bar{x'})(y'_i-\bar{y'})
であり、一次回帰式の誤差はN(0,σ2 )に従う前提（最小二乗推定量の最小分散性は正規性を仮定している）であることから、X'とY'の平均は0なので、
= \frac{1}{n}\sum_i x'_iy'_i = \frac{1}{n}\sum_i\left\{(x_i-\bar{x}-\frac{s_{xz}}{s_z^2}(z_i-\bar{z}))(y_i-\bar{y}-\frac{s_{yz}}{s_z^2}(z_i-\bar{z}))\right\}
= \frac{1}{n}\sum_i\left\{ (x_i-\bar{x})(y_i-\bar{y})-\frac{s_{yz}}{s_z^2}(x_i-\bar{x})(z_i-\bar{z})-\frac{s_{xz}}{s_z^2}(y_i-\bar{y})(z_i-\bar{z})+\frac{s_{xz}s_{yz}}{s_z^4}(z_i-\bar{z})^2 \right\}
=s_{xy}-\frac{s_{xz}s_{yz}}{s_z^2}-\frac{s_{xz}s_{yz}}{s_z^2}+\frac{s_{xz}s_{yz}}{s_z^2}=s_{xy}-\frac{s_{xz}s_{yz}}{s_z^2}
sxy=sxsyrxyなので
=s_xs_yr_{xy}-s_xs_yr_{xz}r_{yz}=s_xs_y(r_{xy}-r_{xz}r_{yz})
となる。同じ要領で、
s_{x'}^2=\frac{1}{n}\sum_i x'^2_i=\frac{1}{n}\sum_i \left\{(x_i-\bar{x}-\frac{s_{xz}}{s_z^2}(z_i-\bar{z}))^2\right\}
=\frac{1}{n}\sum_i \left\{(x_i-\bar{x})^2-2\frac{s_{xz}}{s_z^2}(x_i-\bar{x})(z_i-\bar{z})+ \frac{s_{xz}^2}{s_z^4}(z_i-\bar{z})^2\right\}
=s_{x}^2-2\frac{s_{xz}^2}{s_z^2}+\frac{s_{xz}^2}{s_z^2}=s_{x}^2-\frac{s_{xz}^2}{s_z^2}=s_{x}^2-s_{x}^2r_{xz}^2=s_{x}^2(1-r_{xz}^2)
となる。sy'についてはxとyを入れ替えただけなので省略する。従って、
r_{x'y'}=\frac{s_xs_y(r_{xy}-r_{xz}r_{yz})}{s_x\sqrt{\vphantom{r_{yz}^2}(1-r_{xz}^2)} s_{y}\sqrt{(1-r_{yz}^2)}}=\frac{r_{xy}-r_{xz}r_{yz}}{\sqrt{\vphantom{r_{yz}^2}(1-r_{xz}^2)} \sqrt{(1-r_{yz}^2)}}
である。

F分布は分母の自由度が大きいとχ2分布に近似できる？

「統計検定 ２級 公式問題集 2017～2019年」の解説部分のコメントの中に、

自由度(5,77)のF分布の上側1%点は付表には与えられていないため、近似的に考える。すなわち、残差の自由度77が十分大きいので、F統計量に5を掛けた統計量の分布が自由度5のカイ二乗分布で近似できると考えてよい。
というのがある（残差というのは分散分析の文脈のもので、分母のこと）。最初にこれを読んだ時は理解できた気がして、試験前の復習の為にメモっておいたのだが、試験後にメモ書きを捨てる前にメモを1つ1つ読んでると、このことだけ理解できなかった。

確かに、同書のα=0.05のF分布表の一番下の行（分母の自由度が一番大きい所）とカイ二乗分布表のα=0.05の列を比較すると、F分布表の数値に分子の自由度を掛けた値がカイ二乗分布表の値と大体一致するし、別の本の自由度が∞の行があるF分布表の∞の行と同様に比較すると完全に一致している。

自由度(5,100)のF分布のFを5倍した確率密度と自由度5のχ2分布の確率密度のグラフを比較しても、確かにほぼ一致している。自由度(5,10)のものなどと比較すると、分母の自由度が大きいほどχ2分布に近づくことが見て取れる。

しかし、自由度(m,n)のF分布に従うFは、X,Yがそれぞれ独立で自由度m,nのχ2分布に従うとすると
F = \frac{X / m}{Y / n}
であり、例えばn=100のχ2分布の確率密度は

なので、Yの期待値は1であるがY=1ではないと思った。

自由度mのχ2分布の確率密度関数は
f(x) = \frac{1}{\Gamma(\frac{m}{2})2^{\frac{m}{2}}} x^{\frac{m}{2}-1} e^{-\frac{x}{2}}
で、自由度(m,n)のF分布の確率密度関数は
g(x) = \frac{\Gamma(\frac{m+n}{2})}{\Gamma(\frac{m}{2})\Gamma(\frac{m}{2})} \left(\frac{m}{n}\right)^\frac{m}{2} x^{\frac{m}{2}-1} \left(\frac{m}{n}x+1\right)^{-\frac{m+n}{2}}
だそうなので、n→∞の時にこれらが一致することを確認するしかないか...と思ったが、ネットで調べてみると、大体
\lim_{n\to\infty} Y/n = 1
で片付けられていた。

もう少し丁寧な説明だと、それぞれ独立の、自由度1のカイ二乗分布に従うYiを使って
\frac{Y}{n} = \frac{Y_1 + Y_2 + \cdots + Y_n}{n}
とすると、大数の法則（大数の強法則）より
\lim_{n\to\infty} \frac{Y}{n} = E[Y_1] = 1
なので
\lim_{n\to\infty} mF = \lim_{n\to\infty} \frac{X}{Y/n} = X
とのことである。

統計検定2級 CBT方式試験受験

8/14（土）に統計検定2級のCBT(Computer Based Test)方式の試験を受けた。 その経緯や感想を書き残す。

結果は76点で合格だったが、筆者としては不本意な結果だった。

■筆者のバックグラウンド
  • 47歳
  • 子供の頃から数学好き
  • 大学時代は理系の学部
    統計学に興味があり、教養課程と大学院の修士課程で統計学を計3回履修したが、推定・検定がさっぱり理解できず
  • 大学卒業後はソフトウェア開発の仕事で全く統計学に触れず
  • 35歳の頃にふと思い立って入門書で統計学を復習
    推定・検定をやっと理解して本weblogに色々書いたが、結局初学者に少し毛が生えた程度
  • 46歳の頃に仕事がデータ分析に変わる
    ある時、会議で推定か検定を使った資料が出てきて、質問しようとしたが、すっかり忘れており、「信頼区間」「有意水準」という言葉すら出て来なかった
■受験のきっかけ

今年の4月に会社から、業務分野に関係あるかも知れない試験として

  • Python 3 エンジニア認定データ分析試験
  • AWS認定試験のどれか
  • 情報処理技術者試験のどれか
  • 統計検定2級
が紹介されたこと。実務的には上の2つが合ってそうだったが、統計学に興味があった私の目は「統計検定2級」しか認識しなかった。

通常のPBT(Paper Based Test)試験は今年は大阪会場が無く、来年以降は1級以外PBT方式試験をしないとのことなので、諦めてCBT試験を受けることにした。

■やったこと

4月にたまたま図書館で「統計学の図鑑」（技術評論社）という本を見つけ、借りて復習のつもりでのべ10時間くらいで一通り読んだ。
GWに「統計検定 ２級 公式問題集 2017～2019年」の過去問にトライしたら、時間無制限で少しずつのトライであったが、80%正解したので、目標を優秀「S」判定の基準と言われる90%に置いた。解説を読むと、知らないことが多かったので、「統計学の時間 | 統計WEB」で一通り勉強した。このサイトは大変わかりやすく、ただ感謝である。
それから過去問を2つ解いたが、目標の90%に届きそうに無かったので、もう一度上記サイトを一通り読んだ。
後は過去問をやり続けた。

過去問（1周目）の正答率はこんな感じだった。

過去問やった日正答率所要時間
2019年11月分 5/3 80% 180分
2019年6月分 6/13 86% 150分
2018年11月分 6/19 85% 135分
2018年6月分 7/11 97% 135分
2017年11月分 7/17 85% 100分
2017年6月分 7/22 94% 76分

8月に本番と同じ90分の制限時間を設けて2周目を一通り（6回分）やったが、ある程度問題も答えも覚えてしまっていた為、全て時間内に解け、91〜100%正解した。
最新1回分の試験問題と正解がWebに掲載されているようだが、この投稿を書くまで気付かなかったので、やらなかった。

最初の感想としては、80%正解するのは容易だが90%はかなり難しいと思った。知識不足で何か知らないことが出てきてしまうのは仕方ないとして、知識があって間違えなければ解ける問題も何問か間違えてしまう。問題文の読み損ねや読み間違い（ひっかけ問題にひっかかるのを含む）が特に多いが、電卓操作ミス、マークミスも最後まで出続けた。

勉強時間は、過去問に取り組むのを含めて週末ごとに4〜5時間のペースで、計100時間くらいだったと思う。

過去問に取り組み始めて一番困ったのは、問題文を読むのが辛いことだった。当初は問題文がすっと頭に入って来ず、問題文を見た途端に身体が拒否するようになり、無意識に違うことを考え始めてしまうのが最大の課題だった。筆者は長年、こういう類の文章を素早く次々と読んで理解する必要性が無く、特に数学系の文章問題はほぼ大学以来なので慣れていないのである。
加えて、筆者はこの歳になってから急に身体中にガタが来て、1つの困った症状として、集中できなくなることが増えた。

まだ理解が怪しかったり、記憶が定まらない所があり、もう少し勉強してから受験しようかとも思っていたが、覚えたことを片っ端から忘れ、1ヶ月前に覚えたことは記憶に残ってないような状況が続いたので、肉体的にこれが限界だと思って会社の夏季休暇中に受験することにした。

■本番

最初に全32問と表示されて、34〜35問だと思っていたので、少な目なのだなと油断してしまった。それが原因か、時間管理を誤って大失敗してしまった。

最初の方に確率の問題が多かった。問題のレベルも少し難しいように感じた。確率に関しては2017〜19年の過去問とは全く異なる傾向だったように思う。3つ目の確率の問題で少し手こずって時間を失った後、次の問題がまた確率の問題で、かなり嫌な気分になった。しかも出した答えが選択肢に無く、解き直したりしていると、気付かない内にかなりの時間を使ってしまったようで、その問題は諦めたが最初の7問くらいで30分が過ぎてしまった。
そこからは少し急いだつもりが、12問目を終えた辺りで45分が過ぎてしまい、とんでもないことになったと思った。
そこからは問題文を飛ばし読みで急いで解いたが、26問目を終えて残り12分、そこからも到底1問2分で解けないような問題が続き、最後の大問4つ、選択肢にして6つは問題文をほとんど読まずに回答する羽目になってしまった。

確率に限らず、全体的に2017〜19年の過去問とは出題傾向が違うと感じた。また、推定・検定で時間があっても筆者には解けそうになかった問題も2つくらいあった。

何よりも、問題文の文字数が多かった。CBT試験特有なのだろうか。2017〜19年の過去問では大問が15〜18問で問題数が34〜35なのに対し、大問が32問で問題数が34〜35だったので、前の問題の続きというのがほとんど無く、1問ずつ長文が出てくるのである。
終盤、1問当たり2分くらいしか残っていない状況でも画面いっぱいに文字だらけの問題文が次々と出てきて、読む気が失せてしまった。
最初に出てきた全32問というのが大問のことだと気付かず、問題数が32だと思って油断してしまったが、逆に大問数が倍近いので警戒しないといけなかった。

問題数を勘違いして油断したのと、時間配分の練習が十分でなかったのと、問題文を速く読む練習及び長文を読む集中力が足りなかったのが敗因だった。

■その他、CBT試験について思ったこと

紙の試験と異なり、パソコンの画面を見ながら電卓を操作するのが大変だった。視点を大きく動かすので、先ほど見ていた場所を見失ってしまう。つい今打ち込んだ数字を画面上で探してしまうが、同じ数字が2つあるとわからなくなる。

それから、画面上では小さい文字が判別できない。SxxなのかSxyなのか判別できず、画面を拡大して戻す操作が必要になることがあった。モニターのせいかフォントのせいかわからないが、画面の解像度が低かったようだ。
全体的に、やっぱり紙に印刷された問題より問題文が読みにくいと感じた。画面を凝視すると目が疲れるので、筆者の場合は自然に目を逸らしてしまう。視線の高さが異なるのも一因だろう。

それから、CBT試験は正解が知らされない（PBT試験はWebサイトに正解が発表される）ので、どの問題を間違えたのかがわからないのもデメリットである。

今回、7月下旬にCBT試験の会場を探したら、4月に探した時はあった市内の会場が無くなっていて、さらに新型コロナウィルスの影響か、自宅から行きやすい順に3ヶ所くらいに問い合わせたら、揃って今は休止中と言われた。別のもう1ヶ所は予約不可、満席なら諦めてというあまり力を入れていない感じだったので、やめた。
結局、自宅から行きやすい所を選ぶのを諦め、CBT試験会場として評判が良さそうな難波のパソコン教室を選択した。実際、親切に対応して頂いたし、特に気になることは無く、良い会場だった。
CBT試験は時間帯や会場を選べるとの触れ込みだが、タイミングが悪かったのかも知れないが、会場を選ぶのには結構苦労した。

データサイズが異なる処理時間の差の検定を考える

あるソフトウェアをバージョンアップしたら、処理時間が少し伸びた。
計測結果は、次のような感じだった。

入力データ処理時間(before)処理時間(after)
1300360
2200180
3100120
4500480
5400440
612001300
7250270
8200210
915001800
1010001100

ぱっと見で明らかに悪化しているが、念の為、統計学的に有意な差であることを示しておこうと思って、Excelの分析ツールの「ｔ検定: 一対の標本による平均の検定」（いわゆる対応のあるｔ検定）を使ったら、p値が意外と小さくなかった。例えば上記のデータだと両側検定で0.0706となる。
それでも、有意水準を10%とすると有意差ありと言えないことはないので、不本意ながらその結果を付けて報告した。

そのメールを発信した日の帰り道に、このｔ検定は妥当な計算がなされていないのではないかと疑問に思った。 同じ入力データのセットを使ったのでbeforeとafterで対応があるが、入力データはサイズがまちまちなので、(after-before)の差それぞれは同列に扱えない。(after-before)の単純な差でなく、差分の比率、それぞれ何%増であるかが検定の対象である。 専門的に言うと、ｔ検定は標本の母集団が正規分布に従うのが前提であるが、この(after-before)の差は正規分布に従うはずがないので、ｔ検定を用いるのは不適切だということになるだろう。

もしかして、Excelの分析ツールの「ｔ検定: 一対の標本による平均の検定」はそれも考慮したp値を計算してたりしないだろうか、と期待したが、結果を確かめてみると、一般的な「対応のあるｔ検定」と同様、
t=\frac{\bar{d}}{s/\sqrt{n}}
（dは対応のある数値の差、dの上のバーは平均、sは不偏標準偏差、nは標本データ数）
というｔ値と、それに対応する自由度n-1のｔ分布における外側確率だった。

こういった、差分の比率を検定したいケースはよくあると思うので、何か定番の方法があるのではないかと思って少し探したが、全く見つからなかった。

対応するデータの差分を、正規分布に従うように補正（正規化）してｔ検定すれば良いのだが、よく考えると、一般的には差分をbeforeの値で割ったものが正規分布に従う保証は無いのである。今回のデータが、計算量が入力データのサイズに比例する処理時間であることがわかっているから、差分をbeforeで割ったら大よそ正規分布に従うという仮定を置けるが、もし計算量のオーダーがO(n)でなくO(n2)だったら、差分のオーダーもO(n2)ならbeforeの平方根で割るべきかも知れないし、差分のオーダーが不明なら正規化のしようが無いかも知れない。

今回は結局、対応するデータの差をbeforeとafterの平均で割ることによって正規化し、ｔ検定をやり直した。単純にbeforeの値で割って正規化しても良いだろうが、before→afterとafter→beforeを対称に扱う為、例えば80→100と100→80が相殺されるようにする為にこのようにした。考え過ぎか。
例えば上記のデータだと、このように正規化して計算すると、両側検定のp値が0.0267となる。

統計学復習メモ19: 分散分析の種類

分散分析というと、その名前自体によく「一元配置」「二元配置」とか「対応あり」「対応なし」とか「繰り返しのある」とか「繰り返しのない」とかいう言葉がついて回る。統計学の書籍でも、「一元配置分散分析」と「二元配置分散分析」は項目を分けて説明されることが多い。さらにそれぞれが「対応あり／なし」「繰り返しあり／なし」等で分けられると、目次に「分散分析」がたくさん並ぶことになる。Excelの「分析ツール」のメニューにも分散分析の項目があるが、

  • 分散分析： 一元配置
  • 分散分析： 繰り返しのある二元配置
  • 分散分析： 繰り返しのない二元配置
などと分かれており、試しにやってみようと思っても、それらの違いがわからないとどうしようもない。同じ分散分析でも、それらの選択によっては、結果が的外れな意味を為さないものになる可能性があるのだ。その「分散分析」に付加される言葉の種類の多さに圧倒されて、分散分析は覚えることが多いと思って勉強する気を失ったのは、筆者だけであろうか。

筆者は未だに、どういう時にどの種類の分散分析を使うべきなのかよくわかっていない。 頭の中を整理するために、いつものように体当たり的に、

  • 一元配置（対応なし）
  • 一元配置（対応あり）
  • 二元配置（繰り返しなし）
  • 二元配置（繰り返しあり、対応なし）
  • 二元配置（繰り返しあり、１要因対応あり）
  • 二元配置（繰り返しあり、２要因対応あり）
のそれぞれの適用例を、筆者の職業に関係のあるソフトウェア開発を題材にして、考えてみることにした。

以下のデータは、全て架空のものである。
筆者はExcelを持っていないので、分散分析の計算にはRを使う。Rにも分散分析の関数はoneway.test(), aov(), anova(), lme()など色々あるが、今回は上記の全てをカバーできるaov()を使う。

最もシンプルな、一元配置（対応なし）の例として、次のようなデータを考える。

一元配置（対応なし）の例(1)
要因：モジュールの階層
アプリミドルドライバ
生産性
(steps/人月)		 490
410
590
500
460
690
750
500
770
720
730 		480
650
310
450
530
550
350 		460
370
610
240
500
この生産性はモジュールの階層によって差があると言えるかどうか（差が無いなら滅多にこうはならないこと）を、分散分析で調べる。 
> sample1 <- read.table("anova01.dat", header=TRUE)
> summary(aov(productivity ~ module, data=sample1))
            Df Sum Sq Mean Sq F value  Pr(>F)  
module       2 120939   60470  3.5656 0.04738 *
Residuals   20 339182   16959                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
>
"Df"が自由度、"module"の行の"Sum Sq"が群間平方和（ここではモジュールの階層の違いによる効果）、"Mean Sq"が群間の不偏分散（モジュールの階層の違いによるばらつきの分散値）、Residuals（残差＝群内のばらつき＝全体に共通するはずのばらつき）の行のMean Sqが群内（全群共通）の不偏分散、"F value"が（群間／群内）の分散比、"Pr"（P値、そのF値以上が起こる確率）の右の'*'が有意水準0.05で有意であることを表す。
分散分析の帰無仮説は「全ての標本は同じ分布に従う母集団から得られたもの」→「各群の分布が同じであること」→「群間にばらつきがないこと」なので、分散比が有意であれば、「全体のばらつきに対して群間にばらつきがないとは言えない」、細かいことを言わなければ、群間にばらつきが認められ、各群の分布は同じでないことになる。
この例では、生産性はモジュールの階層によって違いがありそうということになる。

次の例は、同じ観測値を別の要因（観点）で分けたものである。

一元配置（対応なし）の例(2)
要因：所属
A社B社C社
生産性
(steps/人月)		 490
410
460
590
500
370
460
690
 750
480
500
650
770
 610
310
720
450
530
240
550
350
730
500
 
> summary(aov(productivity ~ company, data=sample1))
            Df Sum Sq Mean Sq F value Pr(>F)
company      2  68444   34222  1.7475 0.1998
Residuals   20 391677   19584               
>
測定値が同じであることを明確にする為に、データファイルを同じsample1にしている。'*'のマークとその説明が出力されていないので、所属別では群間に有意な差が無く、このデータでは所属によって生産性にばらつきがあるとは言えないことになる。

もう１つ例を考えた。

一元配置（対応なし）の例(3)
要因：仕様書のファイル形式
.doc.txt.xls.ppt
システムテスト
不具合数(/kstep)		 3
10
8
2
2
5
4
4
6
0
11 		11
7
1
3 		7
5
10
4
10
11 		12
5
10
14
11 
> sample2 <- read.table("anova02.dat", header=TRUE)
> summary(aov(error ~ format, data=sample2))
            Df Sum Sq Mean Sq F value  Pr(>F)  
format       3 113.58  37.861  3.1192 0.04672 *
Residuals   22 267.03  12.138                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
このデータでは、システムテストで見つかったバグ数は、仕様書のファイル形式によって違いがありそうということになる。

測定データを一元配置分散分析（対応なし）する時のポイントを整理する。

  • 測定値は正規分布することを仮定できるものであること
  • 差があることを調べる群は、3群以上あること （2群しか無いなら、その差をt検定するのと変わらない）
    つまり、要因（群の分け方）はその水準（標本が属する群を決めるもの、値やID）が3つ以上あるように選ぶこと
  • 各水準のデータ数は同じでなくても良い
  • 水準はなるべく定量的な数値でないこと、数値であれば、なるべく測定値の従属変数でないことが明確であること
    （要因と測定値に明確な関係があっても、群間のばらつきが十分に大きくないと有意差が検出できないため。また、水準（値の範囲）の取り方によって結果が変わってしまうため。例えば測定値と線形の関係にあるなら、無相関の検定をする方が良い）
最後の1点は、この記事を書くにあたって色々な例を考えてみた上での筆者の感想であり、そういうことを何かで読んだ記憶も無く、必ずしもそうではないかも知れない（温度の範囲などを水準に取る例も見かける）。

一元配置（対応あり）の例として、次のようなデータを考える。
各開発者がその4種類の開発プロセス（作業手順）を経験した時の、それぞれの開発プロセスにおける成績が、次のように得られたとする。

一元配置（対応あり）の例(1)
生産性
(steps/人月)		要因：開発プロセス
WaterfallSpiralIncrementalTDD
対
応
要
因
 開発者A490750450610
開発者B410480530780
開発者C460500240680
開発者D590650550880
開発者E500770350520
開発者F370610730600
開発者G460310500400
開発者H6907206401030
Waterfallは要求分析、設計、実装、テストをこの順で1回ずつ行うやり方、Spiralは設計以降の作業を何回かに分けて設計、実装、テストを繰り返して積み上げ式に開発する方式、Incrementalは要求分析も含めて全体を繰り返す、1サイクル毎に一通り動くものができる機能追加方式、TDDは要求分析や設計をテストプログラム作成に代え、以降の作業はそのプログラムがOKを返すことだけを目的に行う、結果良ければ全て良し方式のことである。
開発者の能力や向き不向きには個人差があることを前提とし、それを計算に入れて、開発プロセスは生産性に影響するかどうかを、分散分析で調べる。 
> sample3 <- read.table("anova03.dat", header=TRUE)
> summary(aov(productivity ~ process + Error(developer), data=sample3))

Error: developer
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  7 337872   48267               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)  
process    3 201184   67061  3.8348 0.02464 *
Residuals 21 367241   17488                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
"Error: developer"の部分は開発者の個人差による（開発者の違いを誤差要因とする）ばらつきに関するものであり、以降の計算はそれが除外されていることを示す。"Error: Within"の部分が、群内変動を誤差と見なして群間変動を検定するものである。
このデータでは、個人差の影響を差し引くと、開発プロセスは生産性に影響する要因と言えることになる。
もし、同じデータを「対応なし」として計算する（開発者の個人差を差し引かない）と、次のように、開発プロセス間に同じ有意差が出ない。 
> summary(aov(productivity ~ process, data=sample3))
            Df Sum Sq Mean Sq F value  Pr(>F)  
process      3 201184   67061   2.663 0.06729 .
Residuals   28 705112   25183                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
実際には各開発プロセスを短期間に経験することは難しく、開発者の経年変化や学習効果もあるので、こういうデータは取りにくいが、例えば開発者がソフトハウスのことであり、ソフトハウス内ではばらつきが無いと仮定すれば、各プロセスを同時に行うことも可能であるし、何より、筆者がこれ以上単純明快な例を思い付かないので、これで良しとする。

次の例は、会社による差（観測対象の個体差）はあるものとして、業務分野の違いが生産性に影響すると言えるかどうかを調べる。

一元配置（対応あり）の例(2)
生産性
(steps/人月)		要因：開発分野
組み込み汎用業務システムゲーム
対
応
要
因
 M社330470450320
N社290440370630
O社450550580750
P社320360470460
Q社370320430480
> sample4 <- read.table("anova04.dat", header=TRUE)
> summary(aov(performance ~ product + Error(company), data=sample4))

Error: company
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4 102420   25605               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)  
product    3  80080 26693.3  4.0332 0.03379 *
Residuals 12  79420  6618.3                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
summary(aov(performance ~ product, data=sample4))
            Df Sum Sq Mean Sq F value Pr(>F)
product      3  80080   26693  2.3487 0.1111
Residuals   16 181840   11365
1つ目のaov()が会社間に「対応あり」として計算した例、2つ目のaov()が「対応なし」として計算した例である。このデータだと、会社間の個体差を考慮に入れると、開発分野は生産性に影響する要因だったと言えることになる。

測定データを一元配置（対応あり）にする時のポイントを整理する。

  • 群内の各測定値が、どの個体から得られた標本であるか（または、誤差要因のどの水準に属する標本であるか）がわかること
  • 個体差によるばらつきはあるものとし、差し引いて考えるものであること
  • 個体に有意差があるかどうかを調べる必要は無いこと
  • 一部のデータが抜けていても計算可能
  • 同水準、同個体のデータは1つでなく複数あっても計算可能（但し全てのセルのデータ数が同じであることが望ましい）
つまり、「対応あり」とは、異なる水準間（上の表の列間）でどれとどれが同じ個体から採取された（一般化すると、同じ誤差要因を持つ）標本であるかの対応が取れることである。

二元配置（繰り返しなし）の例として、次のようなデータを考える。

二元配置（繰り返しなし）の例(1)
生産性
(steps/人月)		要因1：開発プロセス
WaterfallSpiralIncrementalTDD
要因2：
座席の形態
大部屋490700450560
自由席410430530730
パーティション460450240630
小部屋590720550830
自宅500600350470
開発プロセスの違いも作業空間（座席の形態）の違いも生産性に影響し、それらの影響が足し合わせれていると仮定して、それぞれの組み合わせについて1つずつデータを採取し、2要因について同時に、その仮定が正しそうかどうかを調べる。 
> sample5 <- read.table("anova05.dat", header=TRUE)
> summary(aov(productivity ~ process + workspace, data=sample5))
            Df Sum Sq Mean Sq F value  Pr(>F)  
process      3 141255   47085  4.8533 0.01951 *
workspace    4 121420   30355  3.1288 0.05589 .
Residuals   12 116420    9702                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
よって、このデータでは、開発プロセスの違いによる効果は有意水準0.05で有意差あり、作業空間の違いによる効果は有意差なし（有意水準が0.1だと有意差あり）となる。
ちなみに、開発プロセスと作業空間についてそれぞれ、同じデータを一元配置とみなして計算すると、次のように、同じ水準の有意差は出ない。 
> summary(aov(productivity ~ process, data=sample5))
            Df Sum Sq Mean Sq F value  Pr(>F)  
process      3 141255   47085  3.1675 0.05318 .
Residuals   16 237840   14865                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary(aov(productivity ~ workspace, data=sample5))
            Df Sum Sq Mean Sq F value Pr(>F)
workspace    4 121420   30355  1.7671  0.188
Residuals   15 257675   17178
このデータの特徴は、全体のばらつきを、2つの要因のばらつきと共通のばらつきの3つに分解するからこそ、明確になるのである。

次の例は、何らかの開発成果物の欠陥数を、作業期間と納入した曜日との組み合わせのそれぞれについて1つずつ、データを採取したとするものである。

二元配置（繰り返しなし）の例(2)
欠陥数
/kstep		 要因1：開発期間
１週間２週間１ヶ月
要因２:
曜日 		月曜1064
火曜696
水曜864
木曜494
金曜864
土曜1099
日曜N/A69
> sample6 <- read.table("anova06.dat", header=TRUE)
> summary(aov(error ~ time_limit + delivery, data=sample6))
            Df Sum Sq Mean Sq F value  Pr(>F)  
time_limit   2 14.360  7.1798  2.8898 0.09803 .
delivery     6 48.861  8.1435  3.2777 0.04224 *
Residuals   11 27.329  2.4845                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
1 observation deleted due to missingness
> summary(aov(error ~ delivery, data=sample6))
            Df Sum Sq Mean Sq F value  Pr(>F)  
delivery     6 51.217  8.5361  2.8213 0.05523 .
Residuals   13 39.333  3.0256                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
1 observation deleted due to missingness
よって、このデータでは、開発期間の違いによる差も考慮すると、納品日の曜日によって有意水準0.05で有意な差があるという結果になる。
2つ目のaov()は、上の例と同様、一元配置では同じ有意差が出ないことを確認したものである。

二元配置（繰り返しなし）にする時のポイントを整理する。

  • 2つの要因の全ての水準の組み合わせについて、データは1つ
  • どちらの要因も何らか影響していると思われ、それぞれの影響を分離して調べようとしていること
  • 一部のデータが抜けていても計算可能
  • 一元配置（対応あり、繰り返しなし）の代用として行うことも可能
表の形が同じことからもわかるように、一元配置（対応あり）と二元配置は計算方法が同じであり、一元配置の対応を決める要因についてのF検定を省くかどうかが異なるだけである。
Excelの分析ツールに「対応のある一元配置」が無いのは、二元配置で代用できるからだろう。

二元配置（繰り返しあり、対応なし）は、表の形としては、繰り返しのない二元配置の各セルに複数のデータが含まれるだけの違いなので、同じような要因のペアが使える。

二元配置（繰り返しあり、対応なし）の例(1)
生産性
(steps/人月)		要因1：開発プロセス
WaterfallSpiralIncrementalTDD
要因２：
座席の形態
大部屋 570
660
450		 710
480
570		 570
450
620		 520
540
740
自由席 570
500
270		 590
460
670		 510
450
410		 620
350
350
パーティション 300
480
580		 670
500
710		 610
540
540		 430
650
600
小部屋 510
530
530		 660
620
510		 310
550
780		 760
400
590
自宅 510
560
530		 570
750
590		 400
380
520		 600
610
580
繰り返しがある（同じセルに複数のデータがある）と、行の効果、列の効果に加えて、交互作用の効果を分離して計算することが可能で、交互作用を分けるか分けないかで計算結果が変わる。
交互作用とは、2要因の水準の特定の組み合わせにだけに影響する効果のことで、例えばその行とその列は大体高い値なのにそのセルだけやたら低い値が多いということを起こす要因である。
aov()で交互作用を分離して分散分析するには、引数において2つの要因を'*'で繋ぐ。 
> sample7 <- read.table("anova07.dat", header=TRUE)
> summary(aov(productivity ~ process * workspace, data=sample7));
                  Df Sum Sq Mean Sq F value  Pr(>F)  
process            3  98952   32984  2.4676 0.07599 .
workspace          4  65823   16456  1.2311 0.31308  
process:workspace 12  69657    5805  0.4343 0.93963  
Residuals         40 534667   13367                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
"process:workspace"の行が交互作用に関する行で、このデータでは有意差は出なかった。また、交互作用を分離すると、開発プロセスにも作業空間にも有意水準0.05の有意差はなしである。
交互作用を分離せずに計算するには、aov()の引数において、2つの要因を'+'で繋ぐ。 
> summary(aov(productivity ~ process + workspace, data=sample7));
            Df Sum Sq Mean Sq F value  Pr(>F)  
process      3  98952   32984  2.8382 0.04686 *
workspace    4  65823   16456  1.4160 0.24170  
Residuals   52 604323   11622                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
このデータでは、交互作用を分離しなければ、開発プロセスには有意差が認められる（交互作用の効果を分離しない方が有意差が出る）ことがわかる。
ちなみに、繰り返しがないと、'+'でも'*'でも結果は変わらない。

もう１つ例を作る。

二元配置（繰り返しあり、対応なし）の例(2)
結合テストエラー数
/kstep		 要因1：仕様書の配布形態
WordExcelPDFHTML
要因2:
設計書の
フォーマット 		Word 8
6
2
1
7		 5
7
3
9
10		 4
11
5
6
3		 11
7
14
7
9
Excel 11
5
8
13
8		 4
10
6
3
6		 5
9
13
6
10		 17
12
11
5
10
Text 5
6
6
10
11		 6
12
7
9
8		 5
12
8
9
6		 9
8
13
11
3
HTML 12
14
7
13
12		 10
15
6
13
13		 12
7
5
13
7		 4
10
3
8
9		 
> sample8 <- read.table("anova08.dat", header=TRUE)
> summary(aov(error ~ given_spec * design_desc, data=sample8))
                       Df Sum Sq Mean Sq F value  Pr(>F)  
given_spec              3   86.5 28.8333  2.9254 0.04043 *
design_desc             3   17.1  5.7000  0.5783 0.63137  
given_spec:design_desc  9  198.4 22.0444  2.2366 0.03056 *
Residuals              64  630.8  9.8562                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary(aov(error ~ given_spec + design_desc, data=sample8))
            Df Sum Sq Mean Sq F value  Pr(>F)  
given_spec   3   86.5  28.833  2.5384 0.06314 .
design_desc  3   17.1   5.700  0.5018 0.68221  
Residuals   73  829.2  11.359                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1つ目のaov()の結果より、このデータでは仕様書のフォーマットが欠陥数に影響しており、また設計書のフォーマットとの交互作用もある結果となった。 2つ目のaov()の結果は、交互作用を分離しないと、仕様書のフォーマットについても同じ有意差は出ない（１つ前のデータ例とは逆に、交互作用を分離する方が有意差が出る）ことを示している。

二元配置（繰り返しあり、対応なし）にする時のポイントを整理する。

  • 2つの要因の全ての水準の組み合わせについて、複数のデータがある
  • それにより、2つの要因の交互作用が混入する
  • 交互作用によるばらつきを分離して計算するかどうかは場合による
  • 一元配置（対応あり、繰り返しあり）の代用として行うことも可能
一元配置（対応あり）を代用する時に交互作用をどう扱うかが問題になるが、誤差要因と何かの交互作用というのはやはり誤差要因だと考えるなら、交互作用を分離して、誤差要因も交互作用も検定から除外する（aov()なら'*'を使った上で検定中の要因そのもの以外の行を無視する）のが好ましいと思う。

「二元配置（１要因対応あり）」は、「対応のある要因と対応のない要因の二元配置」と書かれることが多いようだが、そう書かれると余計にややこしく見えるのは、筆者だけだろうか。
「対応あり」とは、一元配置の場合と同様、繰り返しのある（セル内に複数ある）データのそれぞれが、他の水準（セル）のどのデータに対応するかがわかる、という意味であり、対応を取るからには、個体差の影響を取り除いて計算する必要があると考えていることになる。
１要因についてのみ対応が取れる状況として、次のような例を考える。

二元配置（繰り返しあり、１要因対応あり）の例(1)
生産性
(steps/人月)		要因1：開発プロセス
WaterfallSpiralIncrementalTDD
要因２：
座席の形態
+ 人 		作業場開発者
大部屋 A630660570340
B530430430400
C390640520530
パーティション D400760500680
E550670290890
F650620280820
小部屋 G470680760410
H530580670930
I850980790980
自宅 J610850530310
K400640230400
L530430470330
開発者間には能力の個人差があるのは普通だから、開発者の違いは、標本の対応を取ってその影響を除外して計算するにふさわしい、誤差要因であろう。
前記のように同じ人が4つの開発プロセスを経験するというのは現実的に多少無理があるが、それはなされたとする。
全ての開発者がこれらの全ての作業空間を経験するのは現実的に不可能であろう。と思うので、そういう席替えはなされず、作業場毎に別の人を選んでデータを採取したとする。
つまり、開発プロセスについては、水準間でどれとどれが同じ人のデータであるかの対応があり、作業空間に関しては、水準間でそういう対応がない、というデータである。
> sample9 <- read.table("anova09.dat", header=TRUE)
> summary(aov(productivity ~ process * workspace + Error(developer), data=sample9))

Error: developer
          Df Sum Sq Mean Sq F value  Pr(>F)  
workspace  3 424492  141497  3.8251 0.05734 .
Residuals  8 295933   36992                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Error: Within
                  Df Sum Sq Mean Sq F value  Pr(>F)  
process            3 163692   54564  3.0262 0.04914 *
process:workspace  9 391275   43475  2.4112 0.04123 *
Residuals         24 432733   18031                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
このように、開発プロセスの違いによる影響には有意差があり、また開発プロセスと作業空間の間には何らかの交互作用があるという結果になった。
もし同じデータを対応なしの二元配置として計算すると、次のように全く異なる結果になる。 
> summary(aov(productivity ~ process * workspace, data=sample9))
                  Df Sum Sq Mean Sq F value   Pr(>F)   
process            3 163692   54564  2.3962 0.086440 . 
workspace          3 424492  141497  6.2140 0.001898 **
process:workspace  9 391275   43475  1.9092 0.086483 . 
Residuals         32 728667   22771                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
これは、データをよく見るとわかるのだが、小部屋に飛び抜けて成績が良い人が居ることに起因しており、作業空間の違いによる効果だと言うにはかなり不適切であろう。

１要因にのみ対応がある例としてもう１つ、次のようなデータを考える。
開発者A〜Fで構成されるグループが、6つの分野の業務を行った時の生産性を測ったという例である。

二元配置（繰り返しあり、１要因対応あり）の例(2)
生産性
(steps/人月)		 要因1：性別 + 人
性別男性女性
開発者ABCDEF
要因２：
業務分野 		組み込み（アプリ） 470720550610580550
業務システム 490530530570310380
生産システム 530520570600520490
PCアプリ 530310460420690650
Javaアプリ 600530500450790590
ゲーム機ソフト 450530400530350530
上の例では対応を決める要因を縦に並べたが、今度は横に並べた。
開発者の能力には当然個人差がある。 全開発者が6つの分野にて仕事をしたので、分野間ではどれとどれが同じ人のデータであるかの対応が取れる。
1人の開発者は男性と女性の両方を経験できないので、男性群と女性群のデータには対応が取れない。 
> sample10 <- read.table("anova10.dat", header=TRUE)
> summary(aov(productivity ~ sex * area - Error(developer), data=sample10))

Error: developer
          Df Sum Sq Mean Sq F value  Pr(>F)  
sex        1 4225.0  4225.0  9.6266 0.03613 *
Residuals  4 1755.6   438.9                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
area       5  77314   15463  1.2919 0.3066
sex:area   5  51892   10378  0.8671 0.5202
Residuals 20 239378   11969               
> summary(aov(productivity ~ sex * area, data=sample10))
            Df Sum Sq Mean Sq F value Pr(>F)
sex          1   4225    4225  0.4205 0.5228
area         5  77314   15463  1.5390 0.2153
sex:area     5  51892   10378  1.0330 0.4209
Residuals   24 241133   10047
1つ目のaov()の出力より、このデータでは、性別による影響に有意差が見られるという結果になる。
2つ目のaov()の出力は、2要因とも対応なしとして計算すると、どの要因、交互作用にも有意差が見られないことを示すものである。

二元配置（１要因のみ対応あり）にする時のポイントを整理する。

  • 標本には個体差があることを前提とする
  • 1つの要因については、複数の水準に同じ個体からの標本があり、水準間でデータの対応が取れる
    例えば、同じ被験者から、各水準の条件下でデータが取られている
  • もう1つの要因については、1つの個体からの標本は1つの水準にしかなく、水準間でデータに対応がない
    例えば、水準を決める条件は同時に発生するので、水準毎に異なる被験者を選ぶ
もちろん、対応を決める要因は、個人や個体に限らず、標本に影響があることが既にわかっている条件であれば何でも良い。

二元配置（２要因対応あり）の標本のデータ構造は、１要因のみ対応ありのものよりも単純である。

二元配置（繰り返しあり、２要因対応あり）の例(1)
欠陥数
(steps/人月)		要因1：開発プロセス
WaterfallSpiralIncrementalTDD
要因２：
座席の形態
+ フェーズ 		作業空間開発フェーズ
大部屋 仕様作成6673
全体設計4787
個別設計68105
テスト設計3264
実装6588
小部屋 仕様作成2635
全体設計5746
個別設計63104
テスト設計7564
実装12567
自宅 仕様作成5494
全体設計9296
個別設計7997
テスト設計5332
実装491011
要するに３次元なのである。ここでは表を立体的に書けないので、対応を決める「開発フェーズ」については列方向に展開している。
開発フェーズが違うと仕事の性質が全く違うので、フェーズの違いは当然ミスの数に影響する。
どんな開発プロセスにも、これくらいのフェーズは存在するので、各フェーズのデータはあり得る。また、どんな作業空間でも、一通りの開発をすれば全てのフェーズを経るので、各フェーズのデータが得られる。従って、異なる開発プロセスのデータ間でも異なる作業空間のデータ間でも、フェーズの違いによる対応が取れる。 
> sample11 <- read.table("anova11.dat", header=TRUE)
> summary(aov(error ~ process * workspace - Error(phase/(process*workspace)), data=sample11))

Error: phase
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4 94.733  23.683               

Error: phase:process
          Df Sum Sq Mean Sq F value  Pr(>F)  
process    3 30.850 10.2833  3.8482 0.03852 *
Residuals 12 32.067  2.6722                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Error: phase:workspace
          Df  Sum Sq Mean Sq F value Pr(>F)
workspace  2  4.9333  2.4667  0.6251 0.5594
Residuals  8 31.5667  3.9458               

Error: phase:process:workspace
                  Df Sum Sq Mean Sq F value Pr(>F)
process:workspace  6  17.20  2.8667  0.5148 0.7912
Residuals         24 133.63  5.5681
２要因対応ありのデータ構造は単純だが、その計算は、コンピューター任せでも複雑である。対応を決める要因の違いによる影響が、要因１との交互作用、要因２との交互作用、要因１と要因２と対応要因との交互作用と多岐に渡って分離して計算されるからである。
上記のaov()の構造式のError()の部分は、その構造がわかりやすいように
  • Error(phase + phase/process + phase/workspace)
  • Error(phase + phase:process + phase:workspace)
  • Error(phase + phase:process + phase:workspace + phase:process:workspace)
といった形で書かれることも多いようだが、ここではシンプルで入力ミスも避けられる
  • Error(phase / (process * workspace))
を採用している。
上の計算結果では開発プロセス間に有意差が出たが、同じデータを対応なしとして計算すると、次のようにどこにも有意差が出ない。 
> summary(aov(error ~ process * workspace, data=sample11))
                  Df  Sum Sq Mean Sq F value Pr(>F)
process            3  30.850 10.2833  1.6904 0.1816
workspace          2   4.933  2.4667  0.4055 0.6689
 process:workspace  6  17.200  2.8667  0.4712 0.8262
Residuals         48 292.000  6.0833

もう１つ例を考える。新規に参入した人がそこの設計業務をマスターするのに何ヶ月かかったかというデータが、次のように得られたとする。

二元配置（繰り返しあり、２要因対応あり）の例(2)
定着期間
（ヶ月）		要因1：設計手法
構造化設計データ指向オブジェクト指向コンポーネント指向
要因２：
開発内容
+ 性別 		開発内容性別
プラットフォーム 男性 8 9 711
女性 3 6 7 7
ミドルウェア 男性 2 4 3 6
女性 9 5 9 8
フレームワーク 男性10 910 9
女性 910 511
アプリケーション 男性11 9 5 5
女性1210 1 5
男性と女性とでは元々差があるはずだという前提で、対応ありの分散分析を行う。 
> sample12 <- read.table("anova12.dat", header=TRUE)
> summary(aov(takes_month ~ design_policy * target - Error(gender/(design_policy*target)), data=sample12))

Error: gender
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals  1 0.03125 0.03125               

Error: gender:design_policy
              Df  Sum Sq Mean Sq F value Pr(>F)  
design_policy  3 23.3438  7.7813  14.647 0.0269 *
Residuals      3  1.5938  0.5313                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Error: gender:target
          Df Sum Sq Mean Sq F value Pr(>F)
target     3 45.844  15.281  0.8886 0.5375
Residuals  3 51.594  17.198               

Error: gender:design_policy:target
                     Df Sum Sq Mean Sq F value Pr(>F)
design_policy:target  9 95.531 10.6146  2.3142 0.1136
Residuals             9 41.281  4.5868               
> summary(aov(takes_month ~ design_policy, data=sample12))
              Df  Sum Sq Mean Sq F value Pr(>F)
design_policy  3  23.344  7.7812  0.9237 0.4422
Residuals     28 235.875  8.4241
このデータでは、男女の違いによる対応ありとすると、設計手法の違いによる有意差が見られ、対応なしとすると、設計手法にも開発内容にもその違いによる差が見られないという結果になった。

統計学復習メモ18: 少数のベルヌーイ試行による発生確率の区間推定

問：ある機械が20個の部品を作ったら、不良品が2個出てきた。この機械が1個作る毎に不良品を発生させる確率は一定だとすると、不良品発生率は信頼度95%の信頼区間として何%〜何%の間と推測できるか。

観測誤差を考えなければ、すなわち点推定量としては、推定確率はもちろん2/20=0.1である。しかし、標本はたった20個である。これをもって、だから2000個作ったら200個の不良品が出ることが予想される、と言うのはあまりにも乱暴だし、直感的にナンセンスであろう。サンプルが少なすぎて信憑性が低く、説得力が無いのである。
信頼性を高めるにはサンプル数を増やせば良いが、例えば以前のエントリーに書いたような計算方法では、推定確率が0.1前後の場合に誤差を±0.05以内にするには、(1.96/0.05)2*0.1*0.9＝約140個のサンプルが必要になる。サンプル数を先に決められない場合や、サンプルはそれしかないという場合は、推定発生率はこの範囲に入る、という信頼区間を求めるしか無い。

観測されたベルヌーイ試行の成功回数（発生回数）からの元の成功率（発生率）の区間推定については、以前のエントリーに書いたように信頼区間の公式がある。
I \in \hat{p}\pm Z\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} ...(1)
（p^は点推定量、Zは標準正規分布N(0,1)に従う値、nはベルヌーイ試行の標本数）

従って、p^=0.1、Z=1.96（右側2.5%点）、n=20とすると、95%信頼区間は0.1±1.96√0.0045 = 0.1±0.131、すなわち-3.1%〜23.1である。（完）

イエース、公式一発、ビバ統計学！
...ちょっと待て、確率なのにマイナスの値？何を間違えたのだろう？

発生率が一定で、１回の試行の結果は発生するかしないかの２通りしか無い試行（ベルヌーイ試行）を複数回行った時の発生回数は二項分布に従う。そして試行回数が十分大きければ二項分布は正規分布で近似できる。試行回数をn、発生率をpとすると、発生回数は平均np、分散np(1-p)の正規分布に（近似的に）従う。観測される発生率（推定確率）は発生回数を試行回数で割ったものだから、平均p、分散p(1-p)/nの正規分布に（近似的に）従う。母分散が既知の場合、正規分布に従う標本の推定値は毋分散÷標本数を分散とする正規分布に従うとするのが区間推定の定跡である（母分散が未知の時は標本分散を使うのでｔ分布に従うとする）から、信頼区間Iは
I \in \hat{\theta}\pm Z\sqrt{\frac{\sigma^2}{n}}
（θ^は点推定量、Zは標準正規分布のZ値、このnは正規分布に従う標本の数なので今回は1）
で与えられる。冒頭の問題については、θ^=p、σ2=p(1-p)/nなので、信頼区間はやはり(1)の公式の通りである。

という訳で、二項分布を正規分布に近似しているのが問題のようである。
次のグラフは、n=20として、p=0.5とp=0.1の二項分布とN(np, p(1-p)/n)の正規分布を比較したものである。

青線が二項分布、赤線が正規分布である。p=0.5だとかなり近いが、p=0.1だと差がある。比率が0.5から遠いほど、二項分布の正規分布への近似は悪くなるのである。
大体、正規分布に近似するということは0より小さい確率や1より大きい確率を認めるということであり、n=20でp=0.1の場合、正規分布の左側6.8%は負の値に入り、無効になってしまうのである。この6.8%を全てp=0の所に加算すれば（信頼区間の下限が負の値なら0に補正すれば）いいだけの話とも言えなくはないが、pの推定値が0に近づくと急激にp=0の確率が上昇することになるので、ちょっと強引であろう。

上記の比率の信頼区間の公式(1)は、np<5またはn(1-p)<5の場合は不適、すなわちベルヌーイ試行の成功回数（発生回数）も失敗回数（発生しなかった回数）も5以上でないと使えないと言われているらしい。視聴率調査や投票率予想なら最低5人が観て／投票していないと、この公式を使うのは適当でないのである。それに基づくと、冒頭の問題もnp=2なのでこの公式は使えない。

ではこの場合はどのようにして信頼区間を求めれば良いのだろうか。それを解決してくれるのがClopper-Pearson methodによる通称Exact Confidence Interval(Exact CI)である。
区間推定の原点に立ち返ると、信頼区間とは、ある推定値がある信頼度で収まる値の範囲であり、まともに計算するなら事前確率、事後確率を考慮して計算すべきものであるが、それに代えて、推定の元になった観測値及びそれより稀な値が発生する確率が（1-信頼度）以下であるような推定値の範囲を信頼区間として計算する、というのを二項分布のパラメーター推定について行うのがClopper-Pearson methodである。冒頭の問題なら、95%の確率でpが0.1の前後のどこまでの範囲に含まれるかを求める代わりに、pが0.1よりどこまで小さければ不良品が2個以上出る確率が2.5%になるか、pが0.1よりどこまで大きければ不良品が2個以下である確率が2.5%になるか、を求めるのである。
という訳で、ベルヌーイ試行の確率の推定値の信頼区間の下限Plbと上限Pubを束縛する式は次のようになる。（αは1-信頼度）
\sum_{x=0}^{k}\pmatrix{n \cr x}P_{ub}^x (1-P_{ub})^{n-x}=\frac{\alpha}{2}

\sum_{x=k}^{n}\pmatrix{n \cr x}P_{lb}^x (1-P_{lb})^{n-x}=\frac{\alpha}{2}
これを計算したものがExact CIで、BetaCDF(x,a,b)をベータ分布の累積密度関数として、その逆関数を用いて、
P_{ub}=1-BetaCDF^{-1}\left(\frac{a}{2},n-k,k+1\right)
P_{lb}=1-BetaCDF^{-1}\left(1-\frac{a}{2},n-k+1,k\right)
と書けるらしい。

Maximaにはベータ分布の累積密度関数の逆関数（β値を求める関数）がquantile_beta()という名前で存在するので、これを使って信頼区間の上限を計算できる。
(%i1) load(distrib);
(%i2) Pub(n,k,a) := 1-quantile_beta(a/2, n-k, k+1);
(%i3) Plb(n,k,a) := 1-quantile_beta(1-a/2, n-k+1, k);
(%i4) Pub(20,2,0.05);
(%o4) 0.317
(%i5) Plb(20,2,0.05);
(%o5) .012349

答：1.23%から31.7%の間

参考文献：Understanding Binomial Confidence Intervals

統計学復習メモ17: 重回帰分析の計算結果の検定

以前に回帰直線を求めることと主成分分析は似ていると書いたが、ちょっと無理があり過ぎたかも知れない。なぜか主成分分析が因子分析より知名度が低いのが気に入らなくて、メジャーなものと結びつけたくて書いたのだが、やっぱり回帰分析と主成分分析は全然違う。回帰分析は目的変数を説明変数で説明しようとするものであり、主成分分析はデータの傾向を見るものである。通常、回帰直線を求める際に最小にする誤差は目的変数軸上の誤差であるが、第一主成分を求める為に最小にする誤差は全変数が構成する空間での誤差である。
筆者は、主成分分析は多変量解析の基本かつ最も重要であると考えており、学生時代には趣味のプログラミングのネタにするほど惚れ込んだが、因子分析はあまり好みでないのである。暴走ついでにさらに暴走すると、筆者の中では、

A群B群C群
因子分析主成分分析回帰分析
LinuxFreeBSDOpenBSD
RubyPythonPerl
WindowsMacUNIX
SleipnirLunascapeFirefox
秀丸サクラエディタEmacs
のような関係があり、因子分析は軟派で軽薄でチャラいのである。


閑話休題。今回は重回帰分析を復習する。一次回帰直線が、xを説明変数、yを目的変数とし、直線y^=a x+bの係数aと切片bをデータから求めたものだとすると、重回帰分析はそのxを複数次元にしたもの、
¥hat{y}=a_1x_1+a_2x_2+¥cdots+a_mx_m+b
の係数aiとy切片bを求めるものである。
便宜上、上のbをa0と書き、
¥hat{y}=XA=¥pmatrix{1 & x_1 & ¥cdots & x_m}¥pmatrix{a_0 ¥cr a_1 ¥cr ¥vdots ¥cr a_m}
と行列の積で表すことにする。

n個のデータの組(yi, x1i, x2i, x3i, ...)(1<=i<=m)が得られた時、誤差をeiとすると、

である。以下、これをY=XA+Eと書く。
yの誤差eiの２乗和が最小になるようにAを求めるとすると、
¥sum_{i=1}^{n}e_i^2=E^TE=(Y-XA)^T(Y-XA)
であり、これはaiの２次関数なので、これを偏微分して0になるAが答であるから、
¥frac{¥partial}{¥partial A}(E^TE)=¥frac{¥partial}{¥partial A}(Y^TY-2(XA)^TY-(XA)^TXA)=-2X^TY+2X^TXA=0
なので
A=(X^TX)^{-1}X^TY
という計算で求められることがわかる。(XTX)が正則でないと逆行列が存在しないという問題はあるが、異なるデータの個数が説明変数の数より少ないとか、中身が0だらけのぜろりとしたデータでない限りは正則だろうから、とりあえず気にしないことにする。

試しに次のサンプルデータのバグ数について計算してみる。

部品バグ数開発規模関係者数疲労度
A70.935
B152.118
C101.249
D71.217
E60.359
F50.693
G70.922
H50.363
I60.358
J70.634
これを"bugs.dat"というタブ区切りのテキストファイルにして、Maximaで 
dat:read_matrix("bugs.dat");
n:10$ /*標本数*/
p:3$ /*説明変数の数*/
Y:submatrix(dat,2,3,4); /*1列目を取り出す（2-4列目を除去）*/
X1:submatrix(dat,1); /*2列目以降を取り出す（1列目を除去）*/
X:map(lambda([x],append([1],x)),X1); /*X1の左端に要素が全て1の列を追加*/
A:invert(transpose(X).X).transpose(X).Y; /*回帰係数*/
fpprintprec:4$
float(A);
を実行すると、次の結果が得られる。
a02.1
a14.612
a20.02896
a30.2436
バグ数=2.1 + 4.612*開発規模 + 0.02896*関係者数 + 0.2436*疲労度、という意味である。
この結果がどれくらい元データに合ってるかを、Yの値で見てみる。 
/*元データ、回帰式による推定値、誤差*/
[Y, X.A, Y-X.A];
元データYYの推定値Y^Y-Y^
77.556-.5558
1513.761.237
109.943.05743
79.369-2.369
65.8210.179
55.859-.8589
76.7960.204
54.388.6115
65.577.4226
75.9291.071

YとY^が十分近いので、これでも満足なのだが、統計学にはこの回帰分析結果をきちんと検定する方法が存在するようである。
その前に一応、よくある指標を計算する。以下、データ数をn、説明変数の数をpとする。

load(descriptive)
load(distrib)
n:10$
p:3$
/*重相関*/
R:cor(addcol(Y,X.A)),numer;
/*決定係数（寄与率）*/
R^2;
/*自由度補正済み決定係数*/
1-(1-R[1,2]^2)*(n-1)/(n-p-1);
重相関0.9361
決定係数0.8762
自由度修正決定係数0.8143
相関係数の類は、値の解釈が難しいので、相関係数同士を比較する以外には使いづらい。一応、自由度修正済み決定係数が0.8以上というのがよく使われる基準であるらしい。
これは置いといて、まず、回帰式全体の有効性を検定する。これには、回帰分散÷残差分散、
¥frac{V_{¥hat{y}}}{V_¥varepsilon}=¥frac{¥frac{1}{p}¥sum_{i=1}^{n}(¥hat{y}_i-¥bar{y})^2}{¥frac{1}{n-p-1}¥sum_{i=1}^{n}(¥hat{y}_i-y_i)^2}
という値が自由度p,n-p-1のF分布に従う、というのが使えるらしい。YとY^の誤差の分散が十分に小さいかどうかを調べようとするものと言えるだろう。 
E: Y-X.A$
/*回帰分散*/
Y_hat: X.A-mean(Y)[1]$
Vy: transpose(Y_hat).Y_hat/p,numer;
/*残差分散*/
Ve: transpose(E).E/(n-p-1);
/*回帰方程式のF値*/
f:Vy/Ve;
/*そのF値以上の発生確率*/
1-cdf_f(f,p,n-p-1);
結果
F値14.16
P(F(3,6)>14.16)0.03951
有意水準を0.05とすると、回帰分散＝残差分散という仮説は棄却されるので、この回帰式は有効ということになる。

次に、回帰係数aiそれぞれの有効性を検定する。これには、回帰係数÷標準誤差＝
¥frac{a_i}{¥sqrt{s^{ii}V_{¥varepsilon}}}
（sijは共変動行列（＝ｎ×共分散行列）の逆行列の要素）
という値が自由度n-p-1のｔ分布に従う、というのを使う。定数項a0については、siiの部分は次のような式になるが、同じようにｔ値が求められる。
s^{00}=¥frac{1}{n}+¥sum_{i=1}^{p}¥sum_{j=1}^{p}¥bar{X_i}¥bar{X_j}s^{ij}

/*共変動行列の逆行列*/
Sinv: invert(cov(X1)*n);
/*標準誤差（係数項）*/
Se[i]:=sqrt(Sinv[i,i]*Ve);
/*標準誤差（定数項）*/
Se[0]:sqrt((1/n + sum(sum(mean(X1)[i]*mean(X1)[j]*Sinv[i,j],j,1,p),i,1,p))*Ve);
for i:0 thru 3 do print(Se[i])$
iaiの標準誤差
01.863
11.023
20.2303
30.1661
/*回帰係数のt値、|t|値がそれより大きい確率*/
t[i]:=A[i+1,1]/Se[i]$
for i:0 thru 3 do print(t[i], 2*(1-cdf_student_t(t[i],n-p-1)))$
iｔ値確率
01.1280.3025
14.5070.004074
20.12570.904
31.4670.1928
有意水準を0.05とすると、a1についてのみ、本当は0である（係数に意味が無い）という仮説が棄却され、有効ということになる。

統計学復習メモ16: 帰無仮説と棄却域

統計学の検定というのは、「帰無仮説」と「対立仮説」を設定して、「帰無仮説」を「棄却」できれば「対立仮説」が「有意」であるとする、というものである。
...というように言葉で書かれるとそんなに難しくなさそうなのだが、ここまでで既に「帰無仮説」だの「有意」だの独特の用語が出てくる上、帰無仮説が棄却できない時は「対立仮説が有意であるとは言えない」だの「差が無いという仮説は棄却されない」だの「差があるとは言えない」だの回りくどい表現が出てきて、混乱する。
棄却して「無」に「帰」するための仮説だから「帰無仮説」と言うんだよ、なんて説明は大して理解の助けにならない。むしろ、それで何か解った気にさせて思考を停止させる、本質の理解を邪魔する雑音にすら思える。
なぜ「差が無いとは言えない」ではないのか。そういう疑問を持ってもう一度理解しようとすると、「第１種の過誤」や「第２種の過誤」や「有意水準」や「危険率」や「検出力」が出てきて、阻止されてしまう。命題の逆が偽（矛盾する）なら命題が真となるような単純な代物ではないのである。

筆者は大学で「検定」を習って15年以上、統計学を復習して1年以上、まだ混乱している。そのため、自己流に整理して、悪あがきしてみる。


仮説として、ある統計量が、下の図の緑色の線で示されるような確率分布に従うと仮定されているとし、標本から得られた統計量がtであるとする。横軸がt、縦軸が発生確率である。

上記の仮説が正しいとすると、t=0.7（青い点の辺り）やt=1.5（紫色の点の辺り）はまだそれなりに発生確率が高いが、t=2.2（赤い点の辺り）は発生確率が低い。従って、t=2.2となるような標本が得られたら、それはかなり稀な事が起こったという事になる。その時、そうではなくて仮説が誤っているのだと考えるのが統計学における検定の考え方である。

そのために、どれくらい稀だと仮説が誤っているとする（仮説を棄却する）かを決める。その確率が有意水準である。何故だかわからないが、0.05や0.01がよく使われるようである。ここでは有意水準を0.05(5%)とする。
次に、統計量がどの範囲だと仮説を棄却するか（棄却域）を、その範囲の発生確率が有意水準に等しくなるように決める。例えば、上図の緑色で塗った部分は、0から遠い、発生確率が全体の5%の範囲である。すなわち、標本のtがこの範囲（大体t<2とt>2）に入れば、仮説を棄却する、という設定である。

注意すべきは、仮説を棄却できなかったからといって仮説が正しいとは言えないことである。仮説が正しいとすると稀にしか起こらない事が起こったかどうか、を調べているだけなので、わかるのは、仮説が誤っているか、何とも言えないか、のどちらかである。調べている仮説を帰無仮説、その反対の仮説を対立仮説として言い換えると、わかるのは、帰無仮説が誤っていて対立仮説が有意であるか、帰無仮説が誤っているとは言えず対立仮説については何とも言えないか、のどちらかである。

上記の設定では、もし標本から得られたtが青や紫の位置の値なら、仮説が棄却できず（検定失敗）、いわゆる玉虫色の決着となるが、もしtが赤い位置の値なら、仮説が間違っているだろうと言える。


以上で検定の計算は可能だが、設定した検定方法が妥当かどうか、また検定結果の意味を考えるためには、他に考える事がある。
帰無仮説にて統計量の確率分布を仮定するということは、標本の統計量にもばらつきがあることを仮定している。もし標本の統計量にばらつきが無いとするなら、棄却域に入るかどうかという以前に、ばらつきがあるとする帰無仮説が誤りになってしまう。
例えばもし、標本の統計量が帰無仮説の通りでないのに、平均的には棄却域に入らない場合は、次のような感じになる。

赤い点が統計量の平均、紫色の線が分布である（見やすさのため、縦に縮めている）。ピンクで塗られている部分が、帰無仮説（緑線）の棄却域に入る部分である。
もし帰無仮説が正しい場合は、紫の線は緑の線に一致し、ピンクの領域は有意水準と同じ5%となる。そのため、帰無仮説が正しくても、5%の確率で帰無仮説が棄却されてしまう（第一種の誤り）。帰無仮説は棄却するために設定する（誤りだろうと思うものを帰無仮説に採る）ものだが、棄却できたと喜んでも、誤って棄却した可能性が残るのである。有意水準は、この第一種の誤りを犯す確率という意味で、危険率とも呼ばれる。

上の図では、緑の線と紫の線はずれており、帰無仮説は誤りである。しかし、ピンクの領域、すなわち帰無仮説が棄却される確率は12.5%しかない。「検出力」が12.5%しか無い、と言う。
次の図は、統計量の平均が棄却域に入るくらいにずれている例である。

ピンクで塗っていないが、ここまでずれると検出力は70%に達する。しかし、逆に言うと、ここまで帰無仮説が誤っていても、30%の確率で棄却されないのである（第二種の誤り）。有意水準を下げて第一種の誤りを犯す確率（緑の領域）を小さくすると、第二種の誤りを犯す確率（藤色の領域）が大きくなる、という関係があることがわかる。


次に、帰無仮説と棄却域の取り方について考える。
通常、帰無仮説は「差が無い」方に取られる。これは、統計量が0から遠い方が「差がある」ことが多いというか、0から遠い方が「差がある」ように統計量を定める方が便利であることが多いからである。例えば標準正規分布やｔ分布やカイ２乗分布に従う統計量は0から遠いと平均と差がある。
そのため、もし帰無仮説を「差がある」方に取ると、統計量の棄却域が0に近い方になり、範囲が狭くなる。標準正規分布やｔ分布のような山形の分布だと、

このように棄却域が狭くなり、統計量が繊細で、計算上不便である。

それよりも、実際の検定統計量は和や積を使って計算されるため、「差がある」は単純ではなく、棄却するのは難しいのである。例えば、和で効いてくる平均値の場合、実際には標本に帰無仮説と全然違う大きなばらつきがあっても、相殺されて帰無仮説の平均値に近くなってしまうことがあり得るし、積で効いてくる場合は標本のどれか1つの確率が（帰無仮説の前提から外れて）0に近いために結果として統計量が0に近くなってしまうことがあり得る。「差がある」を棄却するには色々な可能性を考えないといけないが、統計量が0から遠ければそれだけで「差が無い」を棄却できるので、やりやすいのである。

そのため、上の例のように統計量が標準正規分布やｔ分布などの0中心の左右対称の確率分布に従うと仮定する場合は、棄却域を0から遠い左右両端に定めることが多い（両側検定という）。
カイ２乗分布に従う統計量のように、正の値しか取らない場合は、0から遠い方の（+∞までの）1つの領域を棄却域とする（片側検定という）場合が多いが、それでも、どこかに山がある場合は0に近い方の裾と0から遠い方の裾の両側に棄却域を設ける場合もあるようである。

次の図は、カイ２乗検定（片側検定）の棄却域の例である。

緑の線が、統計量が自由度2のカイ２乗分布に従うという帰無仮説であり、緑で塗られた領域が有意水準（αと書かれることが多い）5%の棄却域である。赤い点のように、標本から得られた統計量がここに入れば、緑の線を棄却する。
実際には自由度4のカイ２乗分布に従う場合は、帰無仮説と標本は次のような関係になる。

藤色の部分が第二種の誤りを犯す範囲（βと書かれることが多い）で、確率は約0.8であり、ピンク色の部分の大きさが検出力であり、(1-β)≒0.2である。

統計学復習メモ15: 幾何分布のパラメーターの最尤推定量

前のメモで間抜けなことをしたついでに、同じ問題についてもう少し考察する。

再現率が低いバグを、(1) Aさんは200回試行して1回発生させ、Bさんは500回試行して1回発生させ、Cさんは1000回試行して1回発生させた場合と、(2) Aさんは200回目で1回発生させ、Bさんは500回目で1回発生させ、Cさんは1000回目で発生させた場合とでは当然数字の意味が異なる。(1)では発生させた後も試行してる可能性があるからである。

発生させたら試行を終了すると考えた場合、発生率をpとすると、k回目で発生する確率は、(k-1)回失敗して1回成功する確率なので、p(1-p)k-1である。簡単にするため、k+1回目で発生する確率をP(k)=p(1-p)kとする。このkの分布には「幾何分布」という名前がついているらしい。幾何級数（等比級数）だからということらしいが、この用語は一般に通じるものなのかどうか疑問だが、今回はこの用語を使う。発生するまでの試行回数の標本X1〜Xnが得られた時、それが幾何分布に従うとすると、発生率pの推定量は何になるだろうか？

前回に続き、最尤推定法を使う。幾何分布の確率をθ、尤度関数（結合確率）をL(θ)とすると、
¥frac{¥partial}{¥partial¥theta}¥log L(¥theta) = ¥frac{¥partial}{¥partial¥theta}¥sum_{i=1}^{n}¥log(¥theta(1-¥theta)^{X_i})
なので、これが0になるθを求める。
¥frac{¥partial}{¥partial¥theta}¥sum_{i=1}^{n}¥log(¥theta(1-¥theta)^{X_i})=¥sum_{i=1}^{n}(¥frac{1}{¥theta}-¥frac{X_i}{1-¥theta})=¥frac{n}{¥theta}-¥frac{1}{1-¥theta}¥sum_{i=1}^{n}X_i=0
で
¥theta=¥frac{n}{n+¥sum_{i=1}^{n}X_i}=¥frac{1}{1+¥bar{X}}
となる。

従って、上記(2)の場合だと、X1=199, X2=499, X3=999なので、p=3/1700と推定される。
(1)の場合は、前のメモの通り、(1/200+1/500+1/1000)/3=1/375が妥当だと思う。

統計学復習メモ14: ポアソン分布のパラメーターの最尤推定量

筆者はある会社のソフトウェア関係の部署に在籍しているが、一般にはソフトウェア開発とは呼ばれない補助的な仕事をしている。調達担当や購買担当という呼び方が最尤推定量であろう。そんな筆者にもお呼びがかかるのが、不再現バグの再現試験である。
バグ発生の報告があったが、開発担当者が解析するために再現させようとしても、再現しない。発生した証拠はある。となると、総出で再現試験である。

そのバグを発生させるべく、100回も200回も同じ操作をしていると、いつまでやらないといけないのだろうと不安になってくる。1万回やっても出ない可能性もあるのである。仮に発生確率が1/1000だったとしても、1000回やったら必ず1回は出る訳ではない。確率が低い事象の発生回数は「小数の法則」ポアソン分布に従う。発生回数（発生確率×試行回数）の平均がλの場合、発生回数がkになる確率は、
P_k=e^{-¥lambda}¥frac{¥lambda^k}{k!}
である。
λ=1のポアソン分布の確率密度
従って、1000回やって平均1回しか起こらないものが1000回やって1回も起こらない確率は1/e≒36.8%もある。1回起こる確率が同じく1/e、従って2回以上起こる確率が(1-2/e)≒26.4%であるが、1回も起こらない確率が36.8%もあるのだ。既に何百回も起こらないでなぜ＋思考になれようか。

自分が毎回その1/eに入るのを不思議に思いながらも、誰かが再現させて解放される。最初に再現させる人は1人しかいないのだから、自分がその1人になる確率は低いので、左脳で考えると当たり前である。複数人でのべ2000回やって1回起こったとすると、結果として発生率は1/2000である。
さて、本題である。この時、本当の発生率も1/2000と推定できるのだろうか？実は本当の発生率が1/5000の時が2000回に1回発生する確率が一番高かったりしないのだろうか？
これはまあ、P(k=1)=λe-λが

こんな形をしてるので、P(k=1)をλで偏微分して=0とすると、(1-λ)e-λ=0でλ=1ということになり、1回発生する確率が最も大きいのはλ=np=1すなわちp=1/2000と力技で求められる。

では、もし、ある人は200回で発生させ、ある人は500回で発生させ、ある人は1000回で発生させたら、その発生率は全て足し合わせて3/1700と推定するものなのだろうか？
問題を整理すると、ポアソン分布のパラメーターλの推定量λ^は何なのだろうか？


今回は、最尤推定法を使うことにする。最尤推定法とは、あるパラメーターθを持つ確率分布を仮定して、結果として標本X1...Xnが得られた時、n個の標本を取ると前記の標本X1...Xnが得られる確率が最も高いθを求める方法である。X1...Xnが得られる確率をP(X1)...P(Xn)とし、X1...Xnが独立（１つ前の標本によって発生率が変わったりしない）とすると、それらが同時に起こる確率はL(θ)=P(X1)×P(X2)×...×P(Xn)である。これが最大になるθ=θ^を求めるのである。さらにL(θ)が単峰（山頂＝極大点が１つしかない）なら、そのθは
¥frac{¥partial}{¥partial¥theta}L(¥theta) = ¥frac{¥partial}{¥partial¥theta}¥prod_{i=1}^{n}P(X_i) = 0
の解である。積の形だと何なので、対数を取って
¥frac{¥partial}{¥partial¥theta}¥log L(¥theta) = ¥frac{¥partial}{¥partial¥theta}¥sum_{i=1}^{n}¥log P(X_i)=0
の解でもあるのであるのである。

ポアソン分布のλについては、
¥frac{¥partial}{¥partial¥lambda}¥log P_k = ¥frac{¥partial}{¥partial¥lambda}¥log e^{-¥lambda¥frac{¥lambda^k}{k!}} = ¥frac{¥partial}{¥partial¥lambda}(-¥lambda+k¥log¥lambda-¥log(k!)) = -1 + ¥frac{k}{¥lambda}
なので、L(λ)が最大になるλを求めると、
¥frac{¥partial}{¥partial¥lambda}¥log L(¥lambda)=¥sum_{i=1}^{n}(-1+¥frac{X_i}{¥lambda})=0
-n+¥frac{1}{¥lambda}¥sum_{i=1}^{n}X_i=0
¥lambda=¥frac{1}{n}¥sum_{i=1}^{n}X_i=¥bar{X}
となるので、標本平均が最尤推定量である。

元の問題に戻って、X1=1/200、X2=1/500、X3=1/1000と置くと、平均λ=(X1+X2+X3)/3=8/3000=1/375となる。

統計学復習メモ13: 有効推定量と不偏性、一致性

我ながら物好きであるという仮説が棄却されなさそうであるが、大学の時に買う羽目になった統計学の教科書を引っ張り出して、推定量というものを勉強している。そこには、推定量の好ましい性質として不偏性と一致性が挙げられ、次に推定量同士を比べる話が出てきて、その後に有効推定量、最尤推定量が出てくる。この教科書のこの流れのために、有効推定量や最尤推定量は不偏性と一致性を満たすんだろうとずっと誤解していた。何せ「有効」「最＊」なのだから。

特に「最尤推定量」の名前の意味は「最ももっともらしい推定量」であり、まるでキャッチフレーズのように「最尤推定」に付け加えられ、そのインパクトによって、これぞ最強の推定量というイメージが刷り込まれてしまう。大学の研究室にも、最尤推定量の意味を理解せず、「推定」という言葉が出てくる度に、馬鹿の一つ覚えのように「最尤推定を使えばいいじゃない」と言っていた教官がいた。実データが取りようが無い（統計量が無い）ものの話をしていても「最尤推定」という単語を持ち出されて、私を含め、学生たちは皆、そのセリフに降参するしか無かった。

余談はさておき、よく読むと、どうやら有効性(efficiency)と不偏性(unbiasedness)、一致性(consistency)とは直接は関係ないようである。例えば標本分散は一致推定量であるが不偏ではない（不偏分散は不偏かつ一致推定量）し、逆に不偏だが一致性が無い場合もある。有効推定量や最尤推定量にも不偏でない（biasedな）ものがあるし、最尤推定量にも一致性が無いものがあったりするらしい。


一致性というのは、数式で書くと、実際のパラメーターをθ、推定量をθ^、Pを確率、nを標本数として
¥exists¥varepsilon ¥lim_{n¥rightarrow¥infty}P(|¥hat{¥theta}-¥theta|>¥varepsilon)=0（または不等号を逆にして=1）
とちょっとアレであるが、要するに標本の数が多ければ推定量が実際の値に近づく性質のことである。意図的に外すかよっぽどチョンボしない限りはそれを満たすのが当たり前な気がするが、推定量が多項式でない場合は気にした方が良いのかも知れない。

不偏推定量であるが一致推定量でない極端な例としては、標本をX1...Xnとして、平均θの推定量θ^=X1やθ^=Xnがある。標本がいくら増えてもその中の１つしか使わなければ精度が上がらないので、当たり前である。標本を全て使っても一致推定量でない例として、こういうのを考えてみた。
¥hat{¥theta}=¥sum_{k=1}^{n}¥frac{X_k}{2^k}+¥frac{X_n}{2^n}
これは、実際の分散をσ2とすると
¥lim_{n¥rightarrow¥infty}Var(¥hat{¥theta})=¥lim_{n¥rightarrow¥infty}¥left(¥sum_{k=1}^{n}¥frac{¥sigma^2}{2^{2k}}+¥frac{¥sigma^2}{2^{2n}}¥right)=¥frac{¥sigma^2}{3}
であり、nを大きくしても分散が0にならないので、θ^はθに収束しない（はず）。（ちなみに、普通にθ^を標本平均(X1+X2+...+Xn)/nとすると、Var(θ^)=σ2/nなので、¥lim_{n¥rightarrow¥infty}Var(¥hat{¥theta})=0となり、θ^はθに限りなく近づく）


有効推定量というのは、クラメル＝ラオの不等式
Var(¥hat{¥theta}) ¥geq ¥frac{1}{n E¥left¥{¥left[¥frac{¥partial}{¥partial¥theta}¥log P(X)¥right]^2¥right¥}} ¥equiv ¥sigma_0^2
の等号が成立する（Var(θ^)が「クラメル＝ラオの下限」σ02になる）ようなθ^、つまり不偏推定量の中で分散が最小になるもののことである。確率の対数をθで偏微分して２乗して平均して逆数にして訳がわからないが、対数が出てくるのは、AとBとCが同時に起こる確率がP(A)×P(B)×P(C)なのでその対数がlogP(A)+logP(B)+logP(C)となるように、標本それぞれが同時に起こる確率が線形和になるためのもので、偏微分は飛ばして、２乗の平均は誤差の２乗和の意味で出てくるのだと思う。

上記の不等式を解く代わりに、
¥hat{¥theta}=¥theta+K¥sum_{i=1}^{n}¥frac{¥partial}{¥partial¥theta}¥log P(X_i)
の右辺がθを含まなくできるK（但しKはXiは含まない）があれば、その時のθ^が有効推定量、という定理が使える。それを使って、X1...Xnが平均μ、分散θの正規分布に従う標本の場合のθの有効推定量を求めてみる。確率密度関数は
P(x)=¥frac{1}{¥sqrt{2¥pi¥theta}}e^{-¥frac{(x-¥mu)^2}{2¥theta}}
なので、
¥frac{¥partial}{¥partial¥theta}¥log P(x)=¥frac{¥partial}{¥partial¥theta}¥left(-¥frac{1}{2}(¥log(2¥pi)+¥log¥theta)-¥frac{(x-¥mu)^2}{2¥theta}¥right)=-¥frac{1}{2¥theta}+¥frac{(x-¥mu)^2}{2¥theta^2}
であり、
¥hat{¥theta}=¥theta+K¥sum_{i=1}^{n}¥left(-¥frac{1}{2¥theta}+¥frac{(x-¥mu)^2}{2¥theta^2}¥right)=¥theta-K¥frac{n}{2¥theta}+K¥frac{1}{2¥theta^2}¥sum_{i=1}^{n}(x-¥mu)^2
なので、K=2θ2/nとするとθが消えることがわかる。従って、
¥hat{¥theta}=¥theta-¥frac{2¥theta^2}{n}¥frac{n}{2¥theta}+¥frac{2¥theta^2}{n}¥frac{1}{2¥theta^2}¥sum_{i=1}^{n}(x-¥mu)^2=¥frac{1}{n}¥sum_{i=1}^{n}(x-¥mu)^2
であり、標本分散がθの有効推定量である。
標本分散は不偏ではないので、正規分布の分散の有効推定量は不偏推定量ではないことがわかる。

統計学復習メモ12: 不偏推定量とは

以前にｔ分布を使うパラメーターの区間推定の方法を覚えたが、この前ちょっと推定の実践のネタを思いついてやってみようとしたら、そのパラメーターが正規分布に従うものではないことに気付き、１手目でつまずいた。これでは勉強した意味が無いと思い、推定というものの基礎を勉強し直すことにした。

区間推定をするにもまず、その区間の基準となる値（推定量）を点推定しないと始まらない。その点推定量が満たすのが好ましい性質として、不偏性と一致性がある。推定量θ^がパラメーターθの不偏推定量であるとは、
E(¥hat{¥theta})=¥theta ...(1)
が成り立つことをいう。

出たE(¥hat{¥theta})=¥theta。学生時代にこのθの上の屋根の意味が理解できなかったことが、私が推定・検定の１手目でつまずいた原因の１つであることは間違いない。そもそもなぜ統計学では推定するパラメーターをθと書くのか。高校の三角関数の授業で習って以来、私にとってはθといえば角度に違いないのであり、パラメーターをθと表すのは、円周率をxと表すとか、y=ax+bのxとyが定数でaとbがグラフの縦軸・横軸であるくらいの混乱の元である。

そもそも、θ^がθの推定値だとしたら、(1)が成り立たないとはどういう状態なのだろうか？それに、θが未知だからθを推定するのに、θ^の平均がθであることをどうやって確認するのか？
...という疑問を持ったのは筆者だけなのだろうか。

これは、θが確率変数のパラメーターであり、その推定量θ^を標本から求める方法を決めた時、仮にθによって決まる母集団からの標本を使った場合、θ^の平均がθに一致するなら、そのθ^は不偏推定量である、という意味らしい。

例題として、X1〜Xnが平均μ、分散σ2の母集団からの標本として、まず標本平均X~=E(X)=1/n ΣXiがμの不偏推定量かどうかを考える。
E(¥bar{X})=E¥left(¥frac{1}{n}¥sum_{i=1}^{n}X_i¥right)=¥frac{1}{n}¥sum_{i=1}^{n}E(X_i)=¥frac{1}{n}n¥mu=¥mu
従ってX~はμの不偏推定量である。

次に、標本分散S^2=¥frac{1}{n}¥sum_{i=1}^{n}(X_i-¥bar{X})^2、不偏分散U^2=¥frac{1}{n}¥sum_{i=1}^{n-1}(X_i-¥bar{X})^2がσ2の不偏推定量かどうかを考える。

ここで、Var(X)=E(X^2)-¥left(E(X)¥right)^2の公式とVar(¥bar{X})=¥frac{1}{n}Var(X)の公式を用いると、
E(X_i^2)=Var(X_i)+¥left(E(X_i)¥right)^2=¥sigma^2+¥mu^2
E(¥bar{X})=Var(¥bar{X})+¥left(E(¥bar{X})¥right)^2=¥frac{¥sigma^2}{n}+¥mu^2
が求まるので、
E(S^2)=¥frac{n-1}{n}¥sigma^2
E(U^2)=¥frac{n}{n-1}E(S^2)=¥sigma^2
となる。従って、S2はσ2の不偏推定量ではなく、U2はσ2の不偏推定量である。

統計学復習メモ11: なぜ特異値分解で主成分が求まるのか

前々項に書いたように、主成分分析は特異値分解を使っても計算できる。なぜだろうか。
それについても答えを得たので、ついでにメモしておく。

ｎ個のサンプルデータを横に並べたものを

(m<n)とし、Pを単位主成分を横に並べたもの

とし、YをXの主成分とすると、
Y=P^TX
である。

特異値分解定理により、Xは何らかの正規直交行列U,Vと、m行m列部分が対角行列でそれ以外が0である行列Σによって、
X=U¥Sigma V^T
と分解される。それを応用すると、
XX^T=(U¥Sigma V^T)(U¥Sigma V^T)^T=U¥Sigma V^TV¥Sigma^TU^T=U¥Sigma¥Sigma^TU^T...(1)
が成り立つ。XXTはm×mの正方行列なので、ΣΣTもm×mの正方行列であり、
...(2)
と置くことができる。（σ1≧σ2≧...≧σmである）

前項に書いたように、Yの共分散行列は
Cov(Y)=P^TCov(X)P
であり、従って
Cov(X)=P\;Cov(Y)P^T
であり、Cov(Y)はCov(X)の固有値を対角成分に持つ行列である。Cov(Y)=Λを

(λ1≧λ2≧...≧λm)と置くと、
...(3)
となる。一方、
Cov(X)=¥frac{1}{n}XX^T
なので、(1)と(2)を使うと、
...(4)
という形にすることができ、(3)と(4)を見比べると、PもUも正規直交行列なので、U=P、
すなわち ¥frac{¥sigma_k^2}{n}=¥lambda_k
の関係があることがわかる。さらに、
Y=P^TX=U^TU¥Sigma V^T=¥Sigma V^T
なので、Xの主成分YはXを特異値分解したものの一部として求められることがわかる。

統計学復習メモ10: なぜ共分散行列の固有ベクトルが単位主成分なのか

前項に書いた通り、主成分分析における主成分の単位ベクトルは、共分散行列の固有ベクトルとして求まる。そのこと自体に昔から興味があったので、主成分分析の復習ついでに考察してみる。

まず、最小２乗法で考えてみる。簡単のために２次元で考える。ｎ個のサンプルデータを

とし、第１主成分の単位ベクトルを

とすると、Xに対応する主成分軸上の第１主成分Yは
Y=P^TX
であり、そのYを元の座標系に戻したものX~は
¥tilde{X}=PY=PP^TX
である。このことは、高校で習った一次変換を思い出してやってみるとわかる。このX~が、Xを第１主成分の軸上に射影したものであり、これとXとの距離が、最小にしたい誤差ということになる。その誤差Eを、Xを直交座標とした場合の距離の２乗とすると、
E=(\tilde{X}-X)^T(\tilde{X}-X)=\sum_{i=1}^n\|\tilde{X}_i-X_i\|^2
であり、p12+p22=1に注意すると、これは
=¥sum_{i=1}^n((p_1(p_1x_{1i}+p_2x_{2i})-x_{1i})^2+(p_2(p_1x_{1i}+p_2x_{2i})-x_{2i})^2)
=¥sum_{i=1}^n(p_2x_{1i}+p_1x_{2i})^2
と計算できる。

このp1とp2の２次関数であるEの最小値ををp12+p22=1の条件下で求めるので、「ラグランジュの未定乗数法」というやつを使う。すると、G=p12+p12-1と置くと、

である。これを整理すると、

となる。これの最初の行列は、よく見ると、Xの共分散行列

の逆行列のスカラー倍であるので、上の方程式は、両辺に左側からCov(X)をかけると

（aはスカラー）となり、a/λをλと置き直すと、

が得られる。従って、固有ベクトルの定義より、PはCov(X)の固有ベクトルであることがわかる。


さて、Xが２次元の場合については辛うじて求まったが、Xが３次元になると、この方法ではEが複雑になり、非常に辛い。筆者はMaximaに式の変形をさせながら３次元についても同じであることを調べようとしたが、E-λGを偏微分した所で変数が目測で200個以上ある状態になり、なんとなく雰囲気はあったが挫折した。この方法ではXをｍ次元に一般化するのはかなり困難か、何か特殊な技術が必要だと思う。

そこで、次に分散最大で考えてみる。こちらは幾分容易である。
ｎ個のサンプルデータを

とし、第１主成分の単位ベクトルを

として、第１主成分

の分散
V=¥frac{1}{n}YY^T=¥frac{1}{n}(P^TX)(P^TX)^T=¥frac{1}{n}P^TXX^TP
を最大にするPを求めることを考える。Xのi行目とj行目の共分散をcov(Xi,Xj)と書くと、
V=¥sum_{i=1}^{m}¥sum_{j=1}^{m}p_ip_jcov(X_i,X_j)
とも書ける。このp1...pmの２次関数であるVの最大値をp12+p22+...+pm2=1の条件下で求めるので、やはりラグランジュの未定乗数法を使う。すると、G=p12+p22...+pm2-1と置くと、

である。これを整理すると、

すなわち

が得られ、PがCov(X)の固有ベクトルであることがわかる。

第２主成分以降についても、それより前の主成分と直交するという条件を加えて、同様の方法で解くと、やはりCov(X)の固有ベクトルであることがわかり、第１主成分から順番に分散最大の固有値が取られていくので、第k主成分はCov(X)のk番目に大きい固有値に対応する固有ベクトルだと求まる。

統計学復習メモ9: 主成分分析の計算方法

主成分分析の話になると大抵、固有ベクトルか特異値分解のどちらかが出てくる。主成分の単位ベクトルが固有ベクトルとして直ちに求められ、主成分の大きさが特異値分解を使って直ちに求められるからだ。

それにしても、行列には固有値や固有ベクトルがあるというのも、統計には分散や共分散があるというのも、理系の大卒なら誰でも知ってることだろうが、「共分散行列の固有ベクトルが主成分の単位ベクトルになる」と言われると、全く関係ない３つのものが突然シンプルに組み合わせられた感じで、衝撃的である。まるで
e^{-i¥pi}=-1
のような神秘的なものを感じるのは筆者だけだろうか。

そのこと自体は、至る所に証明というか説明があり、その方法も幾通りもある。基本的には、あるデータのど真ん中を通る直線として、最小２乗法のように２乗誤差を最小にする直線を求める代わりに、その直線方向の分散が最大になる（その直線上に射影した場合の分散が最大になる）直線を求める、という考え方である。そう言われればそういう考え方もありかな？とも思ったりするし、やっぱり最小２乗法で求める方がいいんじゃないの？とも思ったりするが、結局、主成分は分散最大で求めても最小２乗法で求めても同じ結果になるらしい。


とりあえず、固有ベクトル計算でやってみる。
m次元のデータn個を

（iは標本番号、pはパラメーター番号の意味）とし、共分散行列の固有値を

とし、対応する固有ベクトルを

とすると、W1〜Wmが第１主成分〜第ｍ主成分、λ1〜λmがそれぞれの寄与率となる。
X,Wは行列である。プログラミング用語の「待ち行列」とかのqueueではなく、線形代数のmatrixである。フフン、行列くらい知っとるわい、という顔をしつつ、厄介だなあ、と思わざるを得ないが、主成分分析するのに行列は避けて通れない。
前回使用したデータについて、Maximaに計算させてみる。

Maximaへの入力

load(descriptive);
load(lapack);
fpprintprec:6;
X: read_matrix("data2");
[L,W,dummy]: dgeev(cov(X),true,false);
L;
W;
Maximaの出力（上のをbatch()で読み込んだ）
(%i7) L
(%o7) [.330129, .0531125, .00529463]
(%i8) W
[ .714134 .697955 - .0535787 ]
[ ]
(%o8) [ 0.4882 - .551443 - 0.67644 ]
[ ]
[ .501671 - .456912 .734546 ]
固有値（寄与率）は絶対値が一番大きいのが１つ目、二番目に大きいのが２つ目なので、第一主成分の単位ベクトルは(0.714, 0.488, 0.502)、第二主成分の単位ベクトルは(0.698, -0.551, -0.457)である。３つ目の固有値は小さいので、３つ目の固有ベクトルは無視することにする。
今回はたまたま大きい順に並んでいるが、固有値は絶対値の大きさ順に並ぶとは限らないので注意が必要である。
もしdgeev()が使えない場合は、eigens_by_jacobi()を使う。
(%i9) [L,W]:eigens_by_jacobi(cov(X)); L; W;
(%o10) [.330129, .0531125, .00529463]
[ .714134 - .697955 - .0535787 ]
[ ]
(%o11) [ 0.4882 .551443 - 0.67644 ]
[ ]
[ .501671 .456912 .734546 ]
２つ目の固有ベクトルの符号が逆だが、次に求める主成分の大きさの符号が逆になるだけなので、どちらでも良い。

次に、第一主成分、第二主成分を軸に取った時にXの各標本がどういう値になるか（主成分の大きさ）を求める。主成分の大きさYは
Y=XW
なので、

W: submatrix(W,3); /* 弟三主成分を削除 */
Y: X.W;
で求まる。これによって得られるYをグラフにすると、次のようになる。

前回のグラフと見比べると、視点が反対（前回のグラフ作成の時はengens_by_jacobi()を使ったため）だが、例えば楕円外の点の分布は大体似ていることがわかる。
これによって、今回使った３次元のデータが、情報の欠落を最小限にして、２次元で分析できる。


特異値分解による方法でもやってみる。
今度はXの中の標本は横に並べる。縦方向の１列が１つの標本である。

このXが特異値分解によって
X=U¥Sigma V^T
（Tは転置）と３つの行列（U=左特異行列、Σ=特異値行列、V=右特異行列）の積になる時、Uが主成分の単位ベクトルを縦に並べたもの、Σの対角成分（特異値）が［ｎ×寄与率］の平方根となり、標本の主成分の大きさYは
Y=U^TX=¥Sigma V^T
となる。
Maximaで特異値分解するには、LAPACKのdgesvd()を使うと一発である。（LAPACKが使えないと、かなり厳しそうである）

Maximaへの入力

load(lapack);
fpprintprec:6;
X: transpose(read_matrix("data2"));
[s,U,V]: dgesvd(X,true,true)$
U;
s;
Maximaの出力
(%i7) U
[ - .713141 .698967 - 0.053619 ]
[ ]
(%o7) [ - .488765 - 0.55059 - .676727 ]
[ ]
[ - .502532 - .456395 .734279 ]
(%i8) s
(%o8) [5.77344, 2.30564, .727852]
（Vは100x100の行列なので省略）
行列Uに、先程と大体同じ主成分の単位ベクトルがあることと、３つ目の特異値が小さく、寄与率が低いことがわかる。

Maximaへの入力（続き）

S: zeromatrix(2,100); /* 行列Σ(3x100)の上２行分 */
S[1,1]: s[1]; /* 対角成分１ */
S[2,2]: s[2]; /* 対角成分２ */
Y: S.V;
これによって得られるYをグラフにすると、次のようになる。

さらに反対向きになったが、そういうものである。


さて、固有ベクトル計算による方法と特異値分解による方法をどう使い分けるかだが、一般に、特異値分解を使う方法の方が計算誤差が少ない、と言われる。固有ベクトル計算による方法だと、共分散行列にする時に一気に情報量が落ちるので、感覚的には納得できるが、使用するモジュールやライブラリの精度による所もあるので、一概には言えないと思う（例えばすごい桁数の共分散行列から固有ベクトルを求めれば精度は上がるはず）。

また、ライブラリの特性の問題かも知れないが、eigens_by_jacobi()やLAPACKのdgeev()では固有ベクトルが固有値の絶対値の順に並ぶとは限らないので、運が悪いと固有値を見て固有ベクトルを並び替える必要が生じる。それに対し、LAPACKのdgesvd()は特異値が必ず大きい方から並ぶので、LAPACKを使うなら特異値分解による方法の方が楽かも知れない。

統計学復習メモ8: 回帰直線と主成分分析

高校生の時に、必要も無いのに関数電卓を買った覚えがある。物理の先生に、理系ならずっと使うから早目に買っといて損は無いとか言われたからだと思うが、結局大学で物理実験をやるまでは遊びでしか使わなかった。それでも、元々電卓が好きだった私は、毎日毎日、関数電卓のプラスチックのボタンをカチャカチャと押して遊んでいた。
その電卓には、「一次回帰直線」の係数を求める機能があった。学校では習わなかったが、説明書にy=a+bxのaとbを求める機能と書いてあったので、意味は分かった。来る日も来る日も自室に逃げ込んで電卓と戯れてる内に、「一次回帰直線」がとても身近なものになった。

そんなに一次回帰直線に慣れ親しんだから、主成分分析には敏感に反応した。というのも、主成分分析は繰り返し回帰直線を求めることに似ていると思うからである。

一次回帰直線とは、パラメトリックな２変量の関係を直線で近似するものなので、２次元のデータを１次元で近似することに相当する。２変量をx,yとすると、標本(x1,y1)〜(xn,yn)の関係を直線y=a+bx（a,bは定数）で近似するので、その直線との誤差を無視して(x,y)=(0,a)+(1,b)tと置くと、２変量の標本(xi,yi)が１変量tiで表せるのである。
主成分分析の第一主成分だけを使うことを考えると、多変量の関係を１次元で近似することになるから、k次元の標本(x1i,x2i,...,xki)を(C1,C2,...,Ck)+(p1,p2,...,pk)tiに近似すると、k個の変量を持つ標本(x1i,x2i,...,xki)が１変量tiだけで表せる。(C1,C2,...,Ck)は定数ベクトルで、(p1,p2,...,pk)が第一主成分の単位ベクトルである。第二主成分の単位ベクトルを(q1,q2,...,qk)としてこれも使うとすると、(x1i,x2i,...,xki)≒(C1,C2,...,Ck)+(p1,p2,...,pk)si+(q1,q2,...,qk)tiという感じに、k変量の標本を２変量で近似することになる。

一次回帰直線の最小２乗法による求め方をおさらいする。標本(xi,yi)と直線y=a+bxとの誤差はyi-a-bxiであるが、この誤差の２乗和を最小にするようにa,bを決めたい。誤差の２乗和をEとすると、E=Σ(yi-a-bxi)2で、このEはaの２次関数であり、bの２次関数でもあるから、Eをaで偏微分したものが0になるaと、Eをbで偏微分したものが0になるbを求めれば良い。
\frac{\partial E}{\partial a}=0, \frac{\partial E}{\partial b}=0
これを解くと、
b\frac{\sum xy - \frac{1}{n}\sum x\sum y}{\sum x^2-\frac{1}{n}(\sum x)^2},a=\bar{y} - b\bar{x}
が得られる。bの分子はxとyの共分散、分母はxの分散なので、xの分散をVar(X)、xとyの共分散をCov(X,Y)と表すと、b=Cov(X,Y)/Var(X)とも書ける。実はシンプルである。

念のため、試しにやってみて確認する。例によってMaximaを使う。サンプルデータをdata1とし、

data1: read_nested_list("data1");
datam: read_matrix("data1");
load(descriptive);
b: cov(datam)[1][2] / var(datam)[1];
a: mean(datam)[2] - b * mean(datam)[1];
とすると、b=-.8338664607943058、a=5.848350705139143と求まる。そこで
plot2d([[discrete, data1], a+b*x],[x,0,2],[style,[points],[lines]]);
とすると、次のグラフが作られる。

主成分分析についても、詳細を省いて、とりあえずサンプルデータ(data2)を使って図示してみる。
下のように、標本が原点を中心に楕円盤状に散らばっているとする。（どれかの画像をクリックするとgifアニメが開きます）



少しずつ回転させてみたが、わかりづらいので、大体の分布の形を枠線で表示する。


第一主成分の単位ベクトルは最も誤差の和を小さくする直線の方向、すなわち楕円の長径方向の長さ１のベクトルになる。次の図の青線が第一主成分（をスカラー倍したもの）である。


標本から第一主成分を抜く（第一主成分を0にする）と、下の図のようになる。（回転方向は変更）


つまり、第一主成分と垂直な平面への射影である。元が楕円盤状だったのでこれも楕円状であり、この楕円を最も近似した直線、すなわち長径方向が第二主成分の向きになる。（下図緑線）


３次元なので、第二主成分も抜くと、直線状になる。その直線は第一主成分、第二主成分と直交し、その直線上に第三主成分がある。

空間が何次元でも同様である。最も誤差が小さくなる近似直線の方向の大きさ１のベクトルが第一主成分の単位ベクトルであり、その方向のベクトル成分を抜いた（超平面に射影した）ものの近似直線の方向の大きさ１のベクトルが第二主成分の単位ベクトルであり、その方向のベクトル成分も抜いて...を繰り返すと、全ての主成分が出てくる。

そんな反復計算をせずに主成分分析する方法は、次回まとめる。

統計学復習メモ7: t検定で独立性の検定

前のエントリーではカイ２乗検定を用いてノンパラメトリック（定性的）な独立性の検定を試したが、次はパラメトリック（定量的）な場合をやってみる。独立かどうかを確かめたい２つの属性が定量的な値を持つパラメーター、Ｃ言語で言うとenum型やchar型でなくint型やfloat型の変数である場合に、その２変数間に相関があるかどうかを検定するものである。

なお、統計学では独立、相関という言葉はノンパラメトリックかパラメトリックかで使い分けるのが通常なのかも知れないが、ここでは区別せず、独立である／ないと相関がない／あるは同じ意味で使う。

まず、実際に例題でやってみる。

ある事業分野において、過去に優秀なプログラマーが書いた、読んで理解し易いしその後の仕様変更にも強いと評判の良い、模範的なＣ言語のソースコード一式があったとする。ソースファイル中のコード量（文の数や式の数、a=0;やb=foo(a);は１つ、while(式)はwhileと式とで２つ、for(式１;式２;式３)は４つ等）とコメントの量（/*～*/内の文字数）の関係が以下であった時、コード量とコメント量との間に相関がある無いかを検定したいとする。

ファイル名コード量コメント量
gentoo.h88234
emperor.h98413
macaroni.h485687
rockhopper.h285438
fairy.h184194
magellan.h487174
humboldt.h818344
king.h69139
adelie.h306244
cape.h219222
ファイル名コード量コメント量
gentoo.c5191876
emperor.c5152995
macaroni.c78517
rockhopper.c3131839
fairy.c3922142
magellan.c3691783
humboldt.c3452762
king.c73489
adelie.c7292202
cape.c4452309

統計学で相関と言えば、まず思い付くのは相関係数である。とりあえず、ヘッダファイルとソースファイルに分けて、コード量とコメント量の相関係数を求めてみる。
よく使われる相関係数は、共分散÷それぞれの分散の幾何平均、というもので、式で書くと、２つのパラメーターを持つ標本が(x1,y1)～(xn,yn)の時、
r_{xy}=¥frac {¥sum_{i=1}^{n}(x_i-¥bar{x})(y_i-¥bar{y})} { ¥sqrt{ ¥sum_{i=1}^{n}(x_i-¥bar{x})^2 ¥sum_{i=1}^{n}(y_i-¥bar{y})^2 }}
と表されるものである。（ピアソンの積率相関係数という名前が付いているらしい。）
例によってMaximaに計算させると、ヘッダファイルについてはr≒0.327、ソースファイルについてはr≒0.746であった。

さて、この値をどう考えれば良いのだろう。２つのパラメーターに十分に直線関係があるかどうかの基準として使われるものの１つに、相関係数が0.99以上（強い正の相関）または-0.99以下（強い負の相関）か、というのがある。検量線などの回帰直線の適切さの基準として使われることがあるようだが、そのように基準を決める自由がある場合ならそれでいいだろうが、一般的には（一般論としては）0.99という数字に根拠が無い。では相関係数がどういう値なら何なのだろうか。

そこで、検定統計量
T=¥frac{r_{xy}¥sqrt{n-2}}{¥sqrt{1-r_{xy}^2}}
が自由度n-2のｔ分布に従うという定理を使うのが、統計学の尊い教えである。

仮説を「相関がない」とし、有意水準を0.05（信頼度を95%）とすると、自由度n-2のｔ分布で右側（ｔが大きい方）の棄却域はt≧2.3くらいである。それに対し、ヘッダファイルについてはT≒0.978、ソースファイルについてはT≒3.17なので、ヘッダファイルについては仮説が棄却されず、コード量とコメント量の間に相関があるとは言えない、ソースファイルについては仮説が棄却され、コード量とコメント量の間に（正の）相関がある、ということになる。


結局、この方法では、相関係数はいくらぐらいだと相関があると言えるのだろうか。
まず、相関係数rと検定統計量Tの関係を見ておく。n=5,10,15,20の時のTのグラフは、次のようになる。
r-T
単調増加、点対称である。
Tがｔ分布に従うと、rの分布は次のグラフのようになる。
P(r)
当たり前だが、rが-1（負の相関）や+1（正の相関）に近いのは稀だということになっていることがわかる。
Tがｔ分布の両側5%の境目になるrの絶対値は、以下の値である。

nr
50.878
100.632
150.514
200.444
相関係数の絶対値がこれらより大きければ、相関があるということになる。
思ったより緩い基準である。

統計学復習メモ6: χ2検定で独立性の検定

以前のメモではカイ２乗検定の例としてばらつきの検定（適合度の検定）をやってみたが、もう１つのカイ２乗検定の使われ方として、独立性の検定がある。これもやってみる。

プログラマー35歳定年説というのがある。あるプロジェクトにおいて、
35歳未満で、受けたバグ指摘件数が20件未満のプログラマーが45人、
35歳未満で、受けたバグ指摘件数が20件以上のプログラマーが15人、
35歳以上で、受けたバグ指摘件数が20件未満のプログラマーが25人、
35歳以上で、受けたバグ指摘件数が20件以上のプログラマーが15人、
だったとする。表にすると、

O(a,b)B<20B≧20
age<354515
age≧352515
（Bはバグ票数）である。暗雲が立ちこめてきたが、信頼度を95%として、35歳以上であることとバグ票数20件以上であることに相関はあるだろうか。

まず、age≧35とB≧20が独立の場合の理論値を計算する。age<35が60人、age≧35が40人、B<20が70人、B≧20が30人なので、これだけから考えると、B<20の70人は60:40でage<35とage≧35に分かれるので、それぞれ42人、28人である。同様に、B≧20の30人は18人、12人に分かれる。

E(a,b)B<20B≧20
age<354218
age≧352812
帰無仮説を、相関がない（独立である）と取ると、観測値が理論値から大きくばらついていれば棄却されるので、ピアソンのカイ２乗検定が使える。χ2分布に従う、ピアソンの検定統計量
T=￥sum_{i=1}^n￥frac{(O_i-E_i)^2}{E_i}
を計算すると、T=(45-42)2/42 + (25-28)2/28 + (15-18)2/18 + (15-12)2/12≒1.79となる。このTは、4つの理論値を固定すると、4つの観測値のどれか1つが決まると定まるので、自由度は1である。自由度が1のχ2分布の左側が95%になるχ2の値は3.84なので、今回のTはばらついていない方（χ2が小さい＝理論値の通りに近い方）の95%に入っており、仮説は棄却されず、相関があるとは言えないことになる。

これは、B≧20かどうかで分けたのが明暗を分けた可能性があり、もしB≧25で分けると

O(a,b)B<25B≧25
age<35555
age≧353010
だったりした場合はT=5.23でOuchとなるので、あまりいい例ではないかも知れないが、そこまでは敢えて考えないことにする。


同じ考え方で、2つの属性で2x2に分類したデータから2つの属性の独立性を検定する方法を定式化してみる。1か2かを取る2つの属性をa,bとし、母数がnの、それぞれの属性の組み合わせに属する標本数の観測値を

O(a,b)b=1b=2
a=1x11x12
a=2x21x22
(x11+x12+x21+x22=n)とする。a,bを独立とした場合の期待値は、b=1について(a=1):(a=2)=(x11+x12):(x21+x22)、a=1について(b=1):(b=2)=(x11+x21):(x12+x22)...となるので、
nE(a,b)b=1b=2
a=1(x11+x12)(x11+x21)(x11+x12)(x12+x22)
a=2(x11+x21)(x21+x22)(x12+x22)(x21+x22)
の定数倍である。これらを全て足すと(x11+x12+x21+x22)2=n2なので、これらが期待値のn倍であり、1/nすれば期待値になることがわかる。
E(a,b)b=1b=2
a=1(x11+x12)(x11+x21)/n(x11+x12)(x12+x22)/n
a=2(x11+x21)(x21+x22)/n(x12+x22)(x21+x22)/n
これを、
O(a,b)b=1b=2sum
a=1x11x12a1
a=2x21x22a2
sumb1b2n
という風に行和、列和を取って整理すると、
E(a,b)b=1b=2
a=1a1b1/na1b2/n
a=2a2b1/na2b2/n
となる。検定統計量Tは、各セルの(O-E)2/Eを合計したものなので、
T=¥sum_{i=1}^{2}¥sum_{j=1}^{2}¥frac{(x_{ij}-¥frac{a_ib_j}{n})^2}{¥frac{a_ib_j}{n}}
と書ける。これを自由度1のχ2分布のχ2の値として見れば良いわけである。

aがp個の値、bがq個の値を取る時も同じことが言えるので、
T=¥sum_{i=1}^{p}¥sum_{j=1}^{q}¥frac{(x_{ij}-¥frac{a_ib_j}{n})^2}{¥frac{a_ib_j}{n}}
（aiはa=iの標本の和、bjはb=jの標本の和）
となる。自由度は、aについては(p-1)、bについては(q-1)なので、全体としては(p-1)(q-1)である。

以上を踏まえて、もう１つ例題をやってみる。原理はわかったので、手計算チックな計算は卒業して、次はMaximaに組み込みのカイ２乗検定関数を使おう、と思って関数を探すと、手元のMaximaには入ってなかった。ExcelにはCHI2TEST()という名前でデフォルトで入ってるのだが、Maximaにはデフォルトでは入っていないらしい。上記の計算だけなら大して複雑じゃないので、自分で関数を作ってみる。 

chi2test(x_,row_,col_):=block([i,j,n,a,b,chi2],
    a:map(lambda([i],sum(x_[i][j],j,1,col_)),create_list(i,i,1,row_)),
    b:map(lambda([j],sum(x_[i][j],i,1,row_)),create_list(j,j,1,col_)),
    n:sum(a[i],i,1,row_),
    chi2:sum(sum((x_[i][j]-a[i]*b[j]/n)^2/(a[i]*b[j]/n),j,1,col_),i,1,row_),
    print("chi^2=",chi2),
    print("95% confidence at",quantile_chi2(0.95,(row_-1)*(col_-1))),
    print("90% confidence at",quantile_chi2(0.90,(row_-1)*(col_-1))),
    chi2
)$
試しに
chi2test([[45,25],[15,15]],2,2),numer;
（numerは結果が分数になるのを避けるため）とやってみると、 
chi^2= 1.785714285714286
95% confidence at 3.841458820694166
90% confidence at 2.705543454095465
と出力される。信頼度90%でも、２つの属性が独立である仮説は棄却されないことになる。

次の例は、バグの原因を設計ミス要因かコーディングミス要因かに分けて、さらにそれをプログラマーの年齢層で分けて数えてみたという架空のデータである。

O(a,b)設計ミスコーディングミス
30歳以上35歳未満176195
35歳以上40歳未満85139
40歳以上45歳未満90106
年齢層とバグの原因との間に相関はあるだろうか。

chi2test([[176,195],[85,139],[90,106]],3,2),numer;
これを実行すると、次の出力が得られる。 
chi^2= 5.35085981290286
95% confidence at 5.991464547107982
90% confidence at 4.605170185988094
従って、信頼度を95%とすると、年齢層とバグの原因が独立であるという仮説は棄却されない。
信頼度を95%とすると、である。

統計学復習メモ5: χ2分布とχ2検定

以前のエントリーにて実際にやってみたが、カイ２乗分布とカイ２乗検定に軽く触れてみると、χ2の定義が違うので混乱した。カイ２乗分布で定義されるχ2は
¥chi^2=¥sum_{i=1}^n ¥left(¥frac{x_i-¥mu_i}{¥sigma_i}¥right)^2
（xiは標本値、μiは平均、σiは標準偏差）であるのに対して、カイ２乗検定で使われるχ2は
¥chi^2=¥sum_{i=1}^n¥frac{(O_i-E_i)^2}{E_i}
（Oは観測値(observed)、Eは期待値(expected)または理論値）である。何でそれらをχじゃなくわざわざχ2と置くんだ、ということも少し引っ掛かるが、それは分散のσ2の例もあるので、値の意味的にきっと何か２乗的、２乗系、２乗fulなんだろうと思って通り過ぎるとして、２つのχ2の式は関連がありそうに見えて何と何が対応するのかよくわからない。Σ内の分母は上のが分散なのに対して、下のは期待値である。分散で割るのと期待値で割るのは違うであろう。分母だけ見るからそう思うんだろうかと思ってΣ内全体を見ても、上のはお馴染みの（xが正規分布N(μ,σ2)に従う時の）標準正規分布N(0,1)に従う値の２乗であるし、下のは平均からの差の２乗を平均で割ってて、平均の大きさに比例しそうな、２乗的には見えない、奇妙な値である。まるで１つ期待値で割り忘れたかのような形をしている。

結論から書くと、上の２つの式は、何か対応がありそうな形はしているが、真正面からパーツ毎に対応させてようと考えるのはやめた方が良さそうだ。
落ち着いてWikipediaを眺めていると、「ピアソンのカイ２乗検定」という言葉が目に入った。カイ２乗検定というのも色々ある中で、以前のエントリーに書いた検定の例は、最もポピュラーな、ピアソンのカイ２乗検定というやつらしい。数学的に説明するのは大変難しいが、検定統計量
T=¥sum_{i=1}^n¥frac{(O_i-E_i)^2}{E_i}
がカイ２乗分布に従うことがピアソンによって示されているので、カイ２乗分布を使って、この検定統計量がある値以上になる確率を知ることができるということだ。
Tは、各事象の実際の発生回数と期待値とのずれ（ばらつき）が大きいほど大きくなるので、観測されたTが十分低い確率でしか起こらないかどうか、すなわちそのばらつきが十分低い確率でしか起こらないかどうかを、このピアソンの方法で知ることができる。
逆に、実際の発生回数を真実、期待値を仮説とすると、Tが十分大きければ、もしその仮説が正しいとするとかなり稀なことが実際には起こったということになり、むしろその仮説が間違っていると考える方が合理的であるので、ある仮説を確率的に誤りと見なす（仮説を棄却する）のに使えるということである。

N(0,1)に従う標本の２乗和と、カイ２乗検定の検定統計量の両方がχ2と定義されるから、学ぶ人が混乱するんだ、と思うのは筆者だけだろうか。きっと何か合理的な理由があるんだろうが、それならその理由を先に知っておかないと混乱の元になると思うが、筆者には未だに見つけられていない。

統計学復習メモ4: （思考実験）学力偏差値と合格率

統計解析の入門書を初めから読み返していると、偏差値の説明がちょろっとだけ書かれていた。ああ、あの高校受験や大学受験の模試で出てきたやつだな、と反射的に思った。
筆者が受験生だったのは遠い昔で、入試というものはほとんど意識することは無くなったのだが、それでも私にとって「偏差値」と言えば入試である。というか、これまでの人生で入試関連以外で50が中心の「偏差値」が使われているのを見た記憶が無い。IQテストなんかは点数の分布が正規分布に従っていることを前提にしてるらしいので、これも広い意味では「偏差値」だろうが、一般的な定義では平均値の偏差値が50で、IQの平均点は100なので、狭い意味では偏差値とは言わないと思う。
従って、偏差値という文字を見ると、受験生時代に見た模試の結果の偏差値しか思い浮かばない。

そういえば、偏差値が60とか70とかってどういうことなんだろう？
我ながら貴様本当に理系だったのか？という感じだが、少なくとも私は授業で習わなかったし、試験には出てこなかった。
定義的には何のことはない、偏差値50が平均で、60が+σ、70が+2σだ。正規分布表を見ると、60で大体上位16%、70で上位2.3%らしい。やっぱり、だから何？という感じである。


私は高校時代、偏差値の意味を知らなかったので、模試の偏差値の値を意識することがなく、自分の学力偏差値がどれくらいだったのか覚えていない。しかし、共通一次から名前を変えて間もないセンター試験の点数は覚えていたので、ウェブ上で過去の平均点と最近の標準偏差を調べて（過去の標準偏差は見つからなかった）概算してみると、58くらいだった。自分が受けた大学のは覚えてない（意識すると凹むから）が、東大の偏差値が70だったのをかすかに覚えている。
仮に当時の私の学力の偏差値が58だったとして、仮にもしセンター試験で偏差値70の点を取ったら東大に合格したとすると、私の東大合格率はどれくらいだっただろうか？

そもそもどうやったらそういうのが計算できるのかがさっばりわからなかったので、統計解析の本を眺めながら考えてみた。

上の仮定と命題にもかなり無理があるが、さらに無理な仮定を重ねる。
800点満点で、受験者の点数が正規分布し、平均点が480、標準偏差が120とする。私が必ず偏差値58の点を取ると必ず落ちるので、私の点数も正規分布するとする。標準偏差が不明だが、感覚的な理由で40だったとする。支離滅裂になってきたが、続ける。つまり母集団の点数がN(480,120^2)に従い、私の点数がN(偏差値58の点,40^2)に従い、私の点数が偏差値70の点を超えると東大合格だった、という仮定である。偏差値58とか70とかはN(480,120^2)における話なので、偏差値xの点は480+120*(x-50)/10である。
Maximaに

cdf_normal(480+120*(70-50)/10, 480+120*(58-50)/10, 40),numer;
（私の点数の分布における、偏差値70の点までの累積密度）と入力すると、0.99984という答えが出たので、私の東大合格率は0.00016、すなわち0.016%、6250中1回の確率ということになる。毎日東大を受験しても17年に1回しか合格しない確率だ。

もし私が１浪して偏差値を5上げたとすると、

cdf_normal(480+120*(70-50)/10, 480+120*(63-50)/10, 40),numer;
は0.982と出てくるので、合格率は1.8%である。

...命題と仮定に無理があり過ぎて、話にならないのは、言うまでもない。
とりあえず、こういう計算にはいろんなデータが必要になることがわかった。