メイン

数学 アーカイブ

2008年05月31日

円盤の衝突計算

弾性衝突円盤アプレットの衝突計算について、悲しいことに式を求めるのにかなり時間を費やしたので、忘れないように記録する。これが無いと、自分でもソースコードを理解できなくなってしまう。

まず直線上の衝突について考える。物体1,2の質量をm1,m2、衝突前の速度をv1,v2、衝突後の速度をv1',v2'とすると、完全弾性衝突の場合、運動量保存の法則より
m_1v_1+m_2v_2=m_1v'_1+m_2v'_2
-\frac{v'_1-v'_2}{v_1-v_2}=1
が成り立つ。速度変化に注目してこれを整理すると、
v'_1-v_1=\frac{-2m_2(v_1-v_2)}{m_1+m_2}
v'_2-v_2=\frac{-2m_1(v_2-v_1)}{m_1+m_2}
が得られる。これは、
v_{tmp}=\frac{m_1v_1+m_2v_2}{m_1+m_2}---(1)
と置くと、
v'_1-v_1=2(v_{tmp}-v_1)---(2)
v'_2-v_2=2(v_{tmp}-v_2)---(3)
と整理できる。

次に平面上の衝突について考える。
ball_hit.png
図のように円盤と円盤が衝突する場合の速度変化については、円盤の中心同士を結ぶ線上の方向の
速度が衝突によって変化し、衝突方向でない速度成分は変化しないので、衝突方向の速度V1r、V2rについて、上記の式(1)~(3)により速度の変化を求めれば良い。つまり、 Pを円盤の中心間を結ぶベクトル、V1,V2を円盤1,2の速度ベクトル、v1r,v2rを円盤1,2のP方向の速度、V1,v1rの衝突後の速度をそれぞれV1',v1r'とすると、円盤の速度変化は
\vec{V'_1}-\vec{V_1} = \frac{\vec{P}}{|\vec{P}|}(v'_{1r}-v_{1r})
である。

v1rの計算については、θをPとV1の間の角度とすると、
v_{1r} = |\vec{V_1}|\cos{\theta}
であり、PとV1との内積は
\vec{V_1}\cdot\vec{P} = |\vec{V_1}||\vec{P}|\cos{\theta}
だから、
v_{1r} = \frac{\vec{V_1}\cdot\vec{P}}{|\vec{P}|}
で求められる。同様に
v_{2r} = \frac{\vec{V_2}\cdot\vec{P}}{|\vec{P}|}
となる。

