読者です 読者をやめる 読者になる 読者になる

統計,確率のお勉強

統計学を主に勉強しています。勉強したことをアウトプットしていきます。 (※数式はMathJaxにより描画されています。ロードに少し時間がかかることがあります。)

Study Probability & Statistics

確率統計とアクチュアリーサイエンス

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

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}
$$

となる。

参考文献

特に無し

多変量解析1 多変量分布他

多変量分布

今回は多変量解析です。線形代数の知識が必要になってきて私は少し苦手
です...。
しかし今の時代、1変量でデータ解析なんて殆ど無いでしょうからちゃんとべんきょうしなきゃですなあ。

2変量の場合について

まずは2変量の場合について見ていきましょう。まず、確率変数(random variables = r.v.) X,Y を考えます。\forall x,y \in \mathbb{R} に対して、積分布関数(cumlative distribution function = c.d.f.)は次で定義されます。

$$
F(x,y) = Pr\{X\le x, Y\le y \}
$$

積分布関数が絶対連続(absolutely continuous)であるとき、偏微分がほとんどいたるところで存在し

絶対連続 - Wikipedia


$$
\frac{\partial^2F(x,y)}{\partial x \partial y} = f(x,y)
$$

及び

$$
F(x,y) = \int_{-\infty}^y \int_{-\infty}^x f(u,v)dudv
$$

が成り立つ。


f(x,y) \ge 0
\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(u,v)dudv = 1

p変量

次にp変量の場合を考える。 X_1,X_2,\cdots,X_p をr.v.とする。c.d.f.は

$$
F(x_1,x_2,\cdots,x_p) = Pr(X_1 \le x_1,X_2 \le x_2,\cdots,X_p \le x_p)
$$

F(x_1,x_2,\cdots,x_p) が絶対連続の時、密度関数(density function)は

$$
\frac{\partial^p F(x_1,x_2,\cdots,x_p)}{\partial x_1\partial x_2 \cdots\partial x_p} = f(x_1,x_2,\cdots,x_p)
$$

また、

$$
F(x_1,\cdots,x_p) = \int_{-\infty}^{x_p}\cdots\int_{-\infty}^{x_1}f(u_1,\cdots,u_2)du_1\cdots du_p
$$

周辺分布(Marginal Distribution)

再び2変量で見ていきます。確率変数 X,Y の累積分布関数(c.d.f.)が与えられた時 X周辺分布関数

\begin{eqnarray}
Pr\{X\le x\} &=& Pr\{X\le x,Y \le \infty\} \\
&=& F(x,\infty)
\end{eqnarray}

で与えられ、これを F(x) と表記する。また

\begin{eqnarray}
F(x) &=& \int_{-\infty}^{x}\int_{-\infty}^{\infty}f(u,v)dvdu \\
&=& \int_{-\infty}^{x}f(u)du
\end{eqnarray}

となる。 Y に対しても同様に求めることができる。

さて再びp変量について考えていきます。r.v. X_1,\cdots,X_p のc.d.f.として F(x_1,\cdots,x_p) が与えられたとする。この時、周辺分布は

\begin{eqnarray}
Pr\{X_1\le x_1,\cdots,X_r\le x_r\} &=& Pr\{X_1\le x_1,\cdots,X_r\le x_r,X_{r+1} \le \infty,\cdots,X_p\le \infty\} \\
&=& F(x_1,\cdots,x_r,\infty,\cdots,\infty)
\end{eqnarray}

ここで X_1,\cdots,X_r の周辺密度は

$$
\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}f(u_1,\cdots,u_p)du_{r+1}\cdots du_{p}
$$

で与えられる。

今日はここまで

まだ定義とか書いただけだけどここまでだな...勉強始めたばかりでまだ良く見えてこない...

参考文献

Anderson T.W.(1958)『An Introduction to Multivariate Statistical Analysis』 John Wiley & Sons

最小二乗法(単回帰)

理系の大学に入って、最初の年。物理学などの実験をかせられるところも多いだろう。その時、実験値に対して最小二乗法をしてグラフを書け!みたいなことを言われると思う。
私自身が物理学の実験を行っていた時も最小二乗法を使っていたが、何をしている課さっぱりわからなかった。実験前の授業で前準備として前で物理科の教授が高速で最小二乗法の導出を行っていたが当時はとりあえず実験データを与えられた式に当てはめて、その値をただ使っていただけだった。

ここでは、最小二乗法とは何がしたいのか。そしてその導出を行う。

記号の確認

最小二乗法に入る前に、各記号の定義?を確認しておく。

\begin{eqnarray}
\bar{x} &=& \frac{1}{n}\sum_{i=i}^n x_i \\
\bar{y} & =& \frac{1}{n}\sum_{i=1}^n y_i \\
s_x^2 &=& \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2 \\
s_y^2 &=& \frac{1}{n}\sum_{i=1}^n (y_i - \bar{y})^2 \\
r_{xy} &=& \frac{s_{xy}}{s_x s_y} = \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}}
\end{eqnarray}

上から、 x,y の平均、分散、 xy の相関係数である。

※相関係数は直線的関係を図る尺度である。

