統計・確率のお勉強

統計学を中心に色々勉強するブログ

多項分布

多項分布

基本性質

確率関数
 \displaystyle
f(x_1,\ldots,x_k) = \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k}
期待値
 \displaystyle
E(X_i) = np_i
分散
 \displaystyle
V(X_i) = np_i(1-p_i)
共分散
 \displaystyle
Cov (X_i,X_j) = -np_i p_j

確率関数

1回の試行でk通りの可能な結果A_1,\ldots,A_kのいずれか1つのみが生じ、P(A_i) = p_i(i = 1,\ldots,k)とする。この試行を独立にn回繰り返したときに、A_iが生じる回数をX_iとするとき、X_1,\ldots,X_kの同時分布を多項分布といい、その確率関数は以下で与えられる。
$$
f(x_1,\ldots,x_k) = \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k}
$$

期待値

\begin{eqnarray}
E[X_i] &=& \sum_{i=1}^k x_i \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k} \\
&=& \sum_{i=1}^k \frac{n \cdot (n-1)!}{x_1!\cdots (x_i - 1)! \cdots x_k!}p_i \cdot p_1^{x_1} \cdots p_i^{x_i-1}\cdots p_k^{x_k} \\
&=& np_i \sum_{i=1}^k \frac{(n-1)!}{x_1!\cdots (x_i - 1)! \cdots x_k!} p_1^{x_1} \cdots p_i^{x_i-1}\cdots p_k^{x_k} \\
&=& np_i
\end{eqnarray}

分散

分散を求めるためにまず次の期待値を計算する。
\begin{eqnarray}
E[X_i(X_i - 1)] &=& \sum_{i=1}^k x_i(x_i - 1) \frac{n!}{x_1!\cdots x_k!} p_1^{x_1} \cdots p_k^{x_k} \\
&=& \sum_{i=1}^k \frac{n(n-1) \cdot (n-2)!}{x_1! \cdots (x_i-2)! \cdots x_k!} p_i^2 \cdot p_1^{x_1} \cdots p_i^{x_i-2} \cdots p_k^{x_k} \\
&=& n(n-1)p_i^2 \sum_{i=1}^k \frac{(n-2)!}{x_1! \cdots (x_i-2)! \cdots x_k!} p_i^2 \cdot p_1^{x_1} \cdots p_i^{x_i-2} \cdots p_k^{x_k} \\
&=& n(n-1)p_i^2
\end{eqnarray}
これと、V(X) = E(X(X-1)) + E(X) - E(X)^2であることを用いて分散を求める。
\begin{eqnarray}
V[X_i] &=& E[X_i(X_i-1)] + E[X_i] - E[X_i]^2 \\
&=& n(n-1)p_i^2 + np_i - n^2p_i^2 \\
&=& np_i(1-p_i)
\end{eqnarray}

共分散

期待値、分散のときと同様に計算することで、E(X_i X_j) = n(n-1)p_ip_jが得られるので、

 \displaystyle
\begin{eqnarray}
Cov(X_i, X_j) &=& E(X_i X_j) - E(X_i)E(X_j) \\
&=& n(n-1)p_i p_j - np_i \cdot np_j \\
&=& -np_i p_j
\end{eqnarray}

多項分布の例

サイコロをn回振ったときに、1の目がでる回数をX_1回。2の目がでる回数をX_2回。・・・6の目がでる回数をX_6回とする。また、それぞれの出る確率をp_i(i = 1,2,\ldots,6)とする。この時、確率ベクトル\boldsymbol{X} = (X_1,\ldots, X_6)'は、多項分布\mathrm{Multi}(n,\{p_i\})に従う。


参考文献

岩沢宏和(2012):『リスクを知るための確率・統計入門』,東京図書

「分かりやすい説明」という聞き手(読み手)の怠慢

「説明力」を求められる理系

最近、世の中から、特に理系に求められる能力として説明力がある。Wikipediaのような集合知、そこそこ専門的な知識でもググれば、その内容(理解できるかは別として)をすぐにでも確認することのできる時代に、専門用語をならべて偉そうに話す専門家の存在意義は薄れてきている(専門家が不要な訳ではない)。私自身、このことに対して、強く賛同している。分かりやすく説明できるということは、その分野に対して深い理解を持っていることの証明でもあるし、そのような説明ができる人間になりたいと考えてる。

 

氾濫する超入門書、内容の薄い○○でも分かる△△学

分かりやすさを求めるのは最早、世の中の大きな流れとなっていることが、書店に並ぶ本を見れば分かる。大きな書店に行けば、理工書がコーナーとしてあると思うが、そこに並ぶ本の中に少なからず、超入門系の本が置かれている。○○超入門!とか、文系でも分かる○○学!などといった本たちだ。目を引きやすい、いかにも優しそうなイメージの表紙に、「もう挫折しない!」みたいな帯がまかれていたりする。これらの本が一定数限りのあるならんでいるということは当然それなりに需要があるからならんでいるのだ。つまり、世の中がそのような本を求めているから、数少ない書店のスペースに、内容の薄さに対して、無駄に厚い超入門書が置かれるのである。

 

これらの本が、本当にわかりやすければいいのだ。数式を用いず、統計学ならば、各手法について、どうしてその手法なのか、どうやって使うのか、その結果は何を意味するのかを懇切丁寧に書いてくれるような書籍が本当の「分かりやすい」本なのではないかと思う。しかし、実際はどうだろう。特に私がよく見るのは統計関連のものだが、この手の本はだいたい、Excelの使い方を指南して終わる。なんか「回帰分析で予想!」みたいなことが書いてあって、Excelで一生懸命回帰分析の結果がでるまでの様子を、写真付き(これが、内容は薄く、本は厚くなる理由である)で説明するのだ。このことにどれほどの価値があるというのだろうか。正直、グーグルで検索したほうがよっぽどいい解説が出てくるだろう。そのレベルの無価値な本に、1600円、2000円なんて、どうして払えるだろうか?それ以上に、限りある書店スペースをそれらの本が占有していることが腹正しい限りだ。