直線の場合の式(1)と同様に
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2} \left (=\frac{1}{|\vec{P}|}\frac{m_1<br />
(\vec{V_1}\cdot\vec{P})+m_2(\vec{V_2}\cdot\vec{P})}{m_1+m_2}\right )
と置くと、(2)より
¥vec{V'_1}-¥vec{V_1} = ¥frac{¥vec{P}}{|¥vec{P}|}(v'_{1r}-v_{1r}) = 2¥frac{¥vec{P}}{|¥vec{P}|}(v_{tmp}-v_{1r})
となる。vtmpとv1rの両方に1/|P|が入っているので、計算上の便宜のため、
v_{1r} = \vec{V_1}\cdot\vec{P}
v_{2r} = \vec{V_2}\cdot\vec{P}
と置き直すと、
¥vec{V'_1}-¥vec{V_1}=2¥frac{¥vec{P}}{|¥vec{P}|^2}(v_{tmp}-v_{1r})
¥vec{V'_2}-¥vec{V_2}=2¥frac{¥vec{P}}{|¥vec{P}|^2}(v_{tmp}-v_{2r})
が得られる。2次元のベクトルをX軸成分とY軸成分に分けて
\vec{P}=(P_x,P_y), \vec{V_n}=(V_{nx},V_{ny})
とスカラー表記にすると、
|\vec{P}|^2 = P_x^2+P_y^2, \vec{V_n}\cdot\vec{P}=V_{nx}P_x+V_{ny}P_y
なので、円盤同士の衝突による速度変化は、
\vec{V'_1}-\vec{V_1} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{1r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{1r}) \right )
\vec{V'_2}-\vec{V_2} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{2r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{2r}) \right )
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2}
v_{1r} = V_{1x}P_x+V_{1y}P_y
v_{2r} = V_{2x}P_x+V_{2y}P_y
で求められることがわかる。

2008年07月26日

Maximaでドラゴンカーブ

MathematicaとかMapleとかの数式処理ソフトは、何であんなに高いのだろう。2008/7/26現在、個人向けライセンス価格はMathematicaは40万円超え、Mapleも30万円超えだ。教育機関にも政府機関にも属していない、日曜数学ファンの私にはとても手が出せない。
会社でも研究職とかでなく、ソフトウェア開発の付帯業務をやってるので、数学は一切使っておらず、買ってもらえる見込みは無い。バグ収束曲線が云々で必要とか言っても、そんな統計処理ならExcelで十分と言われてしまうであろう。
従って、MathematicaやMapleに触れる機会が無い。

そんな折、Maximaというフリーの強力な数式処理ソフトがあることを知ったので、今回ヤボ用があって、インストールして使ってみた。
とりあえず、テーラー展開させてみる。

 KURO-BOX> maxima

Maxima 5.10.0 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) taylor(sin(x),x,0,10);
                          3    5      7       9
                         x    x      x       x
(%o1)/T/             x - -- + --- - ---- + ------ + . . .
                         6    120   5040   362880
(%i2) taylor(sqrt(r^2-y^2),y,0,10);
                    2      4      6         8       10
                   y      y      y       5 y     7 y
(%o2)/T/       r - --- - ---- - ----- - ------ - ------ + . . .
                   2 r      3       5        7        9
                         8 r    16 r    128 r    256 r
(%i3)
まあ言いたいことは十分わかるが、表示がなんか残念だ。
そこで、満足のいく表示を得るためにTeXmacsを使ってみる。以下の画面はMaxima 5.11.0とTeXmacs 1.0.6.11を使って実際に表示されたものである。%i??という青い部分が入力した部分だ。(表示環境はFreeBSD 6.3+Xorg 1.4.0)

積分もしてみる。

e^(-x^2)を-∞から+∞まで積分すると√πになるのは、手計算でやると面積分にして極座標変換してと技を尽くして偶然解ける感じの有名な例題だが、Maximaは区分求積による近似値計算ではなく、あっさり解析的に計算したようだ。

せっかくTeXmacsを使ってるので、見た目が派手な連分数を表示してみる。

cf(continued fraction)の出力は連分数のリスト表記になっていて、それをcfdisrepすると連分数表記になる。それをratで実数にすると(無限の連分数でないので元々実数だが)約分され、floatで小数にすると元の√2に限りなく近いことがわかる。
蛇足だが、cfdisrepしないと実数や小数に変換できない。

√3についてもやってみる。

方程式を解いてみる。solveで方程式を解析的に解ける。

解けない場合は、xの解にxが混ざったりする。(例は省略)

allrootsを使うと、多項式の方程式を数値的に計算できる。

指数の方程式を渡すと文句を言われる。

find_rootを使うと、(おそらくニュートン法による)近似計算で、解が1つだけ含まれる範囲の解を計算できる。

x^3=3^xのe~100の間の解は当然3なのだが、誤差が出た。


さて、せっかくフリーソフトなので、参考資料も無料のものを、と思って、図書館でMaximaの本を探したら、無かった。その代わり、あるMathematicaの本で、わずか十数行のコードでドラゴンカーブ(ドラゴン曲線)を描くMathematicaのコードを見つけた。マンデルブロー集合とかペアノ曲線とかのフラクタル図形の簡単なやつだ。Maximaの練習とMathematicaとの比較を兼ねて、その本を参考にMaximaで似たようなコードを書いてみた。

dragon(segments_):= block([p, q, r, s],
 p: map(
  lambda([x], [
   [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
   [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
  ]), segments_),
 q: map(
  lambda([x], [
   {[x[1][1], x[1][2]], [x[1][2], x[1][3]]},
   {[x[2][1], x[2][2]], [x[2][2], x[2][3]]}
  ]), p),
 r: map(lambda([x], listify(x)), flatten(q)),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(r)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 r
)$

p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$
for i:1 thru 10 do p:dragon(p)$
gnuplotが使える環境でこれを実行すると、以下のようなドラゴンカーブが描かれる。TeXmacsは必要ない。順に、L字型に配置した2本の線から始めて、全体を90°回転させて先端に繋げる、というのを1回、2回、…10回繰り返したものである。
アルゴリズムとしては、全ての線をL字に折り曲げる、というのを繰り返せばドラゴンカーブと同じ形になるというのを使っている。必要最小限にコードを簡略化して難解になったし、同様のアルゴリズムの説明は多くあると思うので、ここではコードの説明は省く。

続きを読む "Maximaでドラゴンカーブ" »

2008年08月13日

Maximaで微分方程式

せっかくMaxima使い始めたので、慣れるためにさらにMaximaのマニュアルを読み進めて遊んでると、Maximaは微分方程式も解けることがわかった。難しそうな響きを擁して実際難しい、物理学の基本ツールである。一瞬気が進まなかったが、大したものだと思ったので、敬意を表して試してみた。
すると、微分方程式が見えた気分になるようなグラフが出せることがわかり、面白かったので、ここに書き記すことにした。

まず、コマンドラインで微分させてみる。diff()が微分演算だ。

例によってTeXmacsを使って出力を整形している。(%i??)に続く青い部分が入力した部分、(%o??)に続く黒い部分がMaximaの出力だ。diff()の2つ目の引数は何で微分するか、3つ目の引数(省略可)は何回微分するかを意味する。

次に、desolve()を使ってシンプルな微分方程式\frac{d^2}{dx^2}y+\frac{d}{dx}y+y=0を解かせてみる。

desolve()の2つ目の引数は、どの関数について解くかを示す。diffの前の'(アポストロフィー)は、そのLisp関数を評価(実行&展開)せず、そのまま渡すことを指示するものである(desolve()の場合はこれを省略して評価してしまっても良いようだが、後述のode2()では不可)。y(x)の(x)はyがxの関数であることを示し、これを省略するとdy/dxが0になってしまう。
atvalue()を使って、y,y'の初期値を与えてから再度解かせてみる。

入力行の行末に;の代わりに$を付けると、計算結果を出力しない。
念のため、元の方程式に当てはまるかどうか、検算してみる。

%は直前の計算結果を表す。%o7は見るからに0なのだが、なぜか計算してくれなかったので、ratsimp()で整理させた。

常微分方程式(未知関数が1つだけの微分方程式)はode2()を使って解く方法もある。ode2()を使う場合、初期値は一般解の計算結果に対してic1(),ic2(),bc2()などを使って与える。同じ方程式を解かせてみる。

ode2()の2つ目の引数は常微分方程式の未知関数、3つ目の引数はその関数が従属する変数だ。上の場合は、引数からyがxの関数であることがわかるので、yをy(x)と書かなくても良い。

desolve()とode2()の能力を比較するため、解くには少し難しそうな微分方程式xy+\frac{d}{dx}y=0を解かせてみる。

ode2()は難なく解いた。念のため検算しておく。

それに対して、同じものをdesolve()に解かせると、下のような奇妙な式になった。

おそらく、ラプラス変換してx,yについて解いてラプラス逆変換したもの、という意味なのだろうが、これでは使えないだろう。これでも使える解なのかどうかを一応確認してみる。

やはり解けていない、使えない解であることがわかる。常微分方程式を解くのはdesolve()よりode2()の方が優れているようだ。

しかし、desolve()は複数の未知関数を含む連立微分方程式も解ける。
\left\{ \begin{array}{lll} \frac{d}{dx}f(x) & = & \frac{d}{dx}g(x)+\sin{x} \\ \frac{d^2}{dx^2}g(x) & = & \frac{d}{dx}f(x)-\cos{x} \end{array} \right.

e1とe2が方程式、f(x)とg(x)が未知関数だ。desolve()には、上のように、方程式と未知関数をそれぞれリスト形式で与える。

次に、plotdf()を使って微分方程式のグラフを出してみる。plotdf()に与える微分方程式は、dy/dx=F(x,y)または[dx/dt, dy/dt]=[F(x,y), G(x,y)]に限定されている。

load("plotdf");
plotdf(x);
これを実行すると、下のような画面が開く。(FreeBSD 6.xではTcl/Tkをインストールしていても「wishが無い」というエラーが出たが、wishという名のwish?.?へのリンクを作ればOKだった)

たくさんある矢印は、その地点(x,y)におけるdy/dxの向きと大きさを表している。このグラフ上のどこかをクリックすると、その座標を初期値とするグラフが描かれる。下の図は(0,0)辺りをクリックしたもの。

さらに色々な所をクリックする。

次の2つの図は、

plotdf(sin(x));
を実行したものである。

xとyを含む式も試してみる。

plotdf(sin(x)+cos(y));

放物線を横にした、x=y2を出そうとした。両辺をxで微分すると1=2y(dy/dx)なので、dy/dx=1/(2y)である。よりシンプルに、dy/dx=1/yで試そうとして、

plotdf(1/y);
とすると、原因がよくわからないが、Tclのエラーになった。しかし、次のように小さなオフセットを与えると、何か描画された。2つ目以降の引数に色々なオプションを指定しているが、これらが無くてもエラーにはならない。
plotdf(10^-10+1/y,[trajectory_at,0.5,1],[nsteps,200],[xradius,3],[yradius,6],[xcenter,3],[ycenter,0]);

trajectory_atで初期値を(0.5,1)として描画したものだが、左側が何かおかしい。さらに何ヶ所かクリックしてみると、

y=0の近辺でおかしくなるようだ。1/(y+⊿y)-1/yの計算でもしているのだろうか。

x=cos(t), y=sin(t)のグラフは円になるので、[dx/dt, dy/dt]=[-y,x] (=[-sin(t), cos(t)])のグラフも円になる(図は省略)。それよりもっと場を複雑にさせようと考えて、dx/dtとdy/dtを周期関数にしてみた。次のグラフは[dx/dt, dy/dt]=[sin(y),sin(x)]である。

plotdf([sin(y),sin(x)]);

色を変えながら初期値をクリックした。セル状のグラフが無限に続く。
dx/dtとdy/dtを入れ替えると、次のようなグラフになる。
plotdf([sin(x),sin(y)]);

こんなのも作ってみた。

plotdf([
10^-10 - sin(atan(y/x))*cos(%pi/4*sin(8*atan(y/x))) - cos(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))),
10^-10 - sin(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))) + cos(atan(y/x))*cos(%pi/4*sin(8*atan(y/x)))
]);

[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せない(dx/dtにtが入るとか)方程式も、desolve()で解けるものなら
、いつものplot2d()でグラフにできる。
\left\{\begin{array}{l}\frac{dx}{dt}=3x-7y\\ \frac{dy}{dt}=5x-4y \end{array}\right.
(この方程式は[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せるが、plotdf()ではおかしなグラフになる)

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,1);
Z:desolve([diff(x(t),t)=3*x(t)-7*y(t), diff(y(t),t)=5*x(t)-4*y(t)],[x(t),y(t)]);
plot2d([parametric, x(t), y(t), [t,0,6]],[nticks,200]),Z;

アステロイド曲線\left\{\begin{array}{l} x''(t)=2y'(t)-3x(t)\\ y''(t)=-2x'(t)-3y(t) \end{array}\right.(実はx=cos(t)^3, y=sin(t)^3から逆算したもの)を描いてみる。

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,0);
atvalue(diff(x(t),t),t=0,0);
atvalue(diff(y(t),t),t=0,0);
Z: desolve([
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
],[x(t),y(t)]);
plot2d([parametric,x(t),y(t), [t,0,2*%pi]],[nticks,200]),Z;

[parametric, ...]の部分を、それらのリストにすることにより、複数のグラフを1枚に描ける。

eq: [
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
];
atvalue(x(t), t=0, 1);
atvalue(y(t), t=0, 0);
atvalue(diff(x(t), t), t=0, 0);
atvalue(diff(y(t), t), t=0, 0);
Z1: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,2);
Z2: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,3);
Z3: desolve(eq, [x(t), y(t)]);
L1: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z1;
L2: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z2;
L3: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z3;
plot2d([L1,L2,L3], [nticks,200]);

続きを読む "Maximaで微分方程式" »

2008年08月15日

Maximaでマンデルブロー図形

お盆なので、特にお盆とは関係ないが、ついでに定番のマンデルブロー集合も描いておく。

f(c,z):=c-z^2$
density(x,y):= block([c,z,d,i],
c:x+y*%i,
z:c,
d:0,
for i:0 while (abs(z)<2 and d<40) do (z:f(c,z),d:d+1),
d
)$
plot3d(density,[x,-0.5,1.5],[y,-1.5,1.5],[grid,100,100],
[gnuplot_pm3d,true],[gnuplot_preamble,"set view map;unset surface;set xrange [-0.5:1.5];set yrange [-1.5:1.5];"]);

続きを読む "Maximaでマンデルブロー図形" »

2008年11月09日

統計学復習メモ1: t分布とχ2分布

高校時代は確率が得意中の得意だった私が、大学の統計学の講義にはすぐついていけなくなった。高3の時は、授業で順列と組み合わせの公式を覚えた途端、ほとんど全ての問題が難なく解けた。周りには理解につまずく人がいたが、私にはなぜ理解できないのかがわからなかった。大学2年の時に受けた統計学の講義も、半年くらいは割と楽についていってたが、ある時点から突然理解できなくなった。それは、平坦な道に突如現れた高い壁だった。どちらかというと興味があった科目なので、教科書をじっくりと何度も読んで理解しようとしたが、特に数式や日本語が難しい訳でもないのに、全然頭に入らなかった。なぜ理解できないのかがわからない不思議な感覚だった。
「推定・検定」という名前の関門だった。

結局、大学で統計学の講義を3回くらい取って、毎回1から学んだが、いつも確率分布までは理解できるのに、「推定・検定」で手も足も出ず挫折した。3回目の統計学で悟ったのは、t分布とカイ2乗分布が障害ということだった。

それから幾年月、それ以来統計学に接する機会が無かったが、推定・検定のことがずっと気にかかっていた。大学時代に最も興味があった科目の1つでありながら全く私を寄せ付けず、私に最も挫折感を味わわせた推定・検定を克服しない限り、死んでも死に切れない。今更そんなものを勉強して何の役に立つのかと言われようが関係ないのである。

という訳で先々月ぐらいからリベンジを試みているが、やはり相性が悪い。何となく少し理解した気になっても、2週間もするとすっかり忘れてしまい、先に進まない。後半が推定・検定でそれで終わりの入門書チックな統計解析の本を買って読んでるが、後半のみもう4周目くらいである。私以外のこれを読む人には前半と後半、確率分布と推定・検定の間に落差は無いのだろうか。


無駄話はこの辺にして、勉強したことをメモする。

t分布とは、xを平均μ、分散σ2の正規分布N(μ,σ2)に従う標本、(xの上にバー)を標本平均、nを標本の数、sを標準偏差(不偏分散の平方根)とした時、
t=¥frac{¥overline{x}-¥mu}{s/¥sqrt{n}}
というtが従う分布であり、
χ2分布とは、
¥chi^2=¥sum_{i=1}^n ¥left(¥frac{x_i-¥mu_i}{¥sigma_i}¥right)^2
というχ2が従う分布である、というようなことがどこにでも書いてある。
入門書で無ければ、さらにこれらの分布の確率密度関数が、ガンマ関数を使った容赦ない鮮烈な式で書かれる。
これが推定・検定の本題の前に出てくるから厄介だ。

どちらもにわかにはイメージが掴みにくい数式であり、なぜこういうものが必要かを説明される前に、まずこれを覚えて、と言われるのは結構辛いと思う。
それらの分布の意味を理解する前に、使い道を1つくらい知った方が、最初の理解の足掛かりになるし、トータルとして速く理解できるのではないだろうか。


という訳で、t分布の使用例として、母平均の区間推定を考えてみる。
大勢の人が受けたテストの点が正規分布に従うとする。その中の10人の点数(標本)が
70, 60, 40, 80, 40, 60, 30, 80, 50, 90
だとすると、全体の平均点は95%の確率でどの範囲に入ると推測できるだろうか。

ここで、その範囲が(標本平均±t*√(標本分散÷標本数))で求められるというのが、区間推定のありがたい教えである。

上の標本の場合、平均が60、分散が400、自由度(n-1)=9で信頼区間が95%となるtが2.262なので、推定区間は60±2.262√(400/10) ≒ 60±14.3ということになる。標本の平均は60点だが、95%の確率で全体の平均がこの範囲に入るということだ。

標本平均の分散が(標本分散÷標本数)になるというのは定理としてあるので、t分布のtの値というのは、標本の標準偏差を何倍すれば推定区間の幅になるのかを教えてくれるものと考えられる。


χ2分布の使用例として、ばらつきの片側検定を考えてみる。
プログラマーがそれぞれ2,3,3,4,8人の5つのソフトハウスにソフトウェアの開発を委託してるとし、あるプロジェクトにおいてそれぞれのソフトハウスが発生させたバグの数が15,18,14,19,34だったとする。バグの数がプログラマーの数に比例するという仮定は間違ってる方の5%に入る(95%の確率で間違ってる)だろうか。

ここで、((測定値−期待値)の2乗÷期待値)の和、という値がχ2分布に従うので、χ2分布を見ればその発生確率がわかる、というのが、χ2検定のマジカルかつミラクルな教えである。

上の標本の場合、バグが計100件でプログラマーが計20人なので、1人当たり5件、ソフトハウス別では10,15,15,20,40件となるのが期待値である。従って、測定値のχ2の値は、
(15-10)2/10 + (18-15)2/15 + (14-15)2/15 + (19-20)2/20 + (34-40)2/40 ≒ 4.1
である。自由度(5-1)=4で棄却域(有意水準)が5%となる、χ2分布におけるχ2の値は9.49なので、測定値のχ2の値はそれより小さく、95%の確率で起こり得る範囲内であることがわかる。
従って、バグの数がプログラマーの数に比例するという仮定は間違ってるとは言えない。

χ2分布のχ2の値というのは、同時に起こらない複数の事象の発生回数と各々の回数の期待値がわかる場合に、その期待結果とのずれがどれくらいの確率で起こるものかを調べるのに使える、発生回数のばらつきの尺度だと考えられる。

続きを読む "統計学復習メモ1: t分布とχ2分布" »

2008年11月16日

統計学復習メモ2: t分布と不偏分散

t分布を、正規分布に従う標本の平均が従う分布である、という感じに捉えて、なんとなく理解した気になると、1つ疑問が湧いた。正規分布に従う複数の値の和はやはり正規分布に従うのではなかったか?それなら、値の和を変数の数で割ったもの、つまり平均もt分布でなく正規分布に従うはずではないのか?
x1,x2をそれぞれ正規分布N(μ112),N(μ222)に従う確率変数とすると、x1+x2はN(μ121222)に従う、つまり和の平均は平均の和、和の分散も分散の和となる、というのを正規分布の性質の1つとして教えられている。x1+x2が正規分布に従うなら(x1+x2)/2も正規分布(この場合N(\frac{\mu_1+\mu_2}{2},\frac{\sigma_1^2+\sigma_2^2}{4}))に従うはずである。
というか、正規分布N(μ,σ2)に従うn個の標本の平均はN(μ,σ2/n)に従う、つまりn個の平均を取ると分散が1/nになる、という定理を覚えていないと、t分布の式を読むのが難しいと思う。それからしてもやはり標本の平均は正規分布に従うはずである。

そんなことで混乱したのは私だけだろうか?

混乱の原因は、標本平均、標本分散と母平均、毋分散をごっちゃにしたことだった。母平均をμ、毋分散をσ2、標本をx1...xn、標本平均をm=\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i、標本分散をs_n^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2とすると、mはあくまでN(μ,σ2/n)に従うのであって、N(m,Sn2/n)に従うのではないのである。つまり、μやσが不明な場合はmが従う正規分布も不明なのである。
そこで、mやSn2からμの範囲を推定するのにt分布が出てくる。
(※本エントリーではシステムの都合上、標本平均をmと書いたりxバー(xの上にバー)と書いたりするが、両者は同じものである)

…という感じに混乱を解決するのにも私には結構時間がかかったのだが、その間に私の理解を妨げたものの1つが、不偏分散の概念である。
標本分散が
s_n^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2
であるのに対して、母平均の区間推定に必要になる不偏分散は
s_{n-1}^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2
と教わる。nで割るかn-1で割るかの違いだけである。なぜn-1で割るのかというと、「自由度がn-1だから」と教わる。そんなのは丸覚えすれば引っかかる必要の無いことであるが、私の脳はそこで「???」となり先に進まなくなった。

「推定値の不偏分散の自由度が(n-1)である」理由は、ばらつきの指標の元となる
Σ(xi-m)2
のmにm=1/n Σxiという縛りがあるため、mとx1...xn-1が決まるとxnが決まるので、この2乗和がn-1個分しか動かないから、と説明されることが多い。極端な例としては、n=2の時、本当の平均μがわかってれば
Σ(xi-μ)2
はx1の分もx2の分もばらつきに応じて大きくなるのに対して、
\sum_{i=1}^{n}(x_i-m)^2 = \sum_{i=1}^{n}(x_i-\frac{x_1+x_2}{2})^2
は計算すると
(x1-x2)2/4
となり、標本が2個あるのに、その間の差1個分しかばらつきの値として反映されない。同様の理由により、標本をn個としてもn-1個分しかばらつきの値として加算されないのである。

その説明によって、私は一瞬目から鱗が落ちた気がしたが、まだ納得できなかった。もし標本の数を母集団の大きさと同じにすれば、つまり母集団の全てをサンプルとして拾えば、mが本当の平均になるから、自由度がn-1でも分散は(1/n)Σ(xn-m)2ではないか?
…ここから先はまだ疑問が解けていないが、おそらく推定値の分散と実際のデータの分散は同じ尺度では測れないということのような予感がする。

ところで、2008/11/16時点の日本語のWikipediaの分散の項に、不偏分散の計算、というのが載っている。曰く、
\frac{1}{n} \sum_{i=1}^{n}(x_i - \mu)^2

= \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2+\frac{1}{n}\frac{1}{n} \sum_{i=1}^{n}(x_i-\mu)^2
になるので、
\frac{1}{n} \sum_{i=1}^{n}(x_i - \mu)^2= \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2
となるとのことである。これだ!と思って計算を追ってみると、途中の
\frac{1}{n} \sum_{i=1}^{n}(x_i - \mu)^2
= \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2+\frac{1}{n^2} \{ \sum_{i=1}^{n}(x_i-\mu)^2 +2\sum_{i\neq j}^{n}(x_i-\mu)(x_j-\mu) \}
までは確認できたが、その次の
= \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2+\frac{1}{n}\frac{1}{n} \sum_{i=1}^{n}(x_i-\mu)^2
はわからなかった。というか、少なくともn=2として展開してみる限り、そうはならない。=が≒の間違いなのかも知れないが、よくわからない。

結局、私は次のように理解した。上のWikipediaの計算と同じ方針で、x_i=(x_i-\bar{x})+\bar{x}として両辺の分散を取ると、分散の加法性から
\sigma^2=Var(X-\bar{x})+\frac{\sigma^2}{n}
となる。Var(...)の部分は標本分散なので、これを毋分散σ2について解くと、
\sigma^2 = \frac{n}{n-1}Var(X-\bar{x}) = \frac{n}{n-1}\frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2
= \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2
となる。これが標本から推測される毋分散であり、毋分散と同じ尺度で使える不偏分散ということである。

続きを読む "統計学復習メモ2: t分布と不偏分散" »

2008年11月29日

統計学復習メモ3: 視聴率調査と標本数

TVの視聴率ってどうやって測ってるんだろう?というのはよく聞かれる問いであるし、モニターとして選ばれたいくらかの数の一般家庭に視聴率計測用の機器を取り付けて測っているのだ、というのはよく言われる答えである。そのモニターの数は、一例では関東・関西で600世帯、その他の地域は200世帯、とも言われる。

たったそれだけ?ある番組を20人が見てれば視聴率10%?学校の1クラスより少ない。それではたまたまその番組が好きな人が2人多めや少なめにモニターに入ってれば、すぐ1%くらい変わってしまうじゃないか、そんな値を信用できるのか、そんな値に意味はあるのか?
…とか書くと、貴様本当に大学を理系で卒業したのか?と問われてしまいそうであるが、それでもやはり反射的には疑問に思ってしまうことを自白したい。

このいわゆる標本調査、標本がいくらぐらいだとどれくらいの信頼性があるのだろうか。また、誤差をある程度に抑えるには、どれくらいの標本数が必要なのだろうか。


これを考えるのにも、区間推定の式が使える。母比率の推定区間は[標本比率±t*√(標本比率*(1−標本比率)÷標本数)]で、このtは標準正規分布に従う、と教えられている。数式で書くと、
I \in p\pm t\sqrt{\frac{p(1-p)}{n}} …(1)
(Iは推定区間の下限と上限、pは標本比率、nは標本数)である。信頼度を90%とするとt=1.645であり、n=200で、ある番組を見てたのが20人だった場合、p=0.1ということなので、母比率の推定区間は0.1±1.645*√(0.1*0.9/200) ≒ 0.1±0.035、すなわち6.5%〜13.5%となる。結構大雑把だ。信頼度を95%とするとt=1.96なので、5.8%〜14.2%ともっと幅が広くなる。
この式の±以降の部分が推定区間の幅なので、標本比率を固定すると、誤差をどれくらいにするためにはどれくらいの標本数が必要かを求めることができる。最大誤差をeとすると、
e=t¥sqrt{¥frac{p(1-p)}{n}}
なので、これをnについて解くと、
n=¥left(¥frac{t}{e}¥right)^2p(1-p)
である。例えば、視聴率10%の近辺で、信頼度を90%で視聴率の誤差を1%以内にするためには、n=(1.645/0.01)2*0.1*0.9=2435と求まるので、2400人以上のモニターが必要ということになる。誤差が2%で約600人、3%で約270人である。視聴率が20%近辺だと、さらにハードルが上がり、約4300人で誤差1%、約1100人で誤差2%、約500人で誤差3%となる。視聴率が50%に近づくほど、必要なモニターが増えることがわかる。


ところで、最初の母比率の推定区間の式(1)は、どうしてそうなるのだろうか。

以前のエントリーで、正規分布に従う値の母平均の推定区間が[標本平均±t*√(標本分散÷標本数)]で与えられる、と書いた。数式で書くと、
I \in \bar{X}\pm t\sqrt{\frac{s^2}{n}}
(Iは推定区間の下限と上限、Xは標本、s2はXの不偏分散、nは標本数)という感じであるが、この平方根部分の中身は要するに標本平均の分散ということであるので、ちょっと乱暴であるが、母平均に限らず、標本から求められる推定値からの推定区間が[推定値±t*√(推定値の分散)]で与えられる、と理解していいと思う。
なお、tはt分布に従うが、自由度が大きい場合はt分布が標準正規分布に近似できることになっている。

この推定値が比率の場合、その比率の分散が何になるかが問題だが、これは二項分布の定理から求められる。発生確率がpの場合、n回の試行で発生する回数Xの期待値はE(X)=np、分散はV(X)=np(1-p)だと二項分布の定理で確定しているので、n個の標本中の確率pで存在する何かの個数Xの期待値もE(X)=npであり、分散もV(X)=np(1-p)であることから、存在確率X/nの期待値はE(X/n)=p、分散はV(X/n)=p(1-p)/nだとわかる。これを使って、[推定値±t*√(推定値の分散)]=E(X/n)±t√V(X/n)=(1)となる訳だ。

続きを読む "統計学復習メモ3: 視聴率調査と標本数" »

2008年11月30日

統計学復習メモ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回しか合格しない確率だ。

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

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

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

2008年12月06日

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

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

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

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

2008年12月16日

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

以前のメモではカイ2乗検定の例としてばらつきの検定(適合度の検定)をやってみたが、もう1つのカイ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乗検定が使える。χ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)である。