単回帰,重回帰,説明変数,従属変数

相関係数 r_{xy}\pm 1 に近く正(負)の相関が認められる時、一方の変数を他方の変数の1次関数として表すことができると考えられる。つまり、 yx の関数として考えた時、

$$
y = a + bx
$$

というモデルを考えることができる。このモデルを線形回帰モデル、このモデルを用いた分析を回帰分析という。

この時 x説明変数(独立変数) y従属変数(被説明変数)という。

  • 説明変数が1個 ・・・・・ 単回帰

  • 説明変数が2個以上 ・・・ 重回帰
  • となる。

    最小二乗法(単回帰)

    ここでは説明変数が一つの単回帰の最小二乗法について説明する。
    n 個の観測値 (x_i,y_i) \;\; (i = 1,2,\dots,n) を考える。この時、 n 個の観測値に対して、

    $$
    y_i = a + bx_i
    $$

    という回帰モデルを考える。観測値が一つの直線上に乗ることはありえないので、
    回帰式の誤差を小さくする a,b を推定値とすることになる。

    回帰式の誤差

    $$
    y_i -(a + bx_i)
    $$

    しかし、このまま足してしまうと、値に正負が存在するため、誤差の大きさを図ることができないそのため最小二乗法では上記の式を二乗したものをすべて足し合わせたものを用いる。つまり

    $$
    \mathcal{Q} = \sum_{i=1}^n\{y_i - (a + bx_i)\}^2
    $$

    これは非負値の二次式であり、これを a,b について偏微分した式を 0 とするような a,b\mathcal{Q} を最小化する a,b である。よって以下の連立方程式の解が求める推定値 \hat{a}, \hat{b} である。


    \begin{eqnarray}
    \left\{
    \begin{array}{l}
    \frac{\partial \mathcal{Q}}{\partial a} = -2\sum_{i=1}^n\{y_i - (a+bx_i)\} = 0 \\
    \frac{\partial \mathcal{Q}}{\partial b} = -2\sum_{i=1}^n x_i\{y_i - (a+bx_i)\} = 0 \\
    \end{array}
    \right.
    \end{eqnarray}

    \begin{eqnarray}
    \Leftrightarrow \left\{
    \begin{array}{l}
    \sum_{i=1}^n\{y_i - (\hat{a}+\hat{b}x_i)\} = 0 \\
    \sum_{i=1}^n x_i\{y_i - (\hat{a}+\hat{b}x_i)\} = 0 \\
    \end{array}
    \right.
    \end{eqnarray}

    \begin{eqnarray}
    \Leftrightarrow \left\{
    \begin{array}{l}
    \sum_{i=1}^n y_i = \hat{a}n + \hat{b}\sum_{i=1}^n x_i \\
    \sum_{i=1}^n x_i y_i = \hat{a}\sum_{i=1}^n x_i + \hat{b}\sum_{i=1}^n x_i^2
    \end{array}
    \right.
    \end{eqnarray}

    \begin{eqnarray}
    \Leftrightarrow \left\{
    \begin{array}{l}
    \bar{y}= \hat{a} + \hat{b}\bar{x} \\
    \frac{1}{n}\sum_{i=1}^n x_i y_i = \hat{a}\bar{x} + \hat{b}\frac{1}{n}\sum_{i=1}^n x_i^2
    \end{array}
    \right.
    \end{eqnarray}

    \begin{eqnarray}
    \Leftrightarrow \left\{
    \begin{array}{l}
    \hat{a} = \bar{y} - \hat{b}\bar{x} \\
    \frac{1}{n}\sum_{i=1}^n x_i y_i - \bar{x}\bar{y} = \hat{b}(\frac{1}{n}\sum_{i=1}^n x_i^2 - \bar{x}^2)
    \end{array}
    \right.
    \end{eqnarray}

    \begin{eqnarray}
    \Leftrightarrow \left\{
    \begin{array}{l}
    \hat{a} = \bar{y} - \hat{b}\bar{x} \\
    \hat{b} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n(x_i - \bar{x})^2} = \frac{s_{xy}}{s_{x^2}} = \frac{r_{xy}s_y}{s_x}
    \end{array}
    \right.
    \end{eqnarray}


    したがって推定された回帰直線は

    $$
    \hat{y} = \hat{a} + \hat{b}x = \bar{y} + \hat{b}(x-\bar{x}) = \bar{y} + (\frac{s_{xy}}{s_{x^2}})(x-\bar{x})
    $$

    これを推定回帰直線または標本回帰直線という。これを式変形すると

    $$
    (\frac{\hat{y} - \bar{y}}{s_y}) = r_{xy}(\frac{x-\bar{x}}{s_x})
    $$

    とかける。この形は覚えやすいので、覚えておくといいだろう。

    まとめ

    最小二乗法の考え方の基本は誤差を小さくすることにある。この点だけ気をつけていれば、後の操作は極めて当然のものだと思えるだろう。最後にもう一度推定値を示しておく。

    \begin{eqnarray}
    \left\{
    \begin{array}{l}
    \hat{a} = \bar{y} - \hat{b}\bar{x} \\
    \hat{b} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n(x_i - \bar{x})^2} = \frac{s_{xy}}{s_{x^2}} = \frac{r_{xy}s_y}{s_x}
    \end{array}
    \right.
    \end{eqnarray}

    参考文献

    高澤俊幸・勝野健太郎・服部真・山内恒人(2005)『モデリング』社団法人 日本アクチュアリー会.
    稲垣宣生(2013)『数理統計学』(数学シリーズ)裳華房.

    【デレステ】ありすPが33万円分1000連越えの爆死をしたそうですが、検定してみたいと思う。

    デレステのガシャ確率は本当に正しいのか?消費者には見えないガシャ確率を統計的に考える

    twitterを眺めていたら、最近よく見るソシャゲ爆死記事を見つけました。今回目にしたのは1000連越えのデレステガシャをして目的のキャラ(橘ありす)のSSレアが出なかったということ。

    多分これが目的のヤツです。

    いましたいました。

    提供割合を見てみると...

    SSレアは1.5%、SSレアの種類は26種類。

    ほうほう....

    となると、目的のありすが当たる確率は単純計算で


    \frac{15}{1000}×\frac{1}{26} = \frac{15}{26000}(\simeq 0.0005769)

    と考えられます。

    ここでは、何連して全くでなければ、この示されてる確率が怪しくなってくるか...

    つまりは、何連して目的のキャラが出なかったら運営に文句言えるのか、を、検定していきます。


    先の計算から、運営が示している(?)「はじめての表情」橘ありすが1回のガシャで当たる確率は

    $$
    p_0=\frac{15}{26000}
    $$

    になります。私たちは本当にこの確率なの?もっと低いんじゃないの?と疑ってるわけですので以下のように検定問題を考えます。

    \begin{eqnarray}
    \left\{
    \begin{array}{l}
    H_0 : p = p_0 = \frac{15}{26000} \\
    H_1 : p < p_0
    \end{array}
    \right.
    \end{eqnarray}

    H_0は帰無仮説、H_1は対立仮説になります。これを有意水準5%で検定していきます。

    統計量は

    $$
    T = \frac{\hat{p}-p_0}{\sqrt{\frac{p_0(1-p_0)}{n}}}
    $$

    になります。

    ここで、棄却域を考える。
    標本数は大きいと考えるのでこの統計量は標準正規分布に近似的に従う。
    この検定は片側検定なので、有意水準5%で棄却域

    $$
    T<-z(0.05)
    $$

    となる。つまり、統計量Tが標準正規分布左側5%の値よりも小さいなら帰無仮説は棄却され、運営側が提示している確率より、実際の確率は低いと言える。

    これを用いて、何回引いて出ないならば運営に文句が言えるのかを考える。

    n回ガシャして、1回もありすは出ないので標本比率は\hat{p}=0となる。

    これを用いて実現値T^*

    $$
    T^*= \frac{0-p_0}{\sqrt{\frac{p_0(1-p_0)}{n}}}
    $$

    これに、他の値も代入して先の不等式をnについて解くと

    \begin{eqnarray}
    \frac{0-p_0}{\sqrt{\frac{p_0(1-p_0)}{n}}} &<& -z(0.05) \\
    n &>& z(0.05)^2\frac{1-p_0}{p_0} \\
    \therefore n &>& 4659.3
    \end{eqnarray}

    となり、標本数が4660以上の時、帰無仮説は棄却されます。

    以上から、SSレア各キャラの提供割合が等しいと仮定して計算した場合、上記の回数、4660回ガシャをして目的のカードが出なかったらまず、運営の表示している提供割合は正しいと考えることは統計的にはできません。

    つまり、今回のありすPさんは1000連越えのガシャを回して出ていませんが、この試行回数ではまだ運が悪かったと言わざるを得ません。

    f:id:doratai:20160506225904p:plain

    ありすのガシャ出現確率は上がっているそうですが...それでもほんのすこしでしょう...。4000ガシャは覚悟したほうがいいのかもしれません。

    現状、デレステ運営に文句をいうことはできないということになります。

    さてところで、デレステ運営に文句を言える4660回ガシャをするにはどれくらいお金がかかるのでしょうか....

    すべて有償ジュエルで賄うとして考えます。

    上の画像を見てみると8400個9800円がおそらくいちばん単価が小さいのだと思われます。(PS4の最新ゲーム1つ買える...)
    一回のガシャに必要なスタージュエルは250個。4660回ガシャしたいので

    $$
    250 \times 4660 = 1,165,000
    $$

    11165000個のスタージュエルが必要になります。これを満たす数8400個9800円のスタージュエルのセット数は

    $$
    11165000/8400 \simeq 138.7
    $$

    より139セット買う必要がある。

    つまり、運営に文句を言える4660回ガシャを引くためには

    $$
    9800 \times 139 = 1,362,200
    $$

    より約136万2200円の課金が必要になってきます。(うわあ....)

    100万以下の課金はぬるいと...そういうことですかね...
    目的のキャラを当てようと思って課金するのは完全に泥沼ルートです...

    ソシャゲの闇は深い。