2017-09-20

二項分布の3次のモーメント

統計学の教科書において確率分布を紹介するときに2次のモーメントまでしか載っていないことが多い。

モーメント母関数を用いれば、何次のモーメントでも求めることができるので、何の役に立つかは分からないが、ここでは二項分布の3次のモーメントを求めてみよう。

n次のモーメントの求め方は、モーメント母関数を $n$ 回微分して、 $\xi$ を0と置くのであった。まず、二項分布のモーメント母関数を2階微分をしたものは、

　　 $f^{\prime\prime}(\xi) = np \{ (n-1) pe^{2\xi}(pe^\xi + q)^{n-2} + e^\xi(pe^\xi + q)^{n-1} \}$

である。これをさらに微分すると、

　　 $\begin{split} f^{\prime\prime\prime}(\xi) &= np \{ 2pe^{2\xi}(n-1) (pe^\xi + q)^{n-2} + p^2e^{3\xi}(n-1)(n-2)(pe^\xi +q)^{n-3}\\ \\ &\quad + e^\xi (pe^\xi + q)^{n-1} + pe^{2\xi} (n-1)(pe^\xi + q)^{n-2} \} \end{split}$

となる。そこで、 $\xi =0$ とおき、 $p$ の次数で整理すると、

　　 $\begin{split} f^{\prime\prime\prime}(0) &= np \{ 2p(n−1) + p^2(n−1)(n−2)+ 1 + p(n−1) \} \\ \\ &= 2np^2(n−1)+np^3(n−1)(n−2)+np+np^2 (n−1)\\ \\ &= n(n−1)(n−2)p^3 + 3n(n−1)p^2 + np \end{split}$

となり原点まわりの3次のモーメントが得られる。これを利用し、今度は平均値周りの3次のモーメントを求めてみる。平均値周りのモーメントは、

　　 $m_3 = E((x - \mu)^3)$

であり、これを展開すると、

　　 $m_3 = E(x^3) -3\mu E(x^2) + 2\mu^3$

となる。 $\mu = E(x) = np$ であるので、先に求めた原点まわりのモーメントをこの式に代入することで、 $m_3$ を求めることができる。ちなみに、2次の原点周りのモーメントは、

　　 $E(x^2) = n(n-1)p^2 + np$

である。

まず、それぞれを代入する。

　　 $\begin{split} m_3 &= n(n−1)(n−2)p^3 + 3n(n−1)p^2 + np -3np \{ n(n-1)p^2 + np \} + 2n^3p^3\\ \\ &= np \{ (n−1)(n−2)p^2 + 3(n−1)p + 1 - 3(n-1)np^2 -3np + 2n^2p^2 \} \end{split}$

ごちゃごちゃしていて分かりづらいが、{ }の中を $p$ の次数にそって整理し、各係数を計算していくと、

　　 $\begin{split} m_3 &= np \{ (n−1)(n−2)p^2 + 3(n−1)p + 1 - 3(n-1)np^2 -3np + 2n^2p^2 \} \\ \\ &= np \{ \{(n−1)(n−2)- 3(n-1)n + 2n^2 \}p^2 + \{3(n−1) - 3n \} p + 1 \} \\ \\ &= np (2p^2 - 3p +1)\\ \\ &= np (1 - p)(1 - 2p) \end{split}$

が求まる。

3次のモーメントを標準偏差の3乗 $\sigma^3$ で割ったものは歪度と呼ばれて分布の左右対称性の指標とされたりする。 $q = (1-p)$ とすると分散は $npq$ 、 $\sigma = \sqrt{npq}$ と表せるので、

　　 $\displaystyle 歪度 = \frac{m_3}{\sigma^3} = \frac{npq(1 - 2p)}{(\sqrt{npq})^3} = \frac{npq(1 - 2p)}{npq \sqrt{npq}} = \frac{1 - 2p}{ \sqrt{npq}}$