以上を踏まえて、もう1つ例題をやってみる。原理はわかったので、手計算チックな計算は卒業して、次はMaximaに組み込みのカイ2乗検定関数を使おう、と思って関数を探すと、手元の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%でも、2つの属性が独立である仮説は棄却されないことになる。

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

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%とすると、である。

続きを読む "統計学復習メモ6: χ2検定で独立性の検定" »

2009年01月24日

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

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

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

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

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

ファイル名コード量コメント量
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

統計学で相関と言えば、まず思い付くのは相関係数である。とりあえず、ヘッダファイルとソースファイルに分けて、コード量とコメント量の相関係数を求めてみる。
よく使われる相関係数は、共分散÷それぞれの分散の幾何平均、というもので、式で書くと、2つのパラメーターを持つ標本が(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であった。

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

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

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


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

nr
50.878
100.632
150.514
200.444

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

続きを読む "統計学復習メモ7: t検定で独立性の検定" »

2009年02月06日

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

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

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

一次回帰直線とは、パラメトリックな2変量の関係を直線で近似するものなので、2次元のデータを1次元で近似することに相当する。2変量をx,yとすると、標本(x1,y1)〜(xn,yn)の関係を直線y=a+bx(a,bは定数)で近似するので、その直線との誤差を無視して(x,y)=(0,a)+(1,b)tと置くと、2変量の標本(xi,yi)が1変量tiで表せるのである。
主成分分析の第一主成分だけを使うことを考えると、多変量の関係を1次元で近似することになるから、k次元の標本(x1i,x2i,…,xki)を(C1,C2,…,Ck)+(p1,p2,…,pk)tiに近似すると、k個の変量を持つ標本(x1i,x2i,…,xki)が1変量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変量の標本を2変量で近似することになる。

一次回帰直線の最小2乗法による求め方をおさらいする。標本(xi,yi)と直線y=a+bxとの誤差はyi-a-bxiであるが、この誤差の2乗和を最小にするようにa,bを決めたい。誤差の2乗和をEとすると、E=Σ(yi-a-bxi)2で、このEはaの2次関数であり、bの2次関数でもあるから、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アニメが開きます)



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


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


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


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


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

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

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

続きを読む "統計学復習メモ8: 回帰直線と主成分分析" »

2009年02月08日

Maximaでの固有ベクトル計算について

筆者のMaxima 5.17+sbclの環境だと、標準のeigenパッケージを使って固有ベクトルを求めようとすると、エラーが出て求まらないことが多い。例えば、

load(eigen);
m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
eivects(m);

とやると、異様に時間がかかる上、
algsys failure: the eigenvector(s) for the1th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the2th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the3th eigenvalue will be missing.

というメッセージに続いて、虚数混じりのおぞましい値が固有値として出てきて、固有ベクトルは出力されない。そもそも固有値が間違っている。

固有ベクトルを英語でeigenvectorと言うことをたまたま覚えてたので、Maximaのヘルプでeigenvectorを検索して、出てきたのがeigenパッケージのeigenvectors()/eivects()だったので、それを使ったのだが、それが間違いだったようだ。Maximaのヘルプのeigenvectorsの項にも、固有値が正しく求められないことがあると、はっきり書いてある。対処方法も書いてあるが、とても面倒である。

インターネットで調べてみると、Maximaのメーリングリストのアーカイブに辿り着いた。どうやら2007年4月くらいからこの問題は原因も含めて知られているが、まだ解決していないらしい。

数式を含まない、純粋な数値行列なら、他の方法があることがわかった。
(1) eigens_by_jacobi()
標準関数として存在する。対称行列に限るという制約があるが、例えば共分散行列の固有ベクトルを求めるなら何の問題も無い。先に言ってくれという感じである。

(%i1) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i2) eigens_by_jacobi(m);
(%o2) [[- .7932213800442011, 1.180843735860609, .5243776441835921],
[ .6946378770918785 .5975383492309946 .4005323219261915 ]
[ ]
[ - .2001963037214474 .6953728293011738 - .6902014693159894 ]]
[ ]
[ - .6909411405362675 .3992549930407618 .6026572747810012 ]

eivects()と同様、固有値の行列と固有ベクトルの行列が一緒に返される。

(2) LAPACKパッケージのdgeev()
昔からある最強の線形演算ライプラリであるLAPACKのLisp版がMaximaに同梱されている。最初に使用する時にコンパイルがなされるが、これが環境依存であり、Lispの実行環境によってはコンパイルに失敗してしまうらしい。無事にコンパイルに成功すると、dgeev()で最強の計算結果が得られる。