「今」の入門書と「昔」の入門書の違い

分かりやすさが叫ばれるようになり、入門書のイメージも変わってきている。昔の入門書というのは「事前知識なくとも、読み進めることができる」というのが、「入門」の意味であった。例えば、多変量解析の入門書を読み進めるにあたって、行列の知識がかなり必要になってくる。昔の入門書というのは、このような必要な知識を付録として巻末に書いておいてくれる。また、当然、確率分布の定義から始まり、確率密度関数、期待値、分散共分散行列と、一つ一つ必要な知識を厳密に書かれているものであった。その分野の基本を詰めたようなものが入門書の意味するところであったのだ。それに対して、現在の入門書というのは、全く知識がない人に「なんとなくわかった気にさせる」ことが(超)入門書の意味するところになっている気がする。こういうものがあって、これはこう使います!といった、中身をそぎ落としすぎて、用語とExcelの使い方を説明する、「知ったか」を育成することを目的としているんじゃないかと疑いたくなるようなものだ。そもそも、学問というものはそれぞれの分野に一生をかける研究者がいる以上、本一冊読んだ程度で理解できるわけがない。

 

分かりやすさを求めた結果なくなるもの

分かりやすい説明をするためにはそのトレードオフとして「厳密性」が失われる。このことは多くの場合において重要ではないと考えられるが、それが失われることにより、多くの勘違いを生む結果となる。学問における定理、公式には、それが成り立つための「前提条件」というものが往々にして存在する。当然、説明する側はそのことは十分に理解しているだろうし、当然説明するときも、どこかしらで、その前提に触れていることだろう。しかし、聞き手側はどうだろうか。世の中の多くは、理系的な考え方に触れていない人が多く、そのような簡易な説明を求める層の多くはそのような人たちに占められると考えられる。そうなったときに、でかでかとプレゼンテーションで定理の内容が書かれていたのを見た聞き手が前提条件など覚えているだろうか。多くは覚えてなどいないだろう。関心があるのは、その結果にあるのだから。そうなってくると、当然前提条件が無視された定理の乱用が始まる。無価値な分析が世の中にはびこるようになる。統計等はその手法(Excelでなんかすれば結果出てくるし)よりも解釈のほうが難しい。そして、"Excelで分析した感を出した何か"を一生懸命発表するのだ。その結果が本当に正しいか、解釈があっているのかは二の次である。

 

何事も理解するためにはそれなりの努力が必要

分かりやすい説明を求めることは悪いことではない。ビジネスマンは研究者ではないので、少ない時間で概要をつかむ必要があるだろうし、研究者はそれにこたえられるほうがいいに決まっている。学生は、小難しい説明を永遠とされたらその授業に出る気など起きないし、勉強する意欲を失ってしまうだろう。

 

しかし、何を理解するにも、まっさらな状態で、なんとなく聞いていたのでは、いくらわかりやすい説明でも、理解できるわけがないのだ。その説明している人は数年から数十年をその分野にかけている人だろう。その内容の一部を短時間で、興味を持ってもらえるよう、もしくは、あなたが求めている部分だけに絞って説明してくれているのかもしれない。その努力に敬意を払い、自分なりに、理解する努力する必要がある。教えてもらうといった受け身思考ではよくない。学ぶんだという積極性が、たとえ、ガンガン質問にいくといった積極性までは持てなくとも、心持ちだけは積極でなくてはならない。そうでないと「分かりやすい」を追い求めた結果「内容の薄い」本や説明しか残らなくなるだろう(さすがにそんなことはないと思うが)し、そのような説明を聞いたり、その程度の本を読んだところで何も得るものなどない。そうならないためには、聞き手、読み手側にもそれなりの準備、努力が必要だと考える。

 

 

50年くらい前の本と現在本屋にならんでいる本を比べて、最近理論部分をしっかり扱う本が少ない(無いわけではない)気がして、ちょっとした愚痴でした。なんかすいません。

ベータ関数,ベータ分布

ベータ関数

定義

p,qが正の定数のとき、下記右辺の定積分を、p,qの関数と考え、ベータ関数と呼ぶ。

$$
B(p, q) = \int_0^1 x^{p-1} (1-x)^{q-1} dx \;\;\; (p,q > 0)
$$

ベータ関数とガンマ関数の関係

ベータ関数と、ガンマ関数の間には次の関係がある。
$$
B(p,q) = \frac{\Gamma (p) \Gamma (q)}{\Gamma (p + q)}
$$

ベータ分布

ベータ分布区間(0,1)上の確率分布であり、以下の確率密度関数によって定義される。