となる。分子に注目すると、 $p =0.5$ のときに歪度は0となることがわかる。 $p > 0.5$ では、歪度は負の値をとり、 $p < 0.5$ では歪度は正の値をとる。歪度が0のときに分布は左右対称で、歪度が正のときには右の裾が長く、負のときには左の裾が長い分布が得られる。実際に二項分布の確率分布を確認してみてもその通りになっていることが分かる。

2017-09-19

モーメントについて

調べたこと・考えたこと

モーメントについての学習メモ。

モーメント（積率）という概念がある。これは、確率分布の性質のうちのいくつかを数値で表したものである。平均値や分散というのもモーメントの一種である。

離散値をとる確率変数 $x$ について、その確率分布を $P(x)$ としたときに、 $x^n$ の期待値を $x$ の確率分布の $n$ 次のモーメントという。

　　　 $\acute{m}_n = E(x^n) = \sum_x x^n P(x)$

これは、原点まわりのモーメントと呼ばれたりする。 $n=1$ である1次のモーメントはそのまま定義通りに平均値である。

一般にモーメントが用いられるのは平均値のまわりらしく、確率変数 $x$ と $\mu$ の差のべき指数で求められるモーメントである。

　　　 $m_n = E \big( (x - \mu)^n \big) = \sum_x (x - \mu)^n P (x)$

こういったものを平均値のまわりの $n$ 次のモーメントと呼ぶ。原点まわりと平均値まわりをここでは $\prime$ の有無で区別する。平均値まわりの1次のモーメントを求めると以下のようになる。（ここで、 $\sum_x P(x)$ は確率の定義から1である。※サイコロのすべての目の確率は足すと1。）

　　　 $m_1 = \sum_x (x -\mu)P(x) = \sum_x xP(x) -\mu \sum_x P(x) = \mu - \mu = 0$

平均値まわりの２次のモーメント（つまりは分散）を計算するには、定義式のものより原点まわりのモーメントを用いて次のように計算すると簡単らしい。

　　　 $m_2 = \acute{m}_2 - \mu^2 = E(x^2) - \mu^2$

これは何かというと、分散の簡便な計算式とされる「2乗の平均 - 平均の2乗」にほかならない。

さて、モーメントを求める際には、モーメントを生み出すモーメント母関数というものを用いると便利であるらしい。モーメント母関数とは $\xi$ (グザイとかクシーと読む)についての次のような関数である。

　　　 $f(\xi) = E (e^{\xi x})$

指数関数の級数展開の公式というものがあり、

　　　 $\displaystyle e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dots$

これを用いると、モーメント母関数の式は、

　　　 $\displaystyle f(\xi) = E \big( 1 + \xi x + \frac{\xi^2 x^2}{2!} + \frac{\xi^3 x^3}{3!} + \dots \big)$

このようになる。さらに式を変形すると、

　　　 $\displaystyle f(\xi) = 1 + \xi E(x) + \frac{\xi^2}{2!} E(x^2)+ \frac{\xi^3}{3!} E(x^3) + \dots$

このようになり、これを最初に紹介した原点まわりのモーメントで置き換えると、

　　　 $\displaystyle f(\xi) = 1 + \xi \acute{m}_1 + \frac{\xi^2}{2!} \acute{m}_2+ \frac{\xi^3}{3!} \acute{m}_3 + \dots$

が得られる。展開式の係数に次々と高次のモーメントが現れてくる。この母関数を用いると、面倒な計算を避けてお目当てのモーメントを計算することができるというのだ。

$n$ 次の原点まわりのモーメント $\acute{m}_n$ を求めたいのであれば、 $f(\xi)$ を $\xi$ について $n$ 回微分をし、 $\xi = 0$ とすれば求められる。先に求めた展開式を見れば分かるように、 $\xi = 0$ とすることで、 $n$ よりも右の項はすべて0になるので目当てとするモーメントだけが取り出せるわけである。

二項分布を例にして考えてみよう。二項分布の確率関数は次の式で表される。

　　　 $P(x) = {}_n C _x p^x q^{n-x}$

まずは、これをモーメント母関数に埋め込み、多少の変形を行う。