(%i1) load(lapack);
(%i2) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i3) dgeev(m,true,true);
(%o3) [[- .7932213800442016, 1.180843735860611, .5243776441835922],
[ .6946378770918787 - .5975383492309948 .4005323219261915 ]
[ ]
[ - .2001963037214475 - .6953728293011737 - .6902014693159894 ],
[ ]
[ - .6909411405362672 - .3992549930407617 .6026572747810011 ]

この例は対称行列だが、もちろん非対称の行列でも問題ない。

筆者はgclとsbclで試したが、sbclではコンパイルが終了しなかったが、暴走したMaximaを強制終了すると、次からはコンパイル無しでロードできるようになった。(テスト段階で暴走していた?)

なお、root権限でインストールしたMaximaのLAPACKパッケージをコンパイルするには、root権限が必要である。(UNIX系のOSなら1回rootでMaximaを起動してload(lapack);とやれば良い)

2009年02月15日

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

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

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

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


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

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

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

とすると、W1〜Wmが第1主成分〜第m主成分、λ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 ]

固有値(寄与率)は絶対値が一番大きいのが1つ目、二番目に大きいのが2つ目なので、第一主成分の単位ベクトルは(0.714, 0.488, 0.502)、第二主成分の単位ベクトルは(0.698, -0.551, -0.457)である。3つ目の固有値は小さいので、3つ目の固有ベクトルは無視することにする。
今回はたまたま大きい順に並んでいるが、固有値は絶対値の大きさ順に並ぶとは限らないので注意が必要である。
もし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 ]

2つ目の固有ベクトルの符号が逆だが、次に求める主成分の大きさの符号が逆になるだけなので、どちらでも良い。

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

W: submatrix(W,3); /* 弟三主成分を削除 */
Y: X.W;

で求まる。これによって得られるYをグラフにすると、次のようになる。

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


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

このXが特異値分解によって
X=U¥Sigma V^T
(Tは転置)と3つの行列(U=左特異行列、Σ=特異値行列、V=右特異行列)の積になる時、Uが主成分の単位ベクトルを縦に並べたもの、Σの対角成分(特異値)が[n×寄与率]の平方根となり、標本の主成分の大きさ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に、先程と大体同じ主成分の単位ベクトルがあることと、3つ目の特異値が小さく、寄与率が低いことがわかる。

Maximaへの入力(続き)

S: zeromatrix(2,100); /* 行列Σ(3x100)の上2行分 */
S[1,1]: s[1]; /* 対角成分1 */
S[2,2]: s[2]; /* 対角成分2 */
Y: S.V;

これによって得られるYをグラフにすると、次のようになる。

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


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

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

続きを読む "統計学復習メモ9: 主成分分析の計算方法" »

2009年02月21日

ラグランジュの未定乗数法

統計学を復習していて、どうやらそれを使わないと非常に難しいらしい問題に突き当たったので、ラグランジュの未定乗数法を復習する。
その響きで一瞬怯んでしまうが、前提条件にさえ気を付ければ使う分には何も難しくない、ロピタルの定理のような便利技である。

どこにでも書いてあるので、わざわざここに書くまでもないが、ラグランジュの未定乗数法とは、x1...xnの関数fが束縛条件g=0の下で極値を持つ時、その極点のx1...xnは、f-λg(λは仮置きの定数)をx1...xnで偏微分したものが0になる値として求めることができるというものである。
数式で書くと、

である。ロピタルの定理同様、考えてみるとミラクルである。

例題として、2008年2月1日にフジテレビから放送された「たけしのコマネチ大学数学科」の第75回の問題を取り上げる。

問:

このような円柱と三角錐を組み合わせた鉛筆形の立体があり、体積が決まると、円柱の半径と円柱・円錐の高さは鉛筆形の表面積が最小になるように決まるとする。円錐部分の稜線の長さ(頂点から円錐面に沿って引いた線の長さ)が3の時、円柱部分の半径はいくらか。

解答:
円柱の半径をx、円柱の高さをy、円錐の稜線の長さをzとすると、表面積Sと体積Vは、中学校で習った公式を思い出すと
S=2¥pi x^2+2¥pi xy+¥pi xz
V=¥pi x^2y+¥frac{¥pi x^2}{3}¥sqrt{z^2-x^2}
と求まる。

まずは失敗の例から。普通に考えると、例えばVの式をyについて解いて、Sの式に代入してyを消して、
S=¥frac{2¥left( 3V-¥pi {x}^{2}¥sqrt{{z}^{2}-{x}^{2}}¥right) }{3x}+¥pi xz+¥pi {x}^{2}
このSを最小にするx,zを求めれば良いのだが、よく見ると1つの式に対して変数がx,z,Vと3つある上、例えばxについて偏微分すると
¥frac{¥partial S}{¥partial x}=-¥frac{¥sqrt{{z}^{2}-{x}^{2}}¥,¥left( 6¥,V-3¥,¥pi ¥,{x}^{2}¥,z-6¥,¥pi ¥,{x}^{3}¥right) +2¥,¥pi ¥,{x}^{2}¥,{z}^{2}-4¥,¥pi ¥,{x}^{4}}{3¥,{x}^{2}¥,¥sqrt{{z}^{2}-{x}^{2}}}
(右辺はMaximaの出力)となり、いわゆる強靭な計算力が必要な状況に陥る。ちなみに、Maximaでz,Vについても偏微分させてそれぞれ=0の連立方程式を解かせると、"system too complicated."というエラーを出してgive upする。

そこでラグランジュの未定乗数法である。
Vという条件の下でSの最小値を求める問題なので、

である。(正確にはVは(V-体積)だが、偏微分で体積が消えるので不要)
yについての式を実際に偏微分すると

となり、λ=2/xが得られる。次に、zについての式を偏微分してλを消去すると、

となり、これにz=3を代入すると、x=√5と求まる。
Vもyもわからないままで良い。なんということであろうか。

続きを読む "ラグランジュの未定乗数法" »

2009年03月14日

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

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

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

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

とすると、Xに対応する主成分軸上の第1主成分Yは
Y=P^TX
であり、そのYを元の座標系に戻したものX~は
¥tilde{X}=PY=PP^TX
である。このことは、高校で習った一次変換を思い出してやってみるとわかる。このX~が、Xを第1主成分の軸上に射影したものであり、これとXとの距離が、最小にしたい誤差ということになる。その誤差Eを、Xを直交座標とした場合の距離の2乗とすると、
E=(¥tilde{X}-X)^T(¥tilde{X}-X)
であり、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の2次関数であるEの最小値ををp12+p22=1の条件下で求めるので、「ラグランジュの未定乗数法」というやつを使う。すると、G=p12+p12-1と置くと、

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

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

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

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

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


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

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

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

として、第1主成分

の分散
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の2次関数であるVの最大値をp12+p22+...+pm2=1の条件下で求めるので、やはりラグランジュの未定乗数法を使う。すると、G=p12+p22...+pm2-1と置くと、

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

すなわち

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

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

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

2009年03月20日

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

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

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

