猫も杓子も構造化

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

二項分布の4次のモーメントの導出

前回までで3次のモーメントをやったんだから次は4次だろう(安易)。

nekomosyakushimo.hatenablog.com
nekomosyakushimo.hatenablog.com

原点まわりのモーメントの導出

まず、モーメント母関数を3回微分したのものがこちら。

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


これを、(大変に面倒くさいけど)さらにもう1回微分すると、

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


このような形になる。次に、 \xi =0とおき、次数を揃えながら整理していくと、

  
\begin{split} f^{\prime\prime\prime\prime}(0) &= np \{ 4(n-1)p + 2(n-1)(n-2)p^2 + 3(n-1)(n-2)p^2  + (n-1)(n-2)(n-3) p^3 \\
\\
&\quad + 1 + (n-1) p +  2(n-1) p + (n-1)(n-2) p^2 \}\\
\\
&= np \{ (n-1)(n-2)(n-3) p^3 + [(n-1)(n-2) + 2(n-1)(n-2) + 3(n-1)(n-2) ] p^2 \\
\\
&\quad + [ 4(n-1) + 2(n-1) + (n-1)]p + 1 \}\\
\\
&= np \{ (n-1)(n-2)(n-3) p^3 + 6(n-1)(n-2)p^2 + 7(n-1)p + 1 \} \\
\\
&=  n(n-1)(n-2)(n-3) p^3 + 6n(n-1)(n-2)p^3 + 7n(n-1)p^2 + np 
\end{split}

が求まる。これが二項分布の4次の原点まわりのモーメントである。

平均値まわりのモーメントの導出

次に、平均値まわりでのモーメントを求めてみよう。求める値は、

   \displaystyle \begin{split} m_4 &= E\big( (x - \mu)^4 \big) = E(x^4 -4x^3 \mu + 6x^2\mu^2 -4x\mu^3 + \mu^4 )\\ \\
&= E(x^4)  -4\mu E(x^3) + 6\mu^2 E(x^2) -4\mu^3 E(x) + \mu^4 \\ \\
&= E(x^4)  -4\mu E(x^3) + 6\mu^2 E(x^2) -3\mu^4 \end{split}


である。各項に原点まわりのモーメントを代入すると、

   \begin{split} m_4 &= n(n-1)(n-2)(n-3)p^4 + 6n(n-1)(n-2)p^3 + 7n(n-1)(n-2)p^2 + np \\ \\
& \quad - 4np\{n(n-1)(n-2)p^3 + 3n(n-1)p^2 + np\} + 6n^2p^2 \{n(n-1) + np\} -3n^4p^4
\end{split}


となる。とりあえず強引に全部展開すると、

   \begin{split} m_4 &= n(n-1)(n-2)(n-3)p^4 + 6n(n-1)(n-2)p^3 + 7n(n-1)(n-2)p^2 + np \\ \\
& \quad - 4n^2(n-1)(n-2)p^4 -12n^2(n-1)p^3 -4n^2p^2 + 6n^3(n-1)p^4 +6 n^3p^3 -3n^4p^4
\end{split}


になり、それを pの次数で整理すると、

   \begin{split} m_4 &= (3n^2 -6n)p^4 + (-6n^2 +12n)p^3 + (3n^2 -7n)p^2 + np \end{split}


が得られる。3次までの平均まわりのモーメントでは、 npおよび (1-p)が登場していたことに注意をしつつ、これを因数分解していくと、

   \begin{split} m_4 &= (3n^2 -6n)p^4 + (-6n^2 +12n)p^3 + (3n^2 -7n)p^2 + np \\ \\
&= 3n^2p^4 - 6np^4 -6n^2p^3 + 12np^3 + 3n^2p^2 - 7np^2 + np \\ \\
&=np\{ 3np^3 - 6p^3 -6np^2 + 12p^2 + 3np - 7p + 1 \}\\ \\
&=np\{ 3np^3 - 6np^2 + 3np - 6p^3 + 12p^2 - 6p + 1 - p \}\\ \\
&=np\{ 3np(p^2 - 2p + 1) -6p (p^2 -2p + 1) + 1 - p \}\\ \\
&=np\{ 3np(1 - p)^2 -6p (1 - p)^2 + 1 - p \} \\ \\
&=np\{ (3np -6p)(1 - p)^2 + 1 - p \} \\ \\
&=np(1-p)\{(3np -6p)(1 - p) + 1 \} \\ \\
&=np(1-p)\{1+ 3(n - 2)p(1 - p)\} 
\end{split}


が求まる。これが、平均まわりの4次のモーメントである。非常に煩雑な計算であった。

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

前回の記事で、モーメント母関数を用いた二項分布の3次のモーメントの導出について書いた。今回は、二項分布から抽出したデータの分布を実際に見ながら3次のモーメントと分布の形にについてみてみる。

Rのデフォルトの関数には、歪度を計算するものはないらしい。momentsというパッケージを入れれば、その中にskewness()という歪度を計算する関数があるので今回はそれを使う。

まずは、成功確率 pが異なる二項分布から、それぞれ  n = 1000程度でデータを発生させる。

dat1 <- rbinom(1000,10,prob=0.5 )
dat2 <- rbinom(1000,10,prob=0.25)
dat3 <- rbinom(1000,10,prob=0.75)
dat4 <- rbinom(1000,10,prob=0.10)
dat5 <- rbinom(1000,10,prob=0.90)

それぞれのデータで歪度を求めるとともに、度数分布をプロットしてみる。

まずは、 p= 0.5のものから。

f:id:nekomosyakushimo:20170915225523p:plain

 a_3 = -0.0515で、0に近い。分布もそれに対応して概ね左右対称である。

続いて、 p= 0.25のものを見る。

f:id:nekomosyakushimo:20170915225731p:plain

 a_3 = 0.373である。歪度は正の値だと、右の裾が長くなると言われており、分布からも確認することができる。

今度は、 p = 0.75

f:id:nekomosyakushimo:20170915230344p:plain

 a_3 = -0.382であり、負の値なので左の裾が長くなっている。一つ前のものと比べると絶対値は似たような値であり、分布の歪みの程度が大体同じであることを示している。

より偏りのある p= 0.10で見てみよう。

f:id:nekomosyakushimo:20170915230750p:plain

 a_3 = 0.991であり、より大きな偏りを持ち、右の裾の長い分布であることが分かる。

最後は p= 0.90

f:id:nekomosyakushimo:20170915231059p:plain

 a_3 = - 0.899で、予想通りに左に裾が長い分布が確認できる。

歪度を求めると分布の形についてより詳細に分かるので、分布の形状を確かめるときには算出すると良いでしょう。

二項分布の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のときに分布は左右対称で、歪度が正のときには右の裾が長く、負のときには左の裾が長い分布が得られる。実際に二項分布の確率分布を確認してみてもその通りになっていることが分かる。

モーメントについて

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

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

離散値をとる確率変数 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

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

    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なことから、下から2番目の式の()内は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

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

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

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

続続・自由度について(実験修正編)

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

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

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

続・自由度について(別証明編)

前の記事までで、標本のデータから母分散を推定するには、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/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:色付きヒストグラムを重ねるときには ヒストグラム | を参照した。