　　　 $\begin{equation*} \begin{split} f(\xi) &= \sum_{x=0} ^n e^{\xi x} P(x) \\ &= \sum_{x=0} ^n e^{\xi x} {}_n C _x p^x q^{n-x}\\ &= \sum_{x=0} ^n {}_n C _x (pe^\xi)^x q^{n-x} \end{split} \end{equation*}$

最右辺の式は、二項定理を用いると次のように書くことができる。

　　　 $f(\xi) = (q + pe^\xi)^n$

この式を用いて、 $\acute{m}_1$ を求める。1次のモーメントなので、1回微分を行い、 $\xi = 0$ とおく。（ $(p+q)$ は1、および、eの0乗は1なことから、下から２番目の式の（）内は1である）

　　 $\begin{equation*} \begin{split} \acute{m}_1 &= \frac{df(\xi)}{d\xi}\\ &= npe^\xi(q + pe^\xi)^{n-1} \\ f^{\prime}(0) &= np \end{split} \end{equation*}$

二項分布の平均値である、 $np$ が無事に算出された。続いて、 $\acute{m}_2$ を求めてみよう。まず、2階微分を行う。

$\begin{equation*} \begin{split} \acute{m}_2 &= \frac{d^2f(\xi)}{d\xi^2}\\ &= npe^\xi \big((n-1)(q + pe^\xi)^{n-2} pe^\xi\big) + npe^\xi (q + pe^\xi)^{n-1}\\ &= npe^\xi \big( (n-1)(q + pe^\xi)^{n-2} pe^\xi + (q + pe^\xi)^{n-1} \big)\\ \end{split} \end{equation*}$

そして、 $\xi = 0$ とすると、

$\begin{equation*} \begin{split} f^{\prime\prime}(0) &= np \big( (n-1) p + 1 \big)\\ &= np (np + 1 - p ) \\ &= np (np + q ) \\ &= n^2p^2 + npq\\ \end{split} \end{equation*}$

$\acute{m}_2 = n^2p^2 + npq$ が得られる。あとは、最初の方に述べたように、原点まわりのモーメントから平均の2乗、すなわち $n^2p^2$ を引くことで、

　　　 $m_2 = npq$

が得られる。二項分布の分散も無事に算出された。

このように、モーメント母関数を用いれば様々な確率分布のモーメントを求めることができる。

【参考】
『キーポイント確率・統計』（岩波書店）

2017-09-18

続続・自由度について（実験修正編）

調べたこと・考えたこと

以前書いた記事で不偏分散と標本の分散の比較し、そのときは別々の標本からそれぞれの統計量を計算した。

nekomosyakushimo.hatenablog.com

しかし、よくよく考えると同じ標本からそれぞれの統計量を計算した方が比較の観点からは正しかったのだと気づいた（標本誤差が減るので）。寝ぼけていたのかもしれない。

つまり、

#不偏分散と標本分散の比較
for (i in 1:10000){
  dat <- rnorm(10,100,15)
  dat_mean <- mean(dat)
  sample_var[i] <- sum((dat - dat_mean)^2)/10
  unbias_var[i] <- var(dat)
}

とすることで、純粋に分散の計算方法だけの違いが見れる。期待値を計算して修正式に当てはめると、

mean(unbias_var)
 # [1] 226.4204
mean(sample_var) *10/9
 # [1] 226.4204

このようにぴったりと一致する。

2017-09-16

続・自由度について（別証明編）

調べたこと・考えたこと

前の記事までで、標本のデータから母分散を推定するには、n-1で計算した不偏分散を用いると不偏推定量が得られることについて書いてきた。

nekomosyakushimo.hatenablog.com
nekomosyakushimo.hatenablog.com

今回はそのことについての別の観点からの学習メモ。

標準正規分布 $N (0, 1 )$ からの $n$ 個の独立な標本値の平方和は自由度 $n$ の $\chi^2$ 分布に従うことが知られている。

　　 $\displaystyle \chi^2 = \sum^n_i \big( \frac{x_i - \mu}{\sigma} \big) ^2$

このとき、 $\mu$ ではなく $\bar{x}$ を用いたものは、自由度が1減り n - 1 の $\chi^2$ 分布に従う。