$$
f(x) = \left\{
\begin{array}{cc}
\frac{1}{B(p, q)} x^{p-1} (1-x)^{q-1} & (0 < x < 1) \\
0 & その他
\end{array}
\right.
$$
ベータ分布はBe(p, q)で表す。

以下X \sim Be(p,q)とする。

ベータ分布の平均,分散

\begin{eqnarray}
E[X] &=& \frac{p}{p + q} \\
V[X] &=& \frac{pq}{(p+q)^2 (p+q+1)}
\end{eqnarray}

導出

まずは平均について
\begin{eqnarray}
E[X] &=& \int_0^1 x \frac{1}{B(p,q)} x^{p-1} (1-x)^{q-1} dx \\
&=& \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \int_0^1 x^p (1-x)^{q-1} dx \\
&=& \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \frac{\Gamma(p+1)\Gamma(q)}{\Gamma(p+q+1)} \\
&=& \frac{\Gamma(p+q)}{\Gamma(p+q+1)} \frac{\Gamma(p+1)}{\Gamma(p)} \\
&=& \frac{p}{p+q}
\end{eqnarray}
分散を求めるに当たって、次のモーメントを求める。
\begin{eqnarray}
E[X(X-1)] &=& \int_0^1 x(x-1)\frac{1}{B(p,q)} x^{p-1} (1-x)^{q-1} dx \\
&=& -\frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \int_0^1 x^p (1-x)^q dx \\
&=& -\frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \frac{\Gamma(p+1)\Gamma(q+1)}{\Gamma(p+q+2)} \\
&=& -\frac{pq}{(p+q+1)(p+q)}
\end{eqnarray}
これより、V(X) = E(X(X-1)) + E(X) - E(X^2)から、分散は
\begin{eqnarray}
V[X] &=& E[X(X-1)] + E[X] - E[X]^2 \\
&=& -\frac{pq}{(p+q+1)(p+q)} + \frac{p}{p+q} - (\frac{p}{p+q})^2 \\
&=& \frac{pq}{(p+q)^2(p+q+1)}
\end{eqnarray}

モーメント

ベータ関数のk次モーメントを求める。
\begin{eqnarray}
E[X^k] &=& \int_0^1 x^k \frac{1}{B(p,q)}x^{p-1}(1-x)^{q-1} dx \\
&=& \frac{1}{B(p,q)} \int_0^1 x^{p+q-1}(1-x)^{q-1} dx \\
&=& \frac{B(p+k, q)}{B(p,q)} \\
&=& \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \frac{\Gamma(p+k)\Gamma(q)}{\Gamma(p+q+k)} \\
&=& \frac{\Gamma(p+q)\Gamma(p+k)}{\Gamma(p)\Gamma(p+q+k)}
\end{eqnarray}

p,qが正の整数のときは

 \displaystyle
E [X^k] = \frac{(p+q-1)! (p+k-1)!}{(p-1)!(p+q+k-1)!}
と書くことができる。

ベータ分布の密度関数のグラフ

ベータ分布のグラフは、パラメータごとに以下のようになる。
f:id:doratai:20170528192802p:plain
pythonで描画。以下のサイトをおおいに参考にした。
【Python】scipyとmatplotlibでベータ関数を描画 - 歩いたら休め
ソースコードは以下。

import numpy as np
import scipy.stats
from matplotlib import pyplot as plt

x = np.linspace(0, 1, 1000)
plt.xlim(0,1)
plt.ylim(0,5)
plt.xlabel(r"$x$", fontsize=20, fontname='serif')
plt.ylabel(r"$f(x; p,q)$", fontsize=20,fontname='serif')
plt.title("PDF of Beta Distribution")
params = [[3,9],[6,6],[9,3],[1,1],[1,5],[5,1],[16,16]]
for param in params:
    rv = scipy.stats.beta(param[0], param[1])
    y = rv.pdf(x)
    plt.plot(x,y,'-',lw=2,label=param)
    plt.legend(bbox_to_anchor=(1.05,1), loc='best', borderaxespad=0) #凡例を枠外表示
plt.show()
ベータ関数の特徴

ベータ関数の特徴として、上記のグラフを見ればわかると思うが、パラメータp,qの値によっていろいろな形を取る、ということがある。例えば(p,q) = (1,1)のときは一様分布になっていることがグラフからもわかる。

ベータ分布とベイズ

ベータ分布がよく出てくるのは、ベイズ統計の分野である。ベイズ統計では、事前分布と、事後分布というものを考えるが、その時の事前分布としてベータ分布はよく使われる。グラフでも示した通り、ベータ分布のグラフはパラメータによって非常に柔軟に形を変えることができる。また、グラフとパラメータの対応関係をよく見て欲しい。(p,q) = (3,9)の時、グラフは左に偏り、(p,q) = (9,3)の時、グラフは右に偏っている。つまり、p,qの比がそのままグラフに表現されるのである。このことは、確率的主観を表現する際に都合が良い。予想が6-4であるとして、期待値を0.6とするならば、「私の確信」をp:q = 6:4のベータ分布で表現できるのである。ベイズ主義は確率を「ある事象をどれくらいできるか」の指標と解釈しているため、主観的な確率というものが非常に重要になってくる。その主観的確率を表現する際にベータ分布は非常に都合の良い分布なのだ。

参考文献

松原望,縄田和満,中井検裕(1991):『統計学入門(基礎統計学I)』,東京大学出版会
日本統計学会(2013):『日本統計学会公式認定 統計検定1級対応 統計学』,東京図書

一般逆行列の定義と存在

大学教養レベルで扱う線形代数では、逆行列は「正則行列(非特異行列)」である必要があり、\mathrm{rank}がフルランクであることが逆行列を持つ必要十分条件であった。しかし、行列が特異(逆行列を持たない、フルランクでない)である場合でも、逆行列を持つように、逆行列を拡張した、一般逆行列というものが存在する。統計学の中でも多変量解析などの分野では行列を多用するため、行列の話題というのは非常関心の高いものになる。そのような分野における、一般逆行列の利用は今や日常茶飯事らしいので、定義くらいは知っておきたい。

定義

m \times n行列\boldsymbol{A}の一般逆行列(generalized inverse)とは
$$
\boldsymbol{AGA} = \boldsymbol{A}
$$
を満たす、任意のn \times m行列\boldsymbol{G}のことである。
※一般逆行列の他に擬似逆行列、条件付き逆行列という用語で呼ばれることも多い
実際、\boldsymbol{A}が非特異である場合は\boldsymbol{G} = \boldsymbol{A}^{-1}であるので、\boldsymbol{AGA} = \boldsymbol{AA}^{-1}\boldsymbol{A} = \boldsymbol{A}となっている。

一般逆行列の存在

次に気になるのが、この「一般逆行列」が存在するのかどうかである。結論としては

あらゆる行列は少なくとも1つの一般逆行列を持つ。

これは次の定理で証明される。

定理

\boldsymbol{B}を最大列階数のm \times r行列。\boldsymbol{T}を最大行階数のr \times n行列とする。この時、\boldsymbol{B}は左逆行列\boldsymbol{L}を持ち、\boldsymbol{T}は右逆行列\boldsymbol{R}を持つ。そして、\boldsymbol{RL}\boldsymbol{BT}の一般逆行列である。

証明

\boldsymbol{B}が左逆行列\boldsymbol{L}を持ち、\boldsymbol{T}が右逆行列\boldsymbol{R}を持つことは既知としよう。この時、一般逆行列の定義から
$$
\boldsymbol{BT}(\boldsymbol{RL})\boldsymbol{BT} = \boldsymbol{B}(\boldsymbol{TR})(\boldsymbol{LB})\boldsymbol{T} = \boldsymbol{BT}
$$
である。すなわち、\boldsymbol{RL}\boldsymbol{BT}の一般逆行列である。
今、任意のm \times n行列\boldsymbol{A}を考える。\boldsymbol{A} = \boldsymbol{0}であるならば、明らかに任意のn \times m行列は\boldsymbol{A}の一般逆行列である。\boldsymbol{A} \neq \boldsymbol{0}であるならば、\boldsymbol{A} = \boldsymbol{BT}を満たす最大列階数の行列\boldsymbol{B}と最大列階数の行列\boldsymbol{T}が存在する。したがって、「あらゆる行列は少なくとも1つの一般逆行列を持つ」という結論を得る。


この証明で既知として用いている部分に関しては参考文献を参照ください。

参考文献

David A. Harville,(監訳)伊里正夫(2012) : 『統計のための行列代数(上)』,丸善出版

無限のパラドクス〜数学から見た無限論の系譜〜を読んで

理工系の新書として有名なレーベルの「BLUE BACKS」の本で、足立恒雄先生の著書、「無限のパラドクス」を読みました。

 

この本は、現代の無限論にたどり着くまでの歴史的経緯について非常に簡潔に分かりやすく書かれており、高校レベルの数学の知識があれば難なく読むことが出来ると思います。

 

無限論とは別に読んでいて感動したこと(ライプニッツの夢)

この本の中盤で、ライプニッツ(1646〜1716)の話が出てきます。この本によると、論理学の分野で「AならばB」のように、一般命題に文字をわりあてるのはライプニッツが始めたことだそうです。ライプニッツが夢見たのは「全ての論証を記号化し、算術問題に還元すること」。これを実現すると、裁判などで、データを投入するだけで後は計算機が計算して公平な判決を下してくれる。このような事を夢見ていたそうです。

少し前ならば、なんて夢物語なのだろうと私も思いました。しかし、人工知能、バックデータ解析といった用語がトレンドな現代、ライプニッツの夢見たことは現実味を帯びてきています。将来的にコンピュータに仕事を奪われると予想される職業に「弁護士」が入っているのですから。

 

人工知能が判決を下す。なんか、少し前のノイタミナのアニメで「PSYCHOPATH」という作品が人気を博しましたが、それを思い出します。警官の持つ拳銃型の「ドミネーター」という武器が、対象者の犯罪係数を算出し、その数値を元にその場で刑が執行されます。この犯罪係数が高いとその場で死刑執行も行われます。この世界はコンピュータ(最後にその中身も明かされる)に支配された人間の物語とも言えます。興味がありましたら見てみてください。映画化もされていた...かと思います(多分)。

 

科学が発達し、SF作品なども豊富な現代でこそ、想像のつくような世界を400年前の偉人が夢見ていたというのは想像もつきません。現代に生きる人々の中にも、私たちが理解出来ない様な未来をイメージして研究に取り組む天才もいるのかもしれないと思うと、私たちが夢物語だと思っていることも、500年も経てば現実化してるんでしょうか。最も最近の技術革新は10年前には映画の中の事だったことが当たり前のように行われるような進歩の仕方をしているので、それ以上となるとやはり想像できません(笑)。

 

なんだか、最終的に書評でもなんでもない所に帰着しましたが、興味深い本でしたので是非オススメしたいと思います。

 

また、著者の無限に関する本として別に「無限の果てになにがあるか」というものもありますのでそちらも読んでみてください。共に興味深い内容となっています。

日記的なの「統計の誤用を防ぐ書籍」

本日新宿の紀伊国屋に立ち寄りました。現在紀伊国屋では

 

「ダメな統計学を防ぐための書籍」

 

フェア?をやってるのか、理工書の階に行くと、中央のわかりやすい所にコーナーが設置されています。そこで見つけたフリーペーパーで「ダメな統計学を防ぐための書籍11冊」が紹介されていました(まあ、そのコーナーに置いてある本なんですが)どれも欲しい本ばかりなのですが、なにぶん学生故にお金が無いため買うことが出来ません...

 

あぁ、早く就職してお金が欲しい今日このごろです。

ちなみに、そのフリーペーパーの内容はブログに掲載されてる内容なので、次のリンクをクリックしてくれれば飛びますので(∩´。•ω•)⊃ドゾー

http://id.fnshr.info/2017/01/27/no-sdw/:ダメな統計学 11冊

「ダメな統計学~悲惨なほど完全なる手引書~」レビュー

先月あたりに本屋によった際に見つけ、ずっと気になっていた「ダメな統計学~悲惨なほど完全なる手引書」をついに購入してしまいました。

 

2017年2月時点でAmazon確率・統計カテゴリ1位を獲得したベストセラーです。

 ※上記画像はAmazonにリンクしています

 

内容としては、科学者が陥ってしまいがちな統計の誤用だったり、その原因、なぜ間違いなのかについて書かれています。

 

この本の原著というか、元となっているものは英語版、日本語版共にオンラインで見ることが出来るみたい。

 

英語の場合は「Statistics Done Wrong」

日本語の場合は「ダメな統計学

 

で検索すれば出てくるので検索してください。ってか、下にリンクつけときますのでクリックすれば飛ぶようになってます。

日本語↓

id.fnshr.info

英語↓

Welcome — Statistics Done Wrong

 評価 ★★★★★

アマゾンレビューするなら私は★5つです。

 

実践向けの統計学書があふれかえる世の中で、統計の誤用という視点から、誤用に文句をいうだけでなく、原因を解説し、防ぐ方法まで説明してくれる良書です。

 

統計学という分野は非常に応用範囲の広い分野であるのですが、曖昧さをおおいに含んだ分野でもあります。その為、使い方次第で誤解を生むことも多いです。最近はデータ分析がトレンドなためか、Excelで統計みたいな、実際に使う系の本が沢山出ています。また、科学の世界でも統計は必ず使われますが、統計をキチンと学んで使われていることは少ないと思われます。統計を使う場合、既に確立された手法を用いれば良く、また、Excelのようなソフトには統計の機能がついているため、ツールとしての側面が強い気がします。そのためか統計の曖昧な部分について知らずに使い、勘違いを生む原因になっている気がします。

 

この本では、科学の世界で日常的に使われている統計における、間違いやその原因について分かりやすく解説されています。

 

この本の目次は以下のようになっています。

  1. 統計的有意性入門
  2. 検定力と検定力の足りない統計
  3. 疑似反復:データを賢く選べ
  4. p値と基準率の誤り
  5. 有意性に関する間違った判断
  6. データの二度づけ
  7. 連続性の誤り
  8. モデルの乱用
  9. 研究者の自由:好ましい雰囲気?
  10. 誰もが間違える
  11. データを隠すこと
  12. 何ができるだろうか

ここで解説される原因は、手法的な面だけにとどまりません。研究者が統計を誤ってしまう環境的背景のような面についても解説されています。また、各章の最後には「ヒント」という項目が設けられており、そこではその章で解説されたダメな統計学を回避する方法を教えてくれています。

 

また、ページ下部の注釈が結構面白く、本文よりもそっちが気になってしまったりもします。

 

解説に出てくる論文などは、参考文献のページに網羅されていますのでもっと深く知りたいという人にも親切ですね。参考文献の数は数えてみたら191個ありました(笑)。

対象読者:統計学の経験がある人

数式がガリガリ出てくる本ではなく、ほぼほぼ文章ですので、読みにくいということは全く無いと思います。読み物としての本です。ただ、統計に全く触れたことのない人が読むには、統計用語が頻繁に出てくるため厳しいと思います。

統計学を普段使いしてる人や、統計学を学んでいる学生におすすめだと思います。

 

感想

私は大学で統計学を専門にすべく勉強しています。その中でいつも考えているのは「統計学が専門」というのはどういうことかということです。統計学はどの分野でも使われています。心理や経済学の学生なんかは統計を頻繁に使うと聞きます。物理なんかもそうですね。実際、日本では統計に関する専門書は、経済学者や、心理学者、物理学者が書いているということも多いです。みんな多かれ少なかれ統計を知っているのです。そしたら、統計を専門にするとはどういうことなのか?他の分野で統計をバリバリ使っている人との差は何所に生まれるのかというのをずっと思っていました。この本はその疑問に対する一つの答えを教えてくれた気がします。統計的手法の成り立ち、理論に精通し、それぞれの手法の留意点を知り、誤りを正せるレベルが求められているんだと私は思いました。一言に統計学と言っても幅広いので全てに精通することは難しいのかもしれませんが、自分が専門とする領域を見つけ、その部分だけでもまずは上記のレベルに達することができればと思います。

 

「統計でウソをつく法」は「ダメな統計学」ででてくる関連書籍の一つです。

ガウス積分

統計学を学んでいるとよく出てくるのがガウス積分です。ガウス積分は求め方を知っていれば対して解を得るのは対して難しくは無いですが、いちいち求めるには少し手間です。なので、積分結果を覚えておく。公式的に暗記してしまうと、話がスムーズに進みますので次の公式を導出していきたいとおもいます。A,aは実数定数として
$$
\int_{-\infty}^{+\infty}e^{-A(x-a)^2} dx = \sqrt{\frac{\pi}{A}}
$$

証明)
I = \int_{-\infty}^{+\infty} e^{-A(x-a)^2} dxと置くと、I^2
$$
I^2 = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} e^{-A(x-a)^2 -A(y-a)^2} dxdy = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} e^{-A\{(x-a)^2 + (y-a)^2\}} dxdy
$$
と書くことができる。
ここで
$$
\left\{
\begin{array}{ccc}
x-a & = & r \cos \theta \\
y-a & = & r \sin \theta
\end{array}
\right.
$$
と置く。ヤコビアン
$$
J = \left|
\begin{array}{cc}
\frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\
\frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta}
\end{array}
\right|
= \left|
\begin{array}{cc}
\cos \theta & -r \sin \theta \\
\sin \theta & r \cos \theta
\end{array}
\right| = r
$$
よって
\begin{eqnarray}
I^2 &=& \int_0^{2\pi} \int_0^{+\infty} e^{-A\{(r\cos \theta)^2 + (r\sin \theta)^2\}} |J| dr d\theta \\
&=& \int_0^{2\pi} \int_0^{+\infty} e^{-Ar^2} r dr d\theta \\
&=& \int_0^{2\pi} d\theta \int_0^{+\infty} re^{-Ar^2} dr \\
&=& [\theta]_0^{2\pi} [-\frac{1}{2A} e^{-Ar^2}]_0^{+\infty} \\
&=& 2\pi \cdot \frac{1}{2A} \\
&=& \frac{\pi}{A} \\
\therefore I &=& \sqrt{\frac{\pi}{A}} \;\;\;\;\; \because I > 0
\end{eqnarray}

参考文献

小寺平治(2014):『明快演習 数理統計』,共立出版

多変量解析~同時分布(Joint Distribution)~

久々にこのブログを書きます...前回書いたのはいつだったか...。最近はTexにまとめてるんでこっちのことを完全に忘れてました...

研究室に配属されて、多変量解析の勉強が本格的に始まってきました。まだ、ほとんどやっていないに等しいですが、気が向いた時に学んだことを覚え書きしていこうと思います。

自分が勉強に使っているのは研究室指定の
T.W.Anderson『An Introduction to Multivariate Statistical Analysis』

です。これの流れに沿って勉強をすすめていこうと思います。

積分布関数(1変量の場合)

多変量に入る前に1変量について見ていきます。1変量の累積分布関数(cumulative distribution function)は次のように定義される

 {\displaystyle
F(x) = P(X \le x) = P(\{ \omega \in \Omega : X(\omega) \le x \})
}

分布関数F(\cdot)は次の性質を持つ。

1. 任意のx \in \mathbb{R}^1 に対して 0 \le F(x) < 1でありかつ
{\displaystyle
F(-\infty) \equiv \lim_{x \to -\infty} F(x) = 0, \;\;\; F(+\infty) \equiv \lim_{x \to +\infty} F(x) = 1
}
2.F(x)は単調非減少である. \;\; : \;\; x < y \Leftrightarrow F(x) \le F(y)

3. F(x)は右側連続である.  \;\; : \;\; \lim_{y \to x+0} F(y) = F(x)


指数分布について密度関数と分布関数を見てみる。
指数分布の密度関数は
{\displaystyle
f(x) = \lambda e^{-\lambda x}
}
f:id:doratai:20161129220259j:plain:w300
であり、分布関数は
{\displaystyle
F(x) = 1-e^{-\lambda x}
}
で与えられる。
f:id:doratai:20161129220355j:plain:w300
密度関数は分布関数の微分で定義され、次の関係が成り立つ。
$$
f(x) = \frac{d}{dx}F(x) \Leftrightarrow F(x) = \int_{-\infty}^x f(u)du
$$

また、この他に密度関数は次の性質を満たす。

  • f(x) \ge 0
  •  \int_{-\infty}^{+\infty}f(x) dx = 1

2変量の場合

次は2変量の場合について考える。2つの確率変数X,Yを考える。c.d.f.がすべての実数の組x,yについて次で定義される。
$$
F(x,y) = Pr\{ X \le x, Y \le y\}
$$
ここで考えているF(x,y)は絶対連続の場合。つまり、ほとんど至る所(almost everywhere)で偏導関数が存在する場合を考える。つまり、ほとんど至る所で次が成り立つ。
\begin{eqnarray}
\frac{\partial^2 F(x,y)}{\partial x \partial y} &=& f(x,y) \\
F(x,y) &=& \int_{-\infty}^y \int_{-\infty}^x f(u,v)dudv
\end{eqnarray}
が成り立つものとして考える。ここで非負関数f(x,y) \ge 0XYの密度関数と呼ばれる。この関数は以下の性質を持つ。

  • f(x,y) \ge 0
  • \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(x,y) dxdy = 1

ほとんど至る所(almost everywhere)

ここで直接は関係無いが、少し気になるalmost everywhere について考えてみる。

定義

ほとんど至る所P(\omega)を命題関数とする.\{\omega \in \Omega : P(\omega) = 偽\}\mu-零集合である時、P(\omega)は"ほとんど至る所"で成立する.

ここで\mu-零集合とは、\mu(N) = 0 \;\;(\muは測度)なるN \in \mathcal{F} (\sigma -集合体)\mu-零可測集合といい、これが存在して、A \subset Nなる集合を\mu-零集合という。
ここで、\boldsymbol{\omega} \in \mathbb{R} \times \mathbb{R}とし、F(\cdot)をc.d.f.とする。命題関数を

$$
P(\boldsymbol{\omega}) = \left\{
\begin{array}{cc}
TRUE & if \;\; F(\boldsymbol{\omega})\;has\;a\;partial\;derivative\;at\;the\;point\;\boldsymbol{\omega} \\
FALSE & elsewhere
\end{array}
\right.
$$
で与える。先程述べた、偏導関数がほとんど至る所でん存在するとは、集合
$$
\{ \boldsymbol{\omega} \in \mathbb{R} \times \mathbb{R} :P(\boldsymbol{\omega}) = FALSE\}
$$
\mu-零集合であることを意味する。つまり
$$
\mu(\{ \boldsymbol{\omega} \in \mathbb{R} \times \mathbb{R} : P(\boldsymbol{\omega}) = FALSE\}) = 0
$$
これは、\mathbb{R} \times \mathbb{R}上の点\boldsymbol{\omega} = (x,y)について、偏導関数が存在しない点の集合の測度が0であることを意味する。

p変量の場合

今、p個の確率変数X_1,\cdots,X_pを考える。そのc.d.f.は
\begin{equation}
F(x_1,\ldots,x_p) = Pr\{X_1 \le x_1,\ldots,X_p \le x_p\}
\end{equation}
がすべての実数x_1,\ldots,x_pの集合によって定義される。密度関数はF(x_1,\ldots,x_p)が絶対連続であるならば
\begin{equation}
\frac{\partial^p F(x_1,\ldots,x_p)}{\partial x_1 \ldots \partial x_p} = f(x_1,\ldots,x_p)
\end{equation}
で与えられ、また
\begin{equation}
F(x_1,\ldots,x_p) = \int_{-\infty}^{x_p} \cdots \int_{-\infty}^{x_1} f(u_1,\ldots,u_p) du_1 \ldots du_p
\end{equation}
が成り立つ。p次元ユークリッド空間の任意の可測集合をRとする時、確率変数(X_1,\ldots,X_p)Rに属する確率は
\begin{equation}
Pr\{(X_1,\ldots,X_p) \in R\} = \underset{R}{\idotsint} f(x_1,\ldots,x_p) dx_1\ldots dx_p
\end{equation}
確率要素f(x_1,\cdots,x_p)\Delta x_1 \cdots \Delta x_pはほぼ確率P(x_1 \le X_1 \le x_1 + \Delta x_1,\cdots, x_p \le X_p \le x_p + \Delta x_p)
に等しい。
もしf(x_1,\cdots,x_p)が連続であるならば、同時積率は次で定義される。
\begin{equation}
E(X_1^{h_1}\cdots X_p^{h_p}) = \int_{-\infty}^{+\infty} \cdots \int_{-\infty}^{+\infty} x_1^{h_1} \cdots x_p^{h_p} f(x_1,\cdots,x_p) dx_1\cdots dx_p
\end{equation}

参考図書

T.W. Anderson(2003):『An introduction to Multivariate Statistical Analysis』, John Wiley & Sons
梅垣寿春,塚田真,大矢雅則(2015) : 『測度・積分・確率』,共立出版



書いてから気づいたんですが、似たような記事を前回も書いてあるみたいです...今回のが少し内容が重くなってるんでまあいいかなと...

数理統計学の勝利~ニューヨークタイムズのネイト・シルバーの数理モデル予測が全50州で的中~(外部記事)

統計学がアメリカで政治学者相手に大勝したようですね。
政治学者はどのような思考回路で政治予測や分析をしてるかはわかりませんが、計算機によって膨大なデータを処理することになったこの時代、1人の人間が持つ経験や思考では、もはや上回ることはできないでしょう。


ネイト・シルバー - Wikipedia


以下引用

New York Timesの選挙予測専門家、ネイト・シルバーは昨夜、大統領選の勝敗を全50州で的中させた。 その一方で、いわゆる政治専門家たちの予想はほとんどが外れた。中には笑うしかないような外れ方をした者もいる。

ネイト・シルバーについてはテレビのゲストに呼ばれる政治専門家が口を揃えて「リベラルに偏った見解」と非難してきた。しかしシルバーは今回も彼の作った数理的予測モデルが古臭い専門家の勘や生半可な統計に基づく推測より圧倒的に優れていたことを証明した。

残る疑問は、数理モデルのこれほどの有効性を見た後でもテレビのプロデューサーたちは時代遅れの政治専門家なるものを番組に使い続けるつもりなのかどうかという点だけだ。

シルバーの数理モデルの特長は、どんな政治専門家もとうてい考慮しきれないほど膨大な量の数値を入力として用いるところにある。シルバー・モデルでは各種の世論調査の結果を、規模、質、時期などによって重み付けし、過去の同種の選挙結果と照合される(もちろんそれ以外にもさまざまな高度な統計処理が用いられている)。

今回、予測を100票も外した〔大統領選挙人の総数は538人〕専門家はこんなことを言っていた。曰く、オバマ大統領にはもはや伝えるべきメッセージがない、問題意識がない、 過去4年間の業績に対する説明責任を果たしていない、それを有権者は見ぬいている…。そう書いたのはクリントン大統領の元補佐官、ディック・モリスだが、彼の予測と現実はグランドキャニオンくらいかけ離れていた(ロムニーが325票獲得するというモリスの予測は100票以上外れていた)。

シルバーのアプローチの成功はテレビ局にジレンマを与えている。第一に、この種の議論を理解するためには視聴者に数学の素養が必要だ。仮に古臭い政治評論家を数理統計学者で置き換えたとしても、今度は番組同士で自分たちの予測の優位性を説明するためには面倒な統計学の議論が必要になる。視聴者はそんな議論にはすぐに飽きてしまうだろう。

第2に、ショッキングな選挙予測を報じて視聴率を取りに行けなくなる。視聴率を稼げる意外な結果は、ほとんどの場合不正確な予測だ。ところがシルバー・モデルは多数の世論調査を詳しく分析して平均を出しているので概ね常識的で安定した(番組としては退屈な)結果が出る。

しかしシルバー・モデルが与えるもっとも大きく、破壊的な影響は、伝統的な選挙キャンペーンや政治評論はもはや選挙結果に決定的な影響を与えることはないという事実が明らかになってしまうことだ。シルバー・モデルは選挙の数ヶ月前からオバマの勝利がほぼ確実であると予測していた。選挙は現職有利というのがセオリーであり、ロムニーにはそれを覆すだけのカリスマが欠けており、共和党内でさえそれは意識されていた。また他の重要な要素、景気や失業率に選挙運動は何の影響も与えることができない。選挙を前に景気が上向けば保守系の挑戦者は苦戦を免れない。

つまり、「アメリカ人はもはやオバマのリベラルな社会政策を見放した」云々という「専門家」たちの御託宣はまったく現実とは関係がなかったわけだ。テレビに毎日現れる政治評論家、選挙専門家の発言は大部分がたわごとだった。しかしテレビのプロデューサーたちは派手な党派的な議論、不正確だが一般受けしやすい選挙予測などによって視聴率を稼ごうとする強い動機を持っている。シルバーは今回、ひとつの戦闘には勝ったものの、戦争に勝つのはまだ先のことになりそうだ。

〔日本版〕ネイト・シルバーはアメリカでもっとも注目されている選挙専門家。シカゴ大学経済学部を卒業した後、2002年に、KPMG会計事務所に勤務中にメジャー・リーグ野球選手の統計的評価システムPECOTAをを開発した。これは映画「マネー・ボール」で日本でも知られるようになったセイバーメトリクスをオンライン化したもので、2007年にはPECOTAをBaseballProspectus社に売却して、選挙予測の分野に進出。2008年のブッシュ対ゴアオバマ対マケインの大統領選で50州中49州の勝敗を的中させ、いちやく注目されるようになった。現在、FiveThirtyEightブログはニューヨークタイムズの一部となっている。著書にThe Signal and the Noiseなどがある。)

元記事↓
大統領選でニューヨークタイムズのネイト・シルバーの数理モデル予測が全50州で的中―政治専門家はもはや不要? | TechCrunch Japan


統計学は理系だけでなく文系分野と考えられるところでも非常によく使われています。実務の面での汎用性は現状、"最強の学問"です


多変量解析~多変量正規分布の標準化~

1変量の時の標準化はそんなに苦では無いですよね?ここではp変量の多変量正規分布の標準化をやっていきたいと思います。

まずは多変量正規分布の確認

\boldsymbol{X} \sim N(\boldsymbol{\mu},\Sigma) とする。\Sigma は正定値行列より \Sigma = CC^{\mathrm{T}} なる正則行列 C が存在する。

多変量正規分布


f(\boldsymbol{x}) = \frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}} \exp \{-\frac{1}{2}(\boldsymbol{x-\boldsymbol{\mu}})^{\mathrm{T}}\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}) \}

標準化

\boldsymbol{Y} = C^{-1}(\boldsymbol{X}-\boldsymbol{\mu}) とする。これの逆変換が \boldsymbol{x} = C\boldsymbol{y} + \boldsymbol\mu で与えられる。

ヤコビアンを求めておきます。

$$
J(y_1,\cdots,y_p) = \mod \left|
\begin{array}{ccc}
\frac{\partial x_1}{\partial y_1} & \cdots & \frac{\partial x_1}{\partial y_p} \\
\vdots & \ddots & \vdots \\
\frac{\partial x_p}{\partial y_1} & \cdots & \frac{\partial x_p}{\partial y_p}
\end{array}
\right|
= \mod \left|
\begin{array}{ccc}
c_{11} & \cdots & c_{1p} \\
\vdots & \ddots & \vdots \\
c_{p1} & \cdots & c_{pp}
\end{array}
\right| = \mod |C|
$$

また、\Sigma = CC^{\mathrm{T}} より


C^{-1}\Sigma(C^{\mathrm{T}})^{-1} = I


この時、

|C^{-1}| | \Sigma | | (C^{\mathrm{T}})^{-1} = |I|
\frac{1}{|C|} |\Sigma| \frac{1}{|C|} = 1
|\Sigma| = |C|^2


である。これらから



g(\boldsymbol{y})dy_1\cdots dy_p = f(\boldsymbol{x})dx_1 \cdots dx_p
 = \frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}\exp \{-\frac{1}{2}(C\boldsymbol{y} + \boldsymbol{\mu} - \boldsymbol{\mu})^{\mathrm{T}} \Sigma^{-1} (C\boldsymbol{y} + \boldsymbol{\mu} - \boldsymbol{\mu}) \} J(y_1,\cdots,y_p) dy_1\cdots dy_p
= \frac{1}{(2\pi)^{\frac{p}{2}}|C|} \exp \{-\frac{1}{2} (C\boldsymbol{y})^{\mathrm{T}} \} (CC^{\mathrm{T}})^{-1} (C\boldsymbol{y}) \mod |C| dy_1\cdots dy_p
= \frac{1}{(2\pi)^{\frac{p}{2}}} \exp \{-\frac{1}{2}\boldsymbol{y}^{\mathrm{T}}C^{\mathrm{T}} (C^{\mathrm{T}})^{-1}
 C^{-1}C \boldsymbol{y} \}dy_1\cdots dy_p
= \frac{1}{(2\pi)^{\frac{p}{2}}} \exp \{-\frac{1}{2} \boldsymbol{y}^{\mathrm{T}} \boldsymbol{y} \} dy_1\cdots dy_p
\therefore g(\boldsymbol{y}) = \frac{1}{(2\pi)^{\frac{p}{2}}} \exp \{-\frac{1}{2} \boldsymbol{y}^{\mathrm{T}} \boldsymbol{y} \}



以上から Y はp変量標準正規分布正規分布 N(\boldsymbol{0},I) に従うので変数変換 Y = C^{-1}(\boldsymbol{X} - \boldsymbol{\mu}) は標準化である。

参考文献

特に無し

ベクトル微分


多変量解析を勉強するにあたって、必要になることがあるのがベクトルの微分である。

まずはまとめから


$$
\begin{eqnarray}
\frac{\partial(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta})}{\partial\boldsymbol{\beta}} &=& \boldsymbol{C} \tag{1} \\
\frac{\partial(\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta})}{\partial\boldsymbol{\beta}} &=& (A + A^{\mathrm{T}})\boldsymbol{\beta} \tag{2}
\end{eqnarray}
$$

これらを証明していく。ベクトルの微分を考えていくうえでは、面倒だが、成分を
考えていくことになる。

(1)の証明

$$
\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta} = c_1\beta_1 + c_2\beta_2 + \cdots + c_p\beta_p
$$

より

$$
\frac{\partial(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta})}{\partial\beta_i} = c_i
$$

ただし i = 1,2,\cdots,p

よって

$$
\frac{\partial(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta})}{\partial\boldsymbol{\beta}} = \left(
\begin{array}{c}
\frac{\partial(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta})}{\partial\beta_1} \\
\frac{\partial(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta})}{\partial\beta_2} \\
\vdots \\
\frac{\partial(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{\beta})}{\partial\beta_p}
\end{array}
\right)
= \left(
\begin{array}{c}
c_1 \\
c_2 \\
\vdots \\
c_p
\end{array}
\right)
=\boldsymbol{C}
$$

(2)の証明

$$
\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta} = \sum_{i=1}^p\sum_{j=1}^p a_{ij}\beta_i\beta_j
$$

より

$$
\begin{eqnarray}
\frac{\partial(\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta})}{\partial\beta_k} &=& \sum_{j=1}^p a_{kj}\beta_j + \sum_{i=1}^p a_{ik}\beta_i \\
&=& (a_{k1},\cdots,a_{kp})\boldsymbol{\beta} + (a_{1k},\cdots,a_{pk})\boldsymbol{\beta} \\
&=& \{(a_{k1},\cdots,a_{kp}) + (a_{1k},\cdots,a_{pk})\}\boldsymbol{\beta} \\
&=& (\boldsymbol{a}_k + \boldsymbol{a}_k^{\mathrm{T}})\boldsymbol{\beta}
\end{eqnarray}
$$

ここで

$$
A = \left(
\begin{array}{cccc}
a_{11} & a_{12} & \ldots & a_{1p} \\
a_{21} & a_{22} & \ldots & a_{2p} \\
\vdots & \vdots & \ddots & \vdots \\
a_{p1} & a_{p2} & \ldots & a_{pp}
\end{array}
\right)
\\
A^{\mathrm{T}} = \left(
\begin{array}{cccc}
a_{11} & a_{21} & \ldots & a_{p1} \\
a_{12} & a_{22} & \ldots & a_{p2} \\
\vdots & \vdots & \ddots & \vdots \\
a_{1p} & a_{2p} & \ldots & a_{pp}
\end{array}
\right)
$$

であるから、

$$
\frac{\partial\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta}}{\partial\boldsymbol{\beta}}
= \left(
\begin{array}{c}
\frac{\partial(\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta})}{\partial\beta_1} \\
\frac{\partial(\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta})}{\partial\beta_2} \\
\vdots \\
\frac{\partial(\boldsymbol{\beta}^{\mathrm{T}}A\boldsymbol{\beta})}{\partial\beta_p}
\end{array}
\right)
= (A + A^{\mathrm{T}})\boldsymbol{\beta}
$$

となる。

参考文献

特に無し