(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の正方行列であり、
X=U¥Sigma V^T...(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を特異値分解したものの一部として求められることがわかる。

続きを読む "統計学復習メモ11: なぜ特異値分解で主成分が求まるのか" »

2009年04月29日

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

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

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

出たE(¥hat{¥theta})=¥theta。学生時代にこのθの上の屋根の意味が理解できなかったことが、私が推定・検定の1手目でつまずいた原因の1つであることは間違いない。そもそもなぜ統計学では推定するパラメーターをθと書くのか。高校の三角関数の授業で習って以来、私にとってはθといえば角度に違いないのであり、パラメーターをθと表すのは、円周率を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の不偏推定量である。

続きを読む "統計学復習メモ12: 不偏推定量とは" »

2009年05月10日

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

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

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

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


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

不偏推定量であるが一致推定量でない極端な例としては、標本をX1...Xnとして、平均θの推定量θ^=X1やθ^=Xnがある。標本がいくら増えてもその中の1つしか使わなければ精度が上がらないので、当たり前である。標本を全て使っても一致推定量でない例として、こういうのを考えてみた。
¥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になる)ようなθ^のことである。確率の対数をθで偏微分して2乗して平均して逆数にして訳がわからないが、対数が出てくるのは、AとBとCが同時に起こる確率がP(A)×P(B)×P(C)なのでその対数がlogP(A)+logP(B)+logP(C)となるように、標本それぞれが同時に起こる確率が線形和になるためのもので、偏微分は飛ばして、2乗の平均は誤差の2乗和の意味で出てくるのだと思う。

上記の不等式を解く代わりに、
¥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
であり、標本分散がθの有効推定量である。
標本分散は不偏ではないので、正規分布の分散の有効推定量は不偏推定量ではないことがわかる。

続きを読む "統計学復習メモ13: 有効推定量と不偏性、一致性" »

2009年05月14日

統計学復習メモ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が独立(1つ前の標本によって発生率が変わったりしない)とすると、それらが同時に起こる確率はL(θ)=P(X1)×P(X2)×...×P(Xn)である。これが最大になるθ=θ^を求めるのである。さらにL(θ)が単峰(山頂=極大点が1つしかない)なら、そのθは
¥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となる。

続きを読む "統計学復習メモ14: ポアソン分布のパラメーターの最尤推定量" »

2009年05月16日

統計学復習メモ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が妥当だと思う。

2009年06月06日

Coupon collector's problem

問1:おまけとしてフィギュアか何かが入っているお菓子があるとする。おまけは全m種類あり、お菓子の購入時にはどれが入っているかわからないとし、m種類は同じ確率で入っているとする。このお菓子をn回買うと、おまけが全種類揃う確率はいくらか?


これはTV東京の「ミネルヴァの法則」という昔の番組で出てきた問題のアレンジである。番組での問題は、おまけは10種類だと、100回買うと99%以上の確率で全種類揃うか(○/×)、という感じだったと思う。これの正解発表の時に、全種類揃う確率が99%になるのは66回であると、計算方法抜きで出てきた。

大学を卒業し、数式の類から遠ざかって10年以上、数学力はすっかり鈍っていたとはいえ、大学は理系卒で、高校時代は確率が得意中の得意だった身である。たかがテレビに出てくる確率の問題が解けないはずが無いと思って、その66を計算で求めようとしたが、やり方がわからず、その番組を見た日からのべ10時間以上考えて降参した。

三十半ばで再び数学に目覚めて1年、紙の上で数式をぐりぐりする右手を取り戻して、やっと解けた。


答1:
P¥left( m,n¥right) =¥sum_{k=1}^{m}{¥left( -1¥right) }^{m-k}¥,¥pmatrix{m¥cr k}¥,{¥left( ¥frac{k}{m}¥right) }^{n}
または
P¥left( m,n¥right) =¥sum_{k=0}^{m-1}{¥left( -1¥right) }^{k}¥,¥pmatrix{m¥cr k}¥,{¥left( ¥frac{m-k}{m}¥right) }^{n}
(これらはΣを展開すると同じ式、(m k)が縦に並んでるのは組み合わせ(mCkと書くのと同じ))


解説を試みる。

n個買うと、結果のパターンはmn通りである。この内、m種類揃っているのは、(m-1)種類以下しか無い場合を引いたものである。
(m-1)種類しか無いパターンは、(無い1種類のパターンmC1)×((m-1)種をn回選ぶパターン(m-1)n)、(m-2)種類しか無いパターンは、(無い2種類のパターンmC2)×((m-2)種をn回選ぶパターン(m-2)n)、…のように計算できると簡単なのだが、そうはならない。例えば(m-1)nパターンの中にも(m-2)種類以下しか無いパターンが含まれてしまうからだ。

この問題には「包除原理」(Inclusion-exclusion principle)というのが使えて、正しくは、(m-1)種類以下しか無いパターンの数は
¥pmatrix{m¥cr 1}(m-1)^n - ¥pmatrix{m¥cr 2}(m-2)^n + ¥pmatrix{m¥cr 3}(m-3)^n - ¥ldots - (-1)^{m-1}¥pmatrix{m¥cr m-1}1^n
であり、従ってm種類あるパターンの数はm^n(=mC0 m^n)からこれを引いたもの、
¥pmatrix{m¥cr 0}m^n - ¥pmatrix{m¥cr 1}(m-1)^n + ¥pmatrix{m¥cr 2}(m-2)^n - ¥ldots + (-1)^{m-1}¥pmatrix{m¥cr m-1}1^n
という+と−が交互に出てくる和になる。これを全てのパターンの数m^nで割ってΣ記号で整理したのが上記の答である。

この問題に「包除原理」が使えることを理解するには、勘が冴えてれば何も要らないかも知れないが、勘が鈍い日は次のように考えてはどうだろうか。
包除原理とは、例えばn人の人がいて、A,B,Cの3つの特徴を調べ、1人の人が複数の特徴を持つことができる時、いずれの特徴も持たない人の数が、
全体の人数−(Aの人数+Bの人数+Cの人数)+(AかつBの人数+BかつCの人数+CかつAの人数)−(AかつBかつCの人数)
となるようなことを言うものである。なぜそのようになるかというと、

この図の白い部分が求める人数だが、全体から1つの特徴を満たす人数を引くと2つの特徴を持つ人を2回引くことになり、引き過ぎた分を足すために2つの特徴を満たす人数を足すと3つの特徴を持つ人を3回足すことになり(m個の特徴を持つ人は(m-1)個の特徴を持つ人としてm回数えられるため)、以下同様に全ての特徴を満たす人数まで相殺する必要があるからである。
話を戻して、m種類のおまけが揃うパターンの数を求めるには、おまけ1が欠けてる場合をA、おまけ2が欠けてる場合をB、おまけ3が欠けてる場合をC、…と考えると、AとBが同時に起こることもあるから、A,B,…は上の図のように重なるので、どれも欠けていない場合、すなわち上の図の白い部分の数を数えようとすると、必然的に包除原理のような計算になるのである。


さて、現実的には、何個買った時にコンプリートしてる確率がいくらになるかより、大体何個買ったらコンプリートできるか、すなわち、コンプリートするまでの回数の期待値の方が興味があるのではなかろうか。

問2:おまけが全種類揃うまでのnの期待値は?

この問題は、"Coupon collector's problem"と呼ばれるらしい。日本語訳は不明だが、「食玩問題」とか「クーポン収集問題」とか呼ばれることがあるようだ。
確率分布がP(X)の確率変数Xの期待値E(X)を求める公式はΣ(X P(X))であり、上のP(m,n)より、n回目にm個揃う確率がP(m,n)-P(m,n-1)なので、nの期待値E(N,m)(以降、E(m)と省略する)はΣ(n (P(m,n)-P(m,n-1)))である。コンピューターでそれなりに大きなnまで計算すると正しい値に近づくので、一応正しい式のようであるが、これを解こうとするとひどい目に遭う。
しかし、違う発想で解くと、かなり簡単に求められる。


答2:
E(m)=m¥left(¥frac{1}{1}+¥frac{1}{2}+¥frac{1}{3}+¥ldots+¥frac{1}{m}¥right)


k-1種類まで揃ってから、k種類目が得られるまでに、平均何回かかるかと考える。
まだ得られてないのは(m-k+1)種類なので、1回買ってまだ得られてないのが得られる確率pは、p=(m-k+1)/mである。当たりの確率がpの時、当たりを引くまでの回数は、幾何分布に従う。すなわち、x回目で当たる確率はP(x)=(1-p)^(x-1) pである。このxの平均はΣxP(x)であり、等比級数の和を求める要領で解くと、ΣxP(x)=1/pとなる。従って、k種類まで揃ってからk+1種類目が得られるまでの回数の期待値はm/(m-k+1)である。
E(m)はこれを1種類目からm種類目まで足したものだから、
E(m)=m/m + m/(m-1) + m/(m-2) + ... + m/(m-m+1) = m(1/1 + 1/2 + ... + 1/m)
となる。

続きを読む "Coupon collector's problem" »

2009年06月09日

37%の法則

問:n個の箱にそれぞれ異なる額のお金が入っており、どれか1つの箱のお金をもらえるとする。最大の金額は不明である。1つ箱を開けたら、そのお金をもらって終了する(残りの箱は開けない)か、そのお金をパスするかを決めないといけない。1度パスしたお金はもらえない。最高額を得る確率を最大にするには、どのように箱を選択すれば良いか。


前項と同じく、「ミネルヴァの法則」の問題のアレンジである。番組の問題では、n=20だった。正解は、最初の37%(=7個)の箱は開けても選択せず、次に出てきた、それまでの最高額より大きい金額を選択するというものである。

前項の問題と同様、たかがテレビに出てくる確率の問題が解けないはずが無いと思って、その37%を(中略)降参した。
この度、やっと解けた。


答:最初のn/e個(eは自然対数、小数点以下切り捨て)の箱の中身をパスして、次に出てきた、それまでの最高額より大きいものを選択する。


解説する。
考えられる手順は、最初のr個の箱の中身をパスして、次に出てきた、それまでの最高額より大きいものを選択する、という手順に帰着する。これ以外を考えても無駄であることは容易にわかる(注:ここでの問題では最高額以外は全て外れである)。従って、その適切なrを選ぶのが問題である。

r個をパスして最高額が得られる確率をP(r)とする。最高額がx番目(x>r)にあり、それを選択できる確率は、
・最高額がx番目にある確率=1/n
・1〜r番目の最大より大きな額が(r+1)〜(x-1)番目には無い確率=1〜(x-1)番目の最大が1〜r番目に出てくる確率=r/(x-1)
の積である。P(r)はそれを全てのxについて足し合わせたものなので、
P(r)=¥sum_{x=r+1}^n ¥frac{1}{n}¥frac{r}{x-1}
となる。

P(r)が最大になるrを求める前に、P(r+1)-P(r)を計算して、P(r)の増減を調べてみる。

なので、上から下を引いて、
P(r+1)-P(r)=¥frac{1}{n}(-1+¥frac{1}{r+1}+¥frac{1}{r+2}+¥cdots+¥frac{1}{n-1})...(1)
である。これはよく見るとrが大きくなると項が減る一方で、どこかで0より小さくなる。従って、P(r)が最大になるのは、P(r)が減り始める直前のr、すなわちP(r+1)-P(r)<0となるrである。(1)の右辺の括弧内に注目すると、
¥frac{1}{r+1}+¥frac{1}{r+2}+¥cdots+¥frac{1}{n-1}<1...(2)
となる最小のrとも言える。

さらに、nが十分大きいとして、(2)の左辺を計算し易い積分に近似する。区分求積法より、
¥sum_{k=1}^{n}¥frac{1}{n}f(¥frac{k}{n})¥simeq¥int_0^1f(x)dx
なので、f(x)=1/xとして積分範囲を調節すると、
¥sum_{k=r+1}^{n-1}¥frac{1}{k}¥simeq¥int_¥frac{r}{n}^1¥frac{1}{x}dx
と近似できる。右辺<1を計算すると、-log(r/n)<1であり、r/n>1/eすなわちr>n/eが求まる。n/eは整数ではなく、r>n/eの整数だとP(r)がちょっと減るので、n/eを超えない整数がP(r)を最大にするrである。(nが小さい時は(2)で検算した方が良いかも知れない)


37%とは、1/e(≒0.368)のことであった。
なお、この問題は(classical) secretary problemと呼ばれ、問題の解は37%の法則とか37%ルールとかと呼ばれるらしい。

続きを読む "37%の法則" »

2009年06月28日

秘書問題の考察(2) 期待最高順位による停止戦略

「秘書問題」(secretary problem)というのは、非常に奥深い問題で、数理科学の「最適停止理論」における1つの重要な分野となっており、未解決の問題も少なくないらしい。
前のエントリーに書いた、37%の法則が成り立つ問題は、秘書問題の最も単純なものである。秘書問題に倣って秘書を面接で選ぶという表現に変えると、前のエントリーの問題は、応募者の中の最高の人以外は外れという条件の問題であるが、現実の問題としては、できるだけ最高の人を選ぶことだけを考えればいいとは限らない。誰か1人を選ばないといけない場合、最後の応募者まで採用を見送ると、最後の人はどんな人でも選ばないといけない。その1つ前の、残り1人(まだ面接していない応募者が1人)の時点で、前回の37%の法則の戦略に従うと、今面接している人がこれまでの中の2位であっても、残りの1人が1位である可能性に期待して見送るのであるが、なるべくいい人を選びたい場合は、残り1人で暫定2位であれば、残りの1人が2位以内に入る確率は(応募者が十分多いとすると)極めて低いので、暫定2位の人を採用するのが当然である。
このことからも直感的にわかる通り、1位の人を選ぶ確率を最大化するのとなるべく上位の人を選ぶのとは一般に両立しない。なるべく上位の人を選ぶ、というのを、順位の期待値を最小化する、と定義すると、1位の人を選ぶ確率を最大化するのと順位の期待値を最小化するのを同時に満たす最適停止戦略は存在しないことが数学的に示されているらしい。

という訳で、今回は次のような問題を考える。
・n人の応募者から1人を選ぶ。
・1人ずつ面接し、その場で採用するかどうかを伝える。採用したら、残りの人は面接しない。採用しなければ、その人は2度と選べない。
・応募者の良し悪しは、既に面接した応募者と比べてどうかでしか判断できない。絶対評価で何点というのは無く、あの人と比べてどうか、しかわからない。また、優劣は必ず付けられ、必ず順位付けが可能である。
・結果としてなるべく順位の高い人が選べるようにする。

必ず誰か1人は選ばないといけないということと、最後の1点以外は、前のエントリーの問題と同じである。
この問題は既に解かれていて、nを大きくすると順位の期待値の最小が約3.87に収束する(期待順位がそれより小さい戦略は存在しない)という摩訶不思議な定理になっているのだが、あまりにも難しいので、今回は単純な方法で考えてみる。


xを残りの応募者の数(まだ面接していない人数)とする。
x=0(今面接している人が最後の応募者)の時、その人を選ぶしか無い。
x=1の時、今面接している人を含むn-1人を1位〜n-1位に並べ、残りの人がその間または1位の上、n-1位の下のどこかに入ると考えると、入る場所はnヶ所なので、これを0.5位、1.5位、…、(n-0.5)位とすると、残りの人の順位の期待値はn/2なので、今面接している人の順位がn/2より小さければ選ぶ。(n/2以下なら、としても良い)
x=2の時、同様にこれまでの1位〜n-2位に対して残りの2人が0.5位、1.5位、…、(n-2+0.5)位の(n-1)ヶ所のどこかに入ると考えると、残りの2人のいい方の順位の期待値は、1人目の順位をkとすると
(2人目がkよりいい確率)*(2人目の順位の期待値)+(1 - 2人目がkよりいい確率)*k
であり、1人目の順位がkである確率は1/(n-1)なので、
¥frac{1}{n-1}¥sum_{k=1}^{n-1}¥left¥{¥frac{k-1}{n-1}¥,¥frac{k}{2}+¥left(1-¥frac{k-1}{n-1}¥right)k¥right¥}-¥frac{1}{2}
¥frac{n}{3}-¥frac{1}{2}+¥frac{1}{6-¥frac{6}{n}}
と計算できるので、今面接している人の順位がそれより小さければ選ぶ。

x>2の時、残りのx人の最高順位の期待値を同様に求めるのは難しいので、代わりに一様分布に従うx個の標本の最小値を求めてみる。
標本が0-1の一様分布に従う時、n個の標本の最大値がyより小さい確率は、全ての標本がyより小さい確率だから、ynである。これは最大値がyになる確率の累積密度だから、確率密度関数はそれを微分してn yn-1である。従って、最大値の期待値は、
¥int_{0}^{1}y¥,ny^{n-1}dy=n¥int_{0}^{1}y^ndy=¥frac{n}{n+1}
である。標本を全て(1-標本)にすると、最小値の期待値は1 - n/(n+1) = 1/(n+1)ということになる。
これを使うと、x>2の時、残りのx人の最高順位の期待値は、0.5+(0〜(n-x)の値を取るx個の標本の最小値の期待値)=0.5+(n-x)/(x+1)と概算できるので、今面接している人の順位がこれより小さければその人を選べば良い、と考えられる。


上の計算に沿って、n=10で残りがx人の時、今面接している人が何位以内ならその人を選べば良いかを表にしてみる。

xrank
15
23.17
32.25
41.7
51.33
61.07
70.88

x=1の時(9人目)は暫定5位以内なら採用、x=2の時(8人目)は暫定3位以内なら採用、x=3の時は暫定2位以内なら採用すれば良く、x=7の時(3人目)は暫定1位でもパスすべしということになる。

n=20についてもやってみる。

xrank
110
26.5
34.75
43.7
53
62.5
72.13
81.83
91.6
101.41
111.25
121.12
131
140.9

x=1の時は10位以内なら採用、x=2の時は6位以内なら採用、(中略)x>13の時は1位でもパスすべしということになる。

続きを読む "秘書問題の考察(2) 期待最高順位による停止戦略" »

2009年07月09日

秘書問題の考察(3) 期待順位最小の後ろ向き帰納戦略

前回のエントリーの秘書問題をきちんと解く。

問:以下の条件で、結果としてなるべく順位の高い人が選べるようにするには、どのようにすれば良いか。
・n人の応募者から必ず1人を選ぶ。
・1人ずつ面接し、その場で採用するかどうかを伝える。採用したら、残りの人は面接しない。採用しなければ、その人は2度と選べない。
・応募者は、それまでに面接した人の中で、必ず順位付けされる。応募者の良し悪しは、その暫定順位でしか判断できない。


解:
まず、j番目の応募者の暫定順位(それまでのj人中の順位)がxの時、その応募者が全n人中の何位であるかの期待値を求める。全n人中の順位をyとすると、この期待値は、(順位×確率)の和なので、E(y)=ΣyP(y)のようになる。j人中x位の人がn人中y位である確率は、
1 ...... y ...... n
 ↓
1 ... x ... j
という対応を考えると、j人の組み合わせnCj通りの中で、yの左側(y-1)人の内(x-1)人がxの左側に入って、yの右側(n-y)人の内の(j-x)人がxの右側に入る確率なので、P(y)=(y-1)C(x-1) × (n-y)C(j-x) ÷ nCjである。これを使うと、yの期待値は
¥sum_{y=x}^{n-j+x}y¥frac{¥pmatrix{n-y¥cr j-x}¥,¥pmatrix{y-1¥cr x-1}}{¥pmatrix{n¥cr j}} = ¥frac{n+1}{j+1}¥,x
となるらしい。(参考文献(1)による。この変形は私にはできなかったが、このP(y)は「超幾何分布」の形をしてるので、その知識を使えば簡単な計算で求まるのだと思う。)このE(y)を、以後y(j,x,n)と書く。

次に、最後の応募者の場合から、具体的な手順を考える。j-1番目までの応募者を採用しない時の最終的な採用者の期待順位をR(j,n)とする。
n-1番目までの応募者を採用せず、n番目の応募者を面接する場合、R(n,n)=n番目の応募者の期待順位=1位〜n位の平均=(n+1)/2である。
n-1番目の応募者については、その人が暫定x位の時、その人の最終的な期待順位はy(n-1, x, n)であり、その人を選ばない時の期待順位はR(n,n)である。従って、y(n-1, x, n) ≦ R(n,n)ならその人を採用すれば良いということになる。n-1番目の人がある暫定x位である確率は1/(n-1)なので、min{a,b}をa,bの小さい方とすると、R(n-1,n) = 1/(n-1) * Σ min{ y(n-1, x, n), R(n,n) }と書ける。
これは、n-1番目の応募者に限らず、1≦j<nの全てのjについて同じである。すなわち、j番目の応募者が暫定x位の時、その人の最終的な期待順位はy(j, x, n)であり、その人を選ばない時の期待順位はR(j+1,n)である。従って、y(j, x, n) ≦ R(j+1,n)ならその人を採用すれば良いというのが答である。R(j,n) = 1/j Σ min{ y(j, x, n), R(n,n) }、yを展開すると
¥frac{1}{j}¥sum_{x=1}^{j}min¥left¥{¥frac{n+1}{j+1}x, R(j+1,n)¥right¥}
なので、R(j,n)は末端のR(n,n)から逆向きに帰納的に(backward inductionと言うらしい)計算できる。


n=10の場合、採用する人の期待順位を最小にするために、j番目の人が暫定何位以内ならその人を選べば良いかは、次のようになる。3列目は、その人まで達した時の最終的な期待順位である。


jthresholdR(j,10)
954.28
833.59
723.15
622.89
512.68
412.56
302.56

3人目まではパスし、5人目までは暫定1位なら、7人目までは暫定2位なら、8人目は暫定3位なら、9人目は暫定5位なら採用すれば良く、この戦略による期待順位は約2.56となる。

n=20の場合は次のようになる。


jthresholdR(j,20)
19108.01
1876.62
1755.70
1645.05
1534.56
1434.18
1323.89
1223.64
1123.46
1013.30
913.17
813.06
713.0021
613.0017
503.0017

5人目まではパス、10人目までは1位なら、13人目までは2位なら、15番目までは3位なら(中略)採用すれば良く、期待順位は約3.00となる。

参考文献
(1)「タイミングの数理 最適停止問題」穴太克則 著(朝倉書店)

続きを読む "秘書問題の考察(3) 期待順位最小の後ろ向き帰納戦略" »

2009年07月12日

ソフトウェアリリース問題

秘書問題」の参考文献として図書館から借りた本に、曲がりなりにはソフトウェアに関わっている筆者としては、読んどかないと気になって仕方なくなるに違いない、一際目を引く問題があったので、返却する前に一読した。
「ソフトウェアリリース問題」と呼ばれている、最適停止問題の1つだ。
単純化した概要をメモする。


問1:テストが未実施である、開発中のソフトウェアの全バグ数がBとし、n回目のテストで見つかるバグ数をX(n)とする。BとX(n)の確率分布は既知で、E(X(n))は単調減少とする。見つかったバグは、各回のテスト終了時にコスト0で直ちに直される(バグ修正のコストはテストのコストに含まれる)とする。テスト1回当たりの費用をCt、ソフトウェアリリース後に残されるバグ1件当たりのコストをCb、Ct<<Cbとする時、コストを最小にするには、何回テストしてリリースするのが最善か。


ソフトウェアの開発の仕事をしたことがある人なら1回は考えさせられたことがあるであろう、どこまでテストすればいいのかという問題である。「バグが0になるまでに決まってるだろ!」と思う人が居るかも知れないし、そのように上から言われてる人が居るかも知れない。筆者は過去にそんな理不尽な圧力を何度も目撃した。しかし、現実的にはバグが0であることを証明するのは不可能である。ごくごく一部の例外を除き、はるかに不可能であり、不可能オブ不可能ズであり、1秒でも考えるのが無駄なくらい十分に不可能なのである。バグが無いことが証明されてリリースされるソフトウェアなんて存在しないのであり、どれだけテストして修正しても、ソフトウェアにバグは必ず存在するのである。
現実社会できちんと真っ当なソフトウェア開発の仕事をする人にとっては、バグの数が何らかの確率モデルで仮定されることは至極当然であり、むしろ、そうしないと仕事にならない。


解1:n回テストした時のトータルコストC(n)はn*Ct + (B - ΣX(k))Cbであり、
C(n+1)-C(n) = Ct - X(n+1)Cb
であり、E(X(n+1))は単調減少なので、E(C(n+1)-C(n))は単調増加であり、どこかでテストのコストの平均が残バグのコストの平均を上回る。E(C(n+1)-C(n)) > 0を解くと、E(X(n+1)) < Ct/Cbなので、答えはE(X(n+1)) < Ct/Cbとなるnである。


問2:Bが平均λのポアソン分布に従い、X(n)がその時点のバグ数と確率pをパラメーターとする二項分布に従う時、問1の最適なnはいくらか。


ソフトウェアのバグ数は、ソフトウェアの規模に比例すると考えられることが多い。システムの規模が大きくなると単位コード量当たりの複雑さも増すので、比例するというのはおかしい、と考えるのも自然なことなのだが、現実にはソフトウェア開発の対価はコード量に比例して支払われることが多いため、ソフトウェアは冗長であってもなるべく単純に作られ、単位コード量当たりの複雑さは上がらないのである。コード量が減るように工夫して効率の良いコードを作ると逆に報酬が下がってしまうので、コードは最適化されずに単純で似たようなコードの積み重ねとなり、システムの複雑さがそのままコード量として表れるのである。
従って、ソフトウェアのバグ数が(規模×単位コード量当たりのバグ数)を平均とするポアソン分布に従うと仮定するのは、現実社会に居て違和感を感じない。
1回のシステムテストで見つかるバグの数は、1つ1つのバグについて確率pで見つかると考えると、バグ数×pを平均とする二項分布に従うと仮定するのは、至極妥当であろう。


解2:B(n)をn回目のテスト後のバグ数とすると、1回のテスト当たりに見つからずに残る割合の平均が(1-p)なので、B(n)は平均λ(1-p)nのポアソン分布に従う。X(n+1)の平均はB(n)*pなので、答えはE(X(n+1)) = λp(1-p)n < Ct/Cbとなるnである。


参考文献
「タイミングの数理 最適停止問題」穴太克則 著(朝倉書店)

続きを読む "ソフトウェアリリース問題" »

2009年07月19日

秘書問題の考察(4) 最小期待順位の極限値について

前回のエントリーで説明した、秘書問題(secretary problem)の期待順位を最小にする解に関して、もう少し追究する。


まず、前回は参考文献を鵜呑みにした、j人中x位の人がn人中何位であるかの期待値が
¥sum_{y=x}^{n-j+x}y¥frac{¥pmatrix{n-y¥cr j-x}¥,¥pmatrix{y-1¥cr x-1}}{¥pmatrix{n¥cr j}} = ¥frac{n+1}{j+1}¥,x
となること(左辺が右辺に等しくなること)を導出する。
これの左辺をy(j,x,n)と書く。j人中x位の人は、もう1人加えて(j+1)人とすると、x位か(x+1)位かのどちらかになる。(j+1)人目が上位x位に入ると(x+1)位、入らないとx位のままである。このことから、j人中x位の人の最終順位は(j+1)人中(x+1)位の人の最終順位と(j+1)人中x位の人の最終順位で表され、
y(j,x,n) = x/(j+1) y(j+1,x+1,n) + (j+1-x)/(j+1) y(j+1,x,n)
の関係が成り立つことがわかる。扱いやすいように、jを1つずらして、
y(j-1,x,n)=¥frac{x}{j}¥, y(j,x+1,n) + ¥frac{j-x}{j} y(j,x,n)
としておく。
y(n,x,n)=xは自明だから、そこからy(n-1,x,n)が求まり、再帰的に計算するとy(j,x,n)が求まる。
j=n-1については、
y(n-1,x,n) = ¥frac{x}{n}(x+1) + ¥frac{n-x}{n} x = ¥frac{n+1}{n}x
j=n-2については、
y(n-2,x,n) = ¥frac{x}{n-1}¥frac{n+1}{n}(x+1) + ¥frac{n-1-x}{n-1}¥frac{n+1}{n}x = ¥frac{n+1}{(n-1)n}¥left(x(x+1)+(n-1-x)x¥right)=¥frac{n+1}{(n-1)n}nx=¥frac{n+1}{n-1}x
なので、y(j,x,n) = (n+1)/(j+1) xではないかと予想できる。それを帰納法で確認する。
j=n-1, n-2については上のように成り立つ。あるjについて成り立つとすると、
y(j-1,x,n) = ¥frac{x}{j}¥frac{n+1}{j+1}(x+1) + ¥frac{j-x}{j}¥frac{n+1}{j+1}x = ¥frac{n+1}{j(j+1)}¥left(x(x+1)+(j-x)x¥right)=¥frac{n+1}{j(j+1)}(j+1)x=¥frac{n+1}{j}x
となり、(j-1)についても成り立つので、間違いなくy(j,x,n) = (n+1)/(j+1) xである。


次に、Chowという人の1964年の論文によるらしい、その筋では頻出の、筆者を秘書問題地獄に陥れたミラクルな定理を紹介する。
定理:期待順位最小の最適戦略による期待順位をV(n)とすると、
¥lim_{n¥rightarrow¥infty}V(n)=¥prod_{j=1}^{¥infty}¥left(¥frac{j+2}{j}¥right)^¥frac{1}{j+1}¥simeq3.8695
である。
つまり、応募者が何人に増えようが、千人だろうが1万人だろうが、平均的にその中の3.87位以上の人を選択できる採用戦略が存在するのである。その数字の小ささもさることながら、その数字がnを含まないある定数に収束するというのは驚きではなかろうか。

さらに、期待順位最小の最適戦略を、r(x,n)番目の応募者からは暫定x位以上の人を採用すると書くと、次の式が成り立つ。
¥lim_{n¥rightarrow¥infty}¥frac{r(x,n)}{n}=¥frac{1}{V(¥infty)} ¥prod_{j=1}^{x-1}¥left(¥frac{j+2}{j}¥right)^¥frac{1}{j+1}
r(x,n)/nは面接した応募者の割合だから、nが大きい時、全応募者のどれくらいの割合を面接した時点から何位以上の人を選べば良いかというのが、これで計算できる。x=1〜10について計算すると


xr(x,n)/n
10.2584
20.4476
30.564
40.6408
50.6949
60.735
70.7658
80.7903
90.8101
100.8265

となるので、最初の25.8%の人はパスし、25.8%の人が過ぎたら1位の人を、44.7%の人が過ぎたら2位以上の人を、56.4%の人が過ぎたら3位以上の人を(以下略)採用すれば良いということになる。(但し、これはnが小さい場合は正確な値との差が大きい。前回のエントリーでの計算結果と比較すると、n=20でもかなりずれている。)

参考文献
Optimal Stopping and Applications(最適停止理論の大家、Thomas S. Ferguson教授による解説)Chapter 2.

続きを読む "秘書問題の考察(4) 最小期待順位の極限値について" »

2009年08月02日

正20面体の頂点座標の求め方

正20面体の頂点座標を簡単に求める方法があることがわかった。


同じ長方形3枚をこのように組み合わせて、近い頂点を結ぶと、何か20面体ができる。これの長方形の長辺と短辺の比を変えていくと、2つの長方形から3点取った三角形がどこかで正三角形になる。

その時、20面体の30本の辺の長さが全て同じになるので、その時の3つの長方形の頂点が正20面体の頂点である。

長方形の短辺の長さを2、長辺の長さを2xとし、点A〜Eを

このように取ると、AE=x、ED=(x-1)なので、AD2=x2+(x-1)2であり、DC=1なので、AC2=AD2+DC2=x2+(x-1)2+1である。ABCが正三角形の時、AC=2なので、
x2+(x-1)2+1=4
が成り立ち、これを解くとx=(1+√5)/2(=黄金比)が求まる。

従って、Gを黄金比とすると、
(±1,±G,0)
(0,±1,±G)
(±G,0,±1)
の12点が、1つの正20面体の全ての頂点となる。

続きを読む "正20面体の頂点座標の求め方" »

2009年09月02日

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

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

筆者は大学で「検定」を習って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から遠い方が「差がある」ように統計量を定める方が便利であることが多いからである。例えば標準正規分布やt分布やカイ2乗分布に従う統計量は0から遠いと平均と差がある。
そのため、もし帰無仮説を「差がある」方に取ると、統計量の棄却域が0に近い方になり、範囲が狭くなる。標準正規分布やt分布のような山形の分布だと、

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

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

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

次の図は、カイ2乗検定(片側検定)の棄却域の例である。

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

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

続きを読む "統計学復習メモ16: 帰無仮説と棄却域" »

2009年10月03日

統計学復習メモ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の2乗和が最小になるようにAを求めるとすると、
¥sum_{i=1}^{n}e_i^2=E^TE=(Y-XA)^T(Y-XA)
であり、これはaiの2次関数なので、これを偏微分して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×共分散行列)の逆行列の要素)
という値が自由度n-p-1のt分布に従う、というのを使う。定数項a0については、siiの部分は次のような式になるが、同じようにt値が求められる。
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)))$
it値確率
01.1280.3025
14.5070.004074
20.12570.904
31.4670.1928
有意水準を0.05とすると、a1についてのみ、本当は0である(係数に意味が無い)という仮説が棄却され、有効ということになる。