　　 $\displaystyle \chi^2 = \sum^n_i \frac{(x_i - \bar{x})^2}{\sigma^2}$

この式は

　　 $\displaystyle \skew{0}\hat{\sigma}^2 = \frac{\sum^n_i (x_i - \bar{x})^2}{n -1 }$

を用いると次のように変形できる。まず、分子分母に $n - 1$ をかける。

　　 $\displaystyle \chi^2 = \sum^n_i \frac{(x_i - \bar{x})^2 (n - 1)}{(n - 1) \sigma^2}$

次に、 $\frac{\sum^n_i (x_i - \bar{x})^2}{n -1 }$ を $\skew{0}\hat{\sigma}^2$ に置き換える。

　　 $\displaystyle \chi^2 = \frac{(n - 1) \skew{0}\hat{\sigma}^2}{\sigma^2}$

ところで、自由度 $n$ の $\chi^2$ 分布の平均値はその自由度 $n$ に等しいことが分かっている。前の式の $\chi＾2$ の自由度は n-1であるので、次の式で表せる。

　　 $\displaystyle E(\chi^2) = E \big( \frac{(n - 1)\skew{0}\hat{\sigma}^2}{\sigma^2} \big) = \frac{(n - 1)}{\sigma^2} E (\skew{0}\hat{\sigma}^2) = n - 1$

ここで最右辺とその1つ左の式を変形すると次の式が導ける。

　　 $\displaystyle E(\skew{0}\hat{\sigma}^2 ) = \sigma^2$

つまり、不偏分散の標本分布の期待値は、母分散に一致するということである。これは、不偏分散が母分散の不偏推定量であることの $\chi^2$ 分布を用いた証明である。

【参考】
『基本統計学』（有斐閣）

2017-09-15

続・自由度について（実験編）

調べたこと・考えたこと

この記事は不正確な内容を含みます。下記の新しい方の記事が正確です。（2017/9/18）
nekomosyakushimo.hatenablog.com

ーーーーー

前の記事で、標本から母分散を推定するときに、nではなく自由度であるn-1を用いることについて書いた。今回はそのことをRを用いた簡単な実験で示す。

母平均が100、母分散が225（母標準偏差は15）の正規分布からn=10の標本を抽出して、その分散を求める。この手続きを10000回繰り返して、経験的に標準偏差の標本分布を作ることにする。これを、 n で割って普通に分散を求めた場合と、n-1で割って不偏分散を求めた場合とを比較するのが今回の目的である。

スクリプトは難しいところは何もないが、Rのvar( )関数は不偏分散を求める仕様であるので、普通の分散を求める方では、平均と標本サイズから定義通りの計算をしている。

#標本分散と不偏分散の実験
sample_var <- numeric(10000)
unbias_var <- numeric(10000)

#普通の分散
for (i in 1:10000){
  dat <- rnorm(10,100,15)
  dat_mean <- mean(dat)
  sample_var[i] <- sum((dat - dat_mean)^2)/10
}

#不偏分散
for (j in 1:10000){
  dat <- rnorm(10,100,15)
  unbias_var[j] <- var(dat)
}

得られた結果をヒストグラムにして重ねたのが次のもの*1。

f:id:nekomosyakushimo:20170913235955p:plain

ピンクの方が nで割った普通の分散で、青い方が n-1 で割った不偏分散である。当然のことながら不偏分散の方が平均的に大きくな値を取っている。

それぞれの平均（つまりこの経験的な標本分布における期待値）を求めると以下の通りである。

> mean(sample_var)
[1] 200.6958
> mean(unbias_var)
[1] 223.8554

不偏分散の方は、母分散である225に近い値だが、普通の分散の方はやはり小さく推定してしまっている。普通の分散の方に、n / n-1 をかけると不偏分散による推定量に近づく。

> mean(sample_var) * 10/9
[1] 222.9954

以上、簡易的な実験により経験的に求めた分散の標本分布からも、不偏分散が不偏推定量になっていることが確認できる。

*1:色付きヒストグラムを重ねるときにはヒストグラム | を参照した。