続きを読む "統計学復習メモ17: 重回帰分析の計算結果の検定" »

2009年11月29日

余弦定理、正弦定理の導き方

今週、TVで「余弦定理」という言葉を耳にした。高校でそういう公式を習って名前だけは覚えていたが、中身は雰囲気すら思い出せなかった。
15年以上前にセンター試験を受けた時は完璧に暗記していたはずである。三角関数の加法定理は今でも身体が覚えていた。漢字と同じように書き順で覚えていて、右手が紙に書き始めると、何も考えなくても勝手に手が動いて、筆記体のような加法定理を書き終える。ちなみに我が筆記体では+とxは同じ形であった。
倍角の公式や和積・積和の公式や合成の公式(sinθ+cosθ=sin(θ+α)という変形)を加法定理から無事に導き出せて、正弦定理も何とか思い出すことができたのだが、余弦定理は影も形も思い出せなかった。
仕方なく、インターネットで調べた。

●余弦定理

三角形の3辺の長さをa,b,c、各辺の向かい側の角度をA,B,Cとすると、
a^2=b^2+c^2-2bc¥cos A
b^2=c^2+a^2-2ca¥cos B
c^2=a^2+b^2-2ab¥cos C
が成り立つ。

●正弦定理
同じく、
¥frac{a}{¥sin A}=¥frac{b}{¥sin B}=¥frac{c}{¥sin C}=2R
(Rは外接円の半径)が成り立つ。

そうそう、そういえばそんなのがあった。何でそういうのが成り立つのかがわからないまま、丸暗記した記憶がある。何でそうなるんだろう?と、何度も不思議に思った記憶もある。
これらを初めて教わってから20年、ついに導き方を理解した。

●余弦定理の導出

このようにaに垂線を引くと、
a=b¥cos C+c ¥cos B
の関係があることがわかる。これの両辺を2乗すると、
a^2=b^2¥cos^2 C+c^2¥cos^2 B+2bc¥cos B¥cos C
である。sin2+cos2=1やコサインの加法定理
¥cos(B+C)=¥cos B¥cos C-¥sin B¥sin C
を使って整理すると、

となる。A+B+C=π(180°)だから、
¥cos(B+C)=¥cos(¥pi-A)=-¥cos A
であり、b sinC = c sinB = aへの垂線の長さで末尾の2乗部分は0なので、余弦定理の式が求まる。

●正弦定理の導出

このように同じ外接円を持つ1辺がaの直角三角形を描くと、円周角の定理より、aの向かいの角の角度はやはりAである。また、同じく円周角の定理より、円周角が90°なら中心角が180°なので、斜辺は外接円の中心を通り、長さは直径=2Rに等しい。従って、2R sinA=aである。

続きを読む "余弦定理、正弦定理の導き方" »

2010年01月28日

統計学復習メモ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である。(完)

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

発生率が一定で、1回の試行の結果は発生するかしないかの2通りしか無い試行(ベルヌーイ試行)を複数回行った時の発生回数は二項分布に従う。そして試行回数が十分大きければ二項分布は正規分布で近似できる。試行回数をn、発生率をpとすると、発生回数は平均np、分散np(1-p)の正規分布に(近似的に)従う。観測される発生率(推定確率)は発生回数を試行回数で割ったものだから、平均p、分散p(1-p)/nの正規分布に(近似的に)従う。母分散が既知の場合、正規分布に従う標本の推定値は毋分散÷標本数を分散とする正規分布に従うとするのが区間推定の定跡である(母分散が未知の時は標本分散を使うのでt分布に従うとする)から、信頼区間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

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

2011年06月18日

n個のサイコロの目がk種類になる確率

次のような感じの確率を求めようとしたら、泥沼にはまった。

Q: サイコロをn回振った時、出る目がk種類になる確率は?

A:
k=1の確率は、n回全て同じ目で、その目は6通りあるから、6/(6n)である。
k=2の確率は、場合の数が(n回全てaかbの場合 − n回全てaまたはn回全てbの場合)×((a,b)の組み合わせ)だから、(2n-2)(6C2)/(6n)である。
k=3の確率は、場合の数が(n回全てaかbかcの場合 − n回全てaかbの場合 − n回全てbかcの場合 − n回全てcかaの場合 + n回全てaまたはn回全てbまたはn回全てcの場合)×((a,b,c)の組み合わせ)だから、(3n-3*2n+3)(6C3)/(6n)である。
つまり、k種類の値を固定した時の場合の数が包除原理に従っているので、答は
\frac{1}{6^n}\pmatrix{6 \cr k}\sum_{i=1}^k(-1)^{k-i}\pmatrix{k \cr i}i^n
である。(括弧内に2つの数字が縦に並んでるのは二項係数、(6 k)は6Ckと同じ)


サイコロの目の数6をmに一般化して整理すると、等しい確率でm通りの値を取るn個の変数の値がk種類となる確率P(m,n,k)は
\pmatrix{m \cr k}\sum_{i=1}^k(-1)^{k-i}\pmatrix{k \cr i}\left(\frac{i}{m}\right)^n
となる。

続きを読む "n個のサイコロの目がk種類になる確率" »

2011年07月30日

クロソイド曲線

x'=\cos\left(\frac{\pi t^2}{2}\right)
y'=\sin\left(\frac{\pi t^2}{2}\right)
というシンプルな数式で表されるそうなので、Maximaでその数式でグラフを描いてみようとしたが、残念ながら媒介変数を含む微分方程式をそのまま描画することはできないようだった。
微分方程式をグラフ化する関数としてplotdf()があるが、parametric plotができないようだ。
integrate()で定積分したものをplot2d()でparametric plotすればいいかと思ったが、残念ながら上記の関数のintegrate()は計算されなかった(erf混じりの非数値解になる)。

上記のx',y'を定積分したものはフレネル積分と呼ばれるそうで、Maximaにfresnel_c(), fresnel_s()があったので、今回はそれを使うことで良しとした。

●Maximaへの入力

plot2d([parametric, fresnel_c(t), fresnel_s(t), [t,-6,6]],[x,-1,1],[y,-1,1],[style,lines],[nticks,1000]);

●出力(ブラウザがSVGを表示できない場合→ PNG
Clothoid curve

車のハンドルを一定の速度で回していくとできる曲線だそうだ。

続きを読む "クロソイド曲線" »

2012年01月10日

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

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

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

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

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

以下のデータは、全て架空のものである。
筆者は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にしている。'*'のマークとその説明が出力されていないので、所属別では群間に有意な差が無く、このデータでは所属によって生産性にばらつきがあるとは言えないことになる。

もう1つ例を考えた。

一元配置(対応なし)の例(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:開発期間
1週間2週間1ヶ月
要因2:
曜日
月曜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
要因2:
座席の形態
大部屋 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 
このデータでは、交互作用を分離しなければ、開発プロセスには有意差が認められる(交互作用の効果を分離しない方が有意差が出る)ことがわかる。
ちなみに、繰り返しがないと、'+'でも'*'でも結果は変わらない。

もう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()の結果は、交互作用を分離しないと、仕様書のフォーマットについても同じ有意差は出ない(1つ前のデータ例とは逆に、交互作用を分離する方が有意差が出る)ことを示している。

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

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


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

二元配置(繰り返しあり、1要因対応あり)の例(1)
生産性
(steps/人月)
要因1:開発プロセス
WaterfallSpiralIncrementalTDD
要因2:
座席の形態
+ 人
作業場開発者
大部屋 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 
これは、データをよく見るとわかるのだが、小部屋に飛び抜けて成績が良い人が居ることに起因しており、作業空間の違いによる効果だと言うにはかなり不適切であろう。

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

二元配置(繰り返しあり、1要因対応あり)の例(2)
生産性
(steps/人月)
要因1:性別 + 人
性別男性女性
開発者ABCDEF
要因2:
業務分野
組み込み(アプリ) 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つの水準にしかなく、水準間でデータに対応がない
    例えば、水準を決める条件は同時に発生するので、水準毎に異なる被験者を選ぶ
もちろん、対応を決める要因は、個人や個体に限らず、標本に影響があることが既にわかっている条件であれば何でも良い。


二元配置(2要因対応あり)の標本のデータ構造は、1要因のみ対応ありのものよりも単純である。

二元配置(繰り返しあり、2要因対応あり)の例(1)
欠陥数
(steps/人月)
要因1:開発プロセス
WaterfallSpiralIncrementalTDD
要因2:
座席の形態
+ フェーズ
作業空間開発フェーズ
大部屋 仕様作成6673
全体設計4787
個別設計68105
テスト設計3264
実装6588
小部屋 仕様作成2635
全体設計5746
個別設計63104
テスト設計7564
実装12567
自宅 仕様作成5494
全体設計9296
個別設計7997
テスト設計5332
実装491011
要するに3次元なのである。ここでは表を立体的に書けないので、対応を決める「開発フェーズ」については列方向に展開している。
開発フェーズが違うと仕事の性質が全く違うので、フェーズの違いは当然ミスの数に影響する。
どんな開発プロセスにも、これくらいのフェーズは存在するので、各フェーズのデータはあり得る。また、どんな作業空間でも、一通りの開発をすれば全てのフェーズを経るので、各フェーズのデータが得られる。従って、異なる開発プロセスのデータ間でも異なる作業空間のデータ間でも、フェーズの違いによる対応が取れる。
> 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               
2要因対応ありのデータ構造は単純だが、その計算は、コンピューター任せでも複雑である。対応を決める要因の違いによる影響が、要因1との交互作用、要因2との交互作用、要因1と要因2と対応要因との交互作用と多岐に渡って分離して計算されるからである。
上記の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               

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

二元配置(繰り返しあり、2要因対応あり)の例(2)
定着期間
(ヶ月)
要因1:設計手法
構造化設計データ指向オブジェクト指向コンポーネント指向
要因2:
開発内容
+ 性別
開発内容性別
プラットフォーム 男性 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               
このデータでは、男女の違いによる対応ありとすると、設計手法の違いによる有意差が見られ、対応なしとすると、設計手法にも開発内容にもその違いによる差が見られないという結果になった。

続きを読む "統計学復習メモ19: 分散分析の種類" »

2015年08月10日

フィボナッチ数列の任意の項を高速に計算する方法

そんなプログラミングの問題を見かけて、反射的に思ったのは、一般項を用いることだった。
フィボナッチ数列の漸化式はan=an-1+an-2なので、一般項anを求めることは可能である。

3項を含む漸化式は、特性方程式を用いて解くのが基本であるが、やり方を忘れてしまったので、自己流で解いてみた。
a_{n+1} = a_{n} + a_{n-1}

a_{n+1} + Aa_n = (1+A)(a_n+ Aa_{n-1})
という形にできると、
(1+A)(a_n+ Aa_{n-1}) = (1+A)^2 (a_{n-1}+ Aa_{n-2}) = \cdots = (1+A)^n(a_1+Aa_0)
であることと、a1=1, a0=0であることを用いて
a_{n+1} + Aa_n = (1+A)^n
とできる。このAは、an-1の係数を見ると、(1+A)A=1となるAなので、A=(-1±√5)/2である。従って、

\left\{
\begin{array}{c}
a_{n+1} + \frac{-1+\sqrt{5}}{2} a_n = \left(1 + \frac{-1+\sqrt{5}}{2}\right)^n \\
a_{n+1} + \frac{-1-\sqrt{5}}{2} a_n = \left(1 + \frac{-1-\sqrt{5}}{2}\right)^n
\end{array}
\right.
であり、上の式から下の式を引いて整理すると、
a_n = \frac{1}{\sqrt{5}} \left( \left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n \right) と求まる。
念の為申し上げるが、このanはnが何であろうとも整数である。実にミラクルである。

√5のn乗なんてのがあるのでは、誤差が出まくるだろうとは思ったが、CやJavaでは浮動小数点の精度が低くて無理でも、PythonのmathパッケージにMathematicaやMaxima並の驚異の性能が無いとは限らないので、試しに

#!/usr/bin/python
from math import pow
from math import sqrt

def fib(n):
    return int(1/sqrt(5)*(pow((1+sqrt(5))/2,n)-pow((1-sqrt(5))/2,n)))

an_2 = 0  #a_{n-2}
an_1 = 0  #a_{n-1}
an = 1    #a_{n}
for i in range(1, 100):
    print "n =", i, "a_n =", an, "fib(n) =", fib(i)
    if fib(i) != an:
        print "error!"
    an_2 = an_1
    an_1 = an
    an = an_1 + an_2
とやってみた所、72項目で誤差が出た。
Maximaはn>200000でも正確な整数を弾き出す(下記入出力参照、下3桁が正しいことを確認)が、さすがにPythonでは無理だった。おそらく、Maximaは内部で√5を√5として保存しながら数式処理をしてるからであろう。Pythonの言語仕様上、pow()の結果は一度float値にしないといけないので、正確な整数値を出すには無限の精度が必要になってしまう。

Maximaへの入力とその出力
(%i1) 1/sqrt(5)*(((1+sqrt(5))/2)^n-((1-sqrt(5))/2)^n),n=234567,ratsimp;
(%o1) 179607932890115544888825400217[48962 digits]014673023583001136817746498978

MathematicaやMaximaを使って良いなら、一般項を用いる方法もあるが、今回の問題は、「複数言語選択可能」としながら、MathematicaもMaximaも選択肢に無かった。 フィボナッチ数列のn項目を高速に計算する方法は、他には思い付かなかったのだが、少し調べると、計算量がO(log N)になる簡単な方法があった。

M=\pmatrix{1 & 1 \cr 1 & 0} とすると、Mn-1の{1,1}要素がフィボナッチ数列のn項目だというのが使える。
さらに、M=M*Mという計算を繰り返すことにより、M^(2k)は計算量少なく求められるので、Mn-1をM^(2k)の積とすることによって、高速に計算できる。
例えば、
M5 = 0*M8 · 1*M4 · 0*M2 · 1*M
M10 = 1*M8 · 0*M4 · 1*M2 · 0*M
であるように、2進数にした時の1のビットに対応するM^(2k)を掛け合わせれば良いのである。

従って、例えば次のようにすれば、n=10000でも高速に計算される。(Pythonで扱える整数の上限はメモリが許す限り無限である)

#!/usr/bin/python

class multipliable_2x2matrix(object):
    def __init__(self, a, b, c, d):
        self.a = a;
        self.b = b;
        self.c = c;
        self.d = d;
    def mul(self, b):
        return multipliable_2x2matrix(
            self.a * b.a + self.b * b.c,
            self.a * b.b + self.b * b.d,
            self.c * b.a + self.d * b.c,
            self.c * b.b + self.d * b.d)
    
def fib(n):
    if n == 0: return 0
    m = multipliable_2x2matrix(1, 1, 1, 0)
    f = multipliable_2x2matrix(1, 0, 0, 1)
    n = n - 1
    while n > 0:
        if n & 1 == 1:
            f = f.mul(m)
        n = n >> 1
        m = m.mul(m)
    return f.a

2016年05月04日

"Sum and Product Puzzle"をPythonで解いてみた

昔、多分NetNewsのfj.sci.mathかどこかだと思うが、次のような問題が流れていた。

司会と A さん, B さんがいます.
司会 「今から Joker を除いたトランプ 52 枚から無作為に 2 枚を取り出し, A さんにはその二数の積を, B さんには和を教えます. その数から元の二数が何であるかを考えてください」
A, B 「はい」
司会 「(B さんにわからないよう, A さんに積を教える)」
司会 「(A さんにわからないよう, B さんに和を教える)」
A さん 「教えてもらった数からは元の二数はわかりません」
B さん 「私もわかりません. でも, A さんが『わからない』と言うのはわかっていました」
A さん 「それなら私はわかりました」
B さん 「それなら私もわかりました」
さて, A さん, B さんの聞いた数はそれぞれいくつでしょう?
これだけで答えが1つに特定できることが驚きであり、取ってあった。それをたまたま見つけたので、もう一度やってみた(解説は後述)。
一見どこから手を付ければ良いのかわからないが、「A さんが『わからない』と言うのはわかっていました」からBさんの数字を考えてみると、ほぼ答えまで絞られるので、解くのにそれほど時間はかからない。

NetNewsの元の投稿は探せなかったが、その代わりに、この問題の元ネタが見つかった。 元の問題は1969年に発表されたもので、その後に様々なバリエーションが作られ、一連の問題は"Sum and Product Puzzle"と呼ばれているらしい。参考文献[1]によると、1979年にMartin Gardnerが"the Impossible Puzzle"と呼んで紹介した頃の英語版の1つは次のようなものだったらしい。

Two numbers m and n are chosen such that 2 ≤ m ≤ n ≤ 99. Mr. S is told their sum and Mr. P is told their product. The following dialogue ensues:
Mr. P: I don’t know the numbers.
Mr. S: I knew you didn’t know. I don’t know either.
Mr. P: Now I know the numbers.
Mr. S: Now I know them too.

In view of the above dialogue, what are the numbers?
元の2数の範囲が2〜99なだけで、後は全く同じである。おそらく同じようにして解けるのだが、範囲が広い為に格段に面倒である。
どうせ、部分的にヒューリスティックに絞り込んでいるが、しらみ潰しに探しているのと大差無い。おそらく、整数の範囲を制限しているので、しらみ潰しに探す以外に解法は無いと思う。(例えば「2より大きいあらゆる偶数は素数の和で表せる」というゴールドバッハ予想も、整数の範囲を制限しているので、単独では使えない。m+nとしてあり得る188は7+181 = 31+157 = 37+151 = 61+127 = 79+109であり、99以下の素数の和ではないので、m,nが共に素数である組み合わせが無い。)
プログラミングの勘を取り戻す為のリハビリを兼ねて、プログラムを作ってしらみ潰しに探すことにした。

#!/usr/bin/python

MIN = 2
MAX = 99

def good_factors(p):
    """Generates pairs of factors of p"""
    return [(m,n) for m in range(MIN, MAX+1) for n in range(MIN, MAX+1) if m <= n and m*n == p]

def good_summands(s):
    """Generates pairs of summands of s"""
    return [(m,n) for m in range(MIN, MAX+1) for n in range(MIN, MAX+1) if m+n == s and m <= n]

def fact1(m, n):
    """Mr. P: I don't know the numbers."""
    return len(good_factors(m*n)) > 1

def fact2(m, n):
    """Mr. S: I knew you didn't know."""
    for (m_, n_) in good_summands(m+n):
        if not fact1(m_, n_):
            return False
    return True

# This is not necessary but just writing straightforward
def fact3(m, n):
    """Mr. S: I don't know either."""
    return len(good_summands(m+n)) > 1

def fact4(m, n):
    """Mr. P: Now I know the numbers."""
    fact2_compliers = [(m,n) for (m,n) in good_factors(m*n) if fact2(m,n)]
    return len(fact2_compliers) == 1

def fact5(m, n):
    """Mr. S: Now I know them too."""
    fact4_compliers = [(m,n) for (m,n) in good_summands(m+n) if fact4(m,n)]
    return len(fact4_compliers) == 1

result = [(m,n) for m in range(MIN, MAX+1) for n in range(MIN, MAX+1) if m <= n and fact1(m,n) and fact2(m,n) and fact3(m,n) and fact4(m,n) and fact5(m,n)]

print "result =", result

これを実行すると、(m,n)の正解である(4,13)が出力される。
MIN = 1, MAX = 13にすると、冒頭のトランプの問題の元の2数も計算できる。

但し、心持ち関数型プログラミングっぽくしてみたが、これだととても時間がかかる(Intel Core i5 1.6GHzのCPUでも2分以上)。次のように、good_factors()とgood_summands()が計算結果をキャッシュするように書き換えると、20倍くらい速くなった。

good_factors_dict = {}
def good_factors(p):
    """Generates pairs of factors of p"""
    if not good_factors_dict.has_key(p):
        good_factors_dict[p] = [(m,n) for m in range(MIN, MAX+1) for n in range(MIN, MAX+1) if m <= n and m*n == p]
    return good_factors_dict[p]

good_summands_dict = {}
def good_summands(s):
    """Generates pairs of summands of s"""
    if not good_summands_dict.has_key(s):
        good_summands_dict[s] = [(m,n) for m in range(MIN, MAX+1) for n in range(MIN, MAX+1) if m+n == s and m <= n]
    return good_summands_dict[s]

◆参考文献
[1] "Sum and Product in Dynamic Epistemic Logic"
H. P. van Ditmarsch, et al.
J Logic Computation (2008) 18 (4): 563-588
First published online: December 21, 2007
http://www.cs.otago.ac.nz/staffpriv/hans/sumpro/sumandproduct06.pdf

続きを読む ""Sum and Product Puzzle"をPythonで解いてみた" »

About 数学

ブログ「Weblog on mebius.tokaichiba.jp」のカテゴリ「数学」に投稿されたすべてのエントリーのアーカイブのページです。過去のものから新しいものへ順番に並んでいます。

前のカテゴリは将棋です。

次のカテゴリは社会です。

他にも多くのエントリーがあります。メインページアーカイブページも見てください。

Powered by
Movable Type 3.35