2017-09-14

自由度について

調べたこと・考えたこと

統計関係の概念の中でも理解しづらいものの一つに「自由度」がある。得られたデータから、母分散の不偏推定量を求める時に、nではなくて、n -1　で割るだとか説明されて分かったような分かっていないような気になるアレである。それに関連する学習のメモ。

まず、標本平均の標本分布の話から。ある母集団から標本サイズがnである標本を抽出した際に、その標本平均は様々な値をとる。例えば、母平均が100、母分散が225（母集団標準偏差は15）のIQをとる集団があったとして、その集団から適当に10人ぐらいを抽出して標本平均を求めると、当然のことながら母平均と同じ値になることはあまりない（例えば、102かもしれないし90かもしれない）。ただ、母平均から離れた値はとりづらい（例えば標本平均が60なんかになることは稀である）。このようにして、「標本平均がどのような値を取りやすいか」についての分布を標本平均の標本分布と言う。

この標本平均の標本分布は、平均が母集団と同じμ、分散がσ/nを持つ正規分布に従うことが知られている。

$\displaystyle \mu_\overline{x} = \mu$

$\displaystyle \sigma_\overline{x}^2 = \frac{\sigma^2}{n}$

IQを例で言えば、この標本平均のとる値は、平均100、分散22.5の正規分布に従い、確率的に変動するということである。

ここまでを前提に自由度の話に入る。標本平均の期待値は、母集団の平均と一致している。こういうものを不偏推定量と言う。では、標本の分散の期待値はどうだろうか。結論から言えば、母分散よりも小さい値になることが分かっている。つまり、標本の分散を母分散の推定量として扱うと、母分散を小さく見積もることとなってしまう。

なぜそのようなことが起きるかというと、標本の分散を算出する際に標本平均を使っていることに由来している。標本平均は先に確認したようにそれ自体が標本抽出のたびに変わる確率的なものであり、μを中心に、σ/nの分散で分布している。その関係もあって、母分散は、標本平均を中心とする分布の分散（標本の分散）と標本平均の標本分布における分散の和と等しくなる。このことを表したのが次式である。

$\displaystyle \sigma^2 = \frac{\sum^n_i (x_i - \overline{x})^2}{n} + \frac{\sigma^2}{n}$

右辺の第１項が標本の分散の部分、第２項が標本平均の標本分布の分散である。この式を次のように変形していく。まず、第２項を左辺に移行して通分する。

$\displaystyle \frac{n\sigma^2 - \sigma^2}{n} = \frac{\sum^n_i (x_i - \overline{x})^2}{n}$

左辺のσをくくりだす。

$\displaystyle \frac{\sigma^2(n -1)}{n} = \frac{\sum^n_i (x_i - \overline{x})^2}{n}$

両辺をnで割って（n -1 ）を右辺へ移行して完成。

$\displaystyle \sigma^2 = \frac{\sum^n_i (x_i - \overline{x})^2}{n - 1}$

というわけで、件の n -1 たる数字が導出された訳である。最後の式の右辺で得られた推定量は不偏推定量であり、不偏分散と呼ばれることもある。n - 1というのは、標本の分散から母分散を推定するにあたって、標本平均それ自体のバラツキを勘案し、過小推定を調整したものだと言える。

【参考】
「心理学のためのデータ解析テクニカルブック」（北大路書房）

2017-09-13

連続のための修正とは何か

調べたこと・考えたこと

連続のための修正（連続性の補正）について調べたことのメモ。（まだ理解は不完全なのであまり信用できる記事ではありません。）

やや遠回りではあるが、まずは離散変量での検定について考える。

コインを10回投げて表の出る回数（確率変数）というのは、n=10, p=0.5 の二項分布に従う。それぞれの変数のとる確率はdbinom関数で簡単に求められる。以下はその表および確率分布の図である。（z値はあとで使うので今は無視してよい）

表の回数	確率	z値
0	0.001	-3.162
1	0.010	-2.530
2	0.044	-1.897
3	0.117	-1.265
4	0.205	-0.632
5	0.246	0.000
6	0.205	0.632
7	0.117	1.265
8	0.044	1.897
9	0.010	2.530
10	0.001	3.162

f:id:nekomosyakushimo:20170912222911p:plain

さて、ここでこのコインが本当に0.5の確率で表がでるものなのか（細工がされていないか）調べたいとする。帰無仮説をp=0.5として、例えば、表の回数が（0、1）回あるいは（9、10）回だったらこの帰無仮説を棄却するものとしよう（極端に裏ばかりでる、あるいは極端に表ばかりでるようだったらp=0.5ではないとする両側検定）。通常だと有意水準を定めて、その水準を超えるようだったら帰無仮説を棄却とするだろうが、とりあえずここでは、上記の回数だったら棄却すると設定する。このとき、第１種の誤りαの確率は、それぞれの確率変数のとる確率の和である。すなわち、

coin_prob <- dbinom(0:10,size=10,prob=0.5)
sum(coin_prob[1:2]) +sum(coin_prob[10:11])
[1] 0.02148438

であり、0.021%程度の確率で帰無仮説が真なのに、それを誤って棄却してしまうことになる。

ところで、二項分布はnが大きくなれば正規分布に近づくことが知られている。この性質を利用して、得られた値について検定を行おうという考えがでてくる。実現値を x 、平均を μ 、標準偏差を σ とすると、次のように z へと標準化ができる。

$\displaystyle z = \frac{x - \mu_x}{\sigma_x} = \frac{x - np}{\sqrt{ np (1-p)}}$

この変換をして得られたzの値は、最初に載せた表に記載してある。この z は近似的に、平均0、分散1の標準正規分布に従う。

さて、このzの値を用いて有意確率を計算してみる。表の回数が（0、1）回あるいは（9、10）回というのは, z値が |z| > 2.53 になる確率なので、正規分布が平均に対称なのを利用し、

pnorm(-2.53) * 2
[1] 0.01140625

と求められる。αの確率は、0.011%となり、最初に二項分布を用いて求めたときよりも、第１種の過誤を犯す確率がかなり小さく評価されてしまうこととなった。

これは、もともと二項分布が離散変量に基づく分布にもかかわらず、xを連続変量として扱っていることに原因があるとされている。離散的な整数を連続変量として当てはめているわけだが、それらの離散的な整数と整数の間にも確率があるわけで、この確率を考慮しないと、連続的変量としての分布への当てはまりが悪くなる。

※下の方にその解説がある。
統計学入門−第3章

離散変量をもとに連続変量の分布に近似して検定を行う際には、次に示す連続のための修正（イェーツの修正や半整数補正とも呼ばれる）を施した方が正規分布への近似がよくなるとされる。

$\displaystyle z = \frac{ | x - np | - 0.5 }{\sqrt{ np (1-p)}}$

今の例で考えると、1.5、8.5のところからz値を求めると正規分布により正確に近似されるはずである。上の式をもとにして z を修正し、その値を用いてαの確率を求めると以下のようになる。

(abs(1 - 10*0.5) - 0.5) / sqrt(10*0.5*0.5)
[1] 2.213594
(1 - pnorm(2.213594)) *2
[1] 0.02685672

αの確率は0.027%で、何も考えずにz値に変換した場合に比べて、二項分布で求めていたときの有意水準αに近づいており、第１種の過誤の過小推定の危険性を避けている。（今度は大きく評価するようになったが。）

とりあえずここまで。

【参考】
『心理学のためのデータ解析テクニカルブック』（北大路書房）

心理学のためのデータ解析テクニカルブック

作者: 森敏昭,吉田寿夫
出版社/メーカー: 北大路書房
発売日: 1990/06/01
メディア: 単行本
購入: 3人クリック: 6回
この商品を含むブログ (9件) を見る

猫も杓子も構造化

発達障害、特別支援などについて書いています。最近は心理学関係の内容が多めです。

二項分布の3次のモーメント

モーメントについて

続続・自由度について（実験修正編）

続・自由度について（別証明編）

続・自由度について（実験編）

自由度について

連続のための修正とは何か