default eye-catch image.

分散と標準偏差を計算しやすくする

[mathjax] 分散と標準偏差を計算しやすく変形できる。 いちいち偏差(x_i-bar{x})を計算しておかなくても、2乗和(x_i^2)と平均(bar{x})がわかっていればOK。 begin{eqnarray} s^2 &=& frac{1}{n} sum_{i=1}^n (x_i - bar{x})^2 \\ &=& frac{1}{n} sum_{i=1}^n ( x_i^2 -2 x_i bar{x} + bar{x}^2 ) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 -2 bar{x} sum_{i=1}^n x_i + bar{x}^2 sum_{i=1}^n 1) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 -2 n bar{x}^2 + nbar{x}^2 ) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 - nbar{x}^2 ) \\ &=& frac{1}{n} sum_{i=1}^n x_i^2 - bar{x}^2 \\ s &=& sqrt{frac{1}{n} sum_{i=1}^n x_i^2 - bar{x}^2 } end{eqnarray} 以下みたい使える。平均と標準偏差と2乗和の関係。 begin{eqnarray} sum_{i=1}^n (x_i - bar{x})^2 &=& sum_{i=1}^n x_i^2 - nbar{x}^2 \\ ns^2 &=& sum_{i=1}^n x_i^2 - nbar{x}^2 \\ sum_{i=1}^n x_i^2 &=& n(s^2 + bar{x} ) end{eqnarray}

default eye-catch image.

微分フィルタだけで時系列データの過渡応答終了を検知したい

ある値の近傍を取る状態から別の値の近傍を取る状態へ遷移する時系列データについて、 どちらの状態でもない過渡状態を経て状態が遷移するみたい。 今回知りたいのは理論ではなく、定常状態に遷移したことをいかにして検知するかという物理。 まぁいかようにでも検知できなくもないけれども、非定常状態の継続時間は不明であるという。 なるべく遅れなしに定常状態への遷移を検知したい。 ARモデルとかMAモデルとかの出番ではない。 DeepLearningとかやってる暇はないw 過去の微小時間内に含まれるデータだけから次が定常状態なのかを検知したい。 前状態と後状態のベースラインがどんな値であっても、 それを考慮することなしに過渡応答の後にくる定常状態の先頭を検知できる、という特徴がほしい。 移動平均だとベースラインがケースごとに変わるので扱いづらい。 微分であればベースラインがゼロに揃うのが良いな、と思った。 結局、今度はゼロへの収束を検知する問題に差し代わるだけだけども、 前ラインと後ラインの2つも見ないでも、確実に楽にはなっている。 微分するとノイズが増えるが...。 微分値が上下の閾値を超えたラインから先、 ゼロ収束の判定条件が一定時間継続した最初のデータポイントがソコ。 これは過渡状態の終了検知なだけだから、それとRawデータの変化をみる。

default eye-catch image.

やってみた Markov chain Monte Carlo methods, MCMC , gibbs sampling

[mathjax] マルコフ連鎖モンテカルロ法。2変量正規分布からGibbs Samplingする方法を考えてみました。 式を流してもよくわからないので、行間ゼロで理解できるまで細切れにして書いてみます。 Gibbs Sampling 無茶苦茶一般的に書かれている書籍だと以下みたいになっている。 ステップごとに(theta=(theta_1,theta_2,cdots,theta_n))の各要素を、 その要素以外の要素の条件付き確率でランダムに発生させる。 begin{eqnarray} theta^0 &=& (theta_1^0,theta_2^0,cdots,theta_n^0) \\ theta_i^{t+1} &sim& p(theta_i|theta_1^{t+1},cdots,theta_{i-1}^{t+1},theta_{i+1}^t,cdots,theta_n^t) end{eqnarray} 2次元空間であれば、ある点を決める際に、 片方の変数(theta_1)を固定して条件付き確率(P(theta_2|theta_1))を最大にするパラメタ(hat{theta_2})を決める。 その時点で((theta_1,hat{theta_2}))が決まる。 次は変数(hat{theta_2})を固定して条件付き確率(P(theta_1|hat{theta_2}))を最大にするパラメタ(hat{theta_1})を決める。 2次元正規分布であれば、片方の確率変数を固定したときの条件付き確率が1変数の正規分布になるから、 これがすごくやりやすい。 Gibbs Samplingは、このように条件付き確率を計算できる必要がある。 2変量正規分布の条件付き確率 2変量正規分布の片方の確率変数を固定すると1変数の正規分布が出てきます。 山の輪切りにして出てきた正規分布の母数を調べたいのですが、導出が無茶苦茶大変そうです。 出てきた正規分布の母数について式変形すると出てくるようです。こちらを参考にさせて頂きました。 この関係式を利用してひたすら式変形していきます。 begin{eqnarray} f(x|y) &=& frac{f(x,y)}{f(x)} end{eqnarray} 今、確率変数(y=y\')を固定し確率変数(x)の正規分布を考えます。 2変量正規分布の分散共分散行列が(sum)であるとします。 begin{eqnarray} sum &=& begin{pmatrix} sigma_{xx} & sigma_{xy} \\ sigma_{yx} & sigma_{yy} end{pmatrix} end{eqnarray} 出てくる正規分布の平均は以下となります。 確率変数が(x)なのに固定した(y)が出てくるのがポイント。 begin{eqnarray} mu\' = bar{x} + frac{sigma_{xy}}{sigma_{yy}} (y\'-bar{y}) end{eqnarray} また、出てくる正規分布の分散は以下となります。 begin{eqnarray} sigma\'^2 &=& sigma_{xx} - frac{sigma_{xy}}{sigma_{yy}} sigma_{yx} end{eqnarray} Python実装例 実際にPythonコードにしてみた図です。2変量正規分布の条件付き確率さえわかってしまえば あとはコードに落とすだけです。 import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats n_dim = 2 def gibbs_sampling(mu, sigma, sample_size): samples = [] start = [0, 0] samples.append(start) search_dim = 0 for i in range(sample_size): search_dim = 0 if search_dim == n_dim-1 else search_dim + 1 prev_sample = samples[-1][:] s11 = sigma[search_dim][search_dim] s12 = sigma[search_dim][search_dim - 1] s21 = sigma[search_dim -1 ][search_dim] s22 = sigma[search_dim - 1][search_dim - 1] mu_x = mu[search_dim] mu_y = mu[search_dim-1] _y = prev_sample[search_dim - 1] new_mean = mu_x + s12/float(s22) * (_y - mu_y) new_sigma = s11 - s12/float(s22) * s21 sample_x = np.random.normal(loc=new_mean, scale=np.power(new_sigma, .5), size=1) prev_sample[search_dim] = sample_x[0] samples.append(prev_sample) return np.array(samples) nu = np.ones(2) covariance = np.array([[0.5, 0.5], [0.5, 3]]) samples = gibbs_sampling(nu, covariance, 1000) fig, ax1 = plt.subplots(figsize=(6, 6)) ax1.scatter(sample[:, 0], sample[:, 1], marker=\"o\", facecolor=\"none\", alpha=1., s=30., edgecolor=\"C0\", label=\"Samples\" ) 実行結果

default eye-catch image.

最尤推定とベイズの定理とMAP推定

[mathjax] 最尤推定とMAP推定とベイズの定理は繋がっていたので、 記憶が定かなうちに思いの丈を書き出してみるテスト。俯瞰してみると面白い。 あるデータ達(x)が観測されていて、それらは未知のパラメータを持つ確率分布から発生している。 観測されたデータ達(x)を使って、それらのデータを発生させたモデルのパラメータを推定したい。 確率密度関数の中に2つの変数があって。 片方を定数、片方を確率変数として扱うことで2通りの見方ができる。 例えば(n)回のコイントスで(k)回表が出る確率が(theta)だとしてベルヌイ分布の確率密度関数は (k)と(theta)のどちらが確率変数だとしても意味がある。 begin{eqnarray} f(k;theta)=theta^k (1-theta)^{n-k} end{eqnarray} 表が出る確率(theta)が定数だと思って、確率変数(x)の確率密度関数と思う。 単なる確率変数(x)の確率密度関数の中に(theta)という定数がある。尤度。 尤度は確率変数(x)の確率密度関数!。 begin{eqnarray} p(X=x|theta) end{eqnarray} 尤度(p(x|theta))を最大にする(theta)を推定するのが最尤推定。 begin{eqnarray} newcommand{argmax}{mathop{rm arg~max}limits} hat{theta} = argmax_{theta} p(X=x|theta) end{eqnarray} 事後確率と事前確率には関係があって以下のようになる。ベイズの定理。 begin{eqnarray} p(theta|x)=frac{p(x|theta) p(theta)}{p(x)} end{eqnarray} ちなみに、(p(x))は以下のようにしておくとわかりやすい。 同時確率と周辺確率の関係。表を書いて縦、横がクロスするところが同時確率だとして、 縦、横いずれかの方向に同時確率を足し合わせる操作にあたるらしい。 なにか確率変数が独立でなければならない、というのは気にしない。 begin{eqnarray} p(x) = int p(x,theta) dtheta end{eqnarray} なので以下みたいに書き直せる。最後の比例のところは...。 左辺は事後確率分布、右辺は尤度と事前確率分布の積!!。 begin{eqnarray} p(theta|x) &=& frac{p(x|theta) p(theta)}{p(x)} \\ &=& frac{p(x|theta) p(theta)}{int p(x,theta) dtheta} \\ &propto& p(x|theta) p(theta) end{eqnarray} (p(theta))は確率変数(theta)の確率分布。尤度(p(x|theta))は(theta)に対して定数。(x)に対して変数。 ということで、右辺は確率分布(p(theta))を尤度(p(x|theta))を使って変形した確率分布。 で、左辺の(p(theta|x))は右辺の(p(x|theta) p(theta))を定数倍した確率分布。 データを観測していない状態で立てた(p(theta))があって、 観測したデータを使って求めた尤度(p(x|theta))が得られたことで、 左辺の(p(theta|x))が得られた、という状況。 (p(theta|x))は確率変数(theta)の確率分布なので、 最尤推定とベイズの定理を俯瞰してみると、最尤推定が点推定である一方で、 ベイズの定理では確率分布が得られるという具合で異なる。 (観測値が極端なデータだったとき、最尤推定は極端な推定結果が得られるだけだけれども、 ベイズの定理で得られる事後確率分布は確率分布なので様子がわかる..??) 事後確率分布を最大化する(theta_{MAP})を求めるのがMAP推定。(点推定) begin{eqnarray} hat{theta_{MAP}} = argmax_{theta} p(theta|x) end{eqnarray} 尤度をかけて得られた事後分布と同じ形になる便利な分布があって、 観測データ達の分布と対応して決まっている(共役事前分布)。 ベルヌイ分布の共役事前分布はベータ分布。

default eye-catch image.

正規分布に従う確率変数の二乗和はカイ二乗分布に従うことを実際にデータを表示して確かめる

以前、\"正規分布に従う確率変数の二乗和はカイ二乗分布に従うことの証明\"という記事を書いた。 記事タイトルの通り、正規分布に従う確率変数の二乗和はカイ二乗分布に従う。 [clink url=\"https://ikuty.com/2018/08/01/chi-square-distribution\"] 実際にデータを生成して確かめてみる。 まずは、scipi.stats.chi2.pdfを使って各自由度と対応する確率密度関数を書いてみる。 \"pdf\" ってPortableDocumentFormatではなく、ProbabilityDensityFunction(確率密度関数)。 import numpy as np import matplotlib.pyplot as plt from scipy import stats # 0から8まで1000個のデータを等間隔で生成する x = np.linspace(0, 8, 1000) fig, ax = plt.subplots(1,1) linestyles = [\':\', \'--\', \'-.\', \'-\'] deg_of_freedom = [1, 2, 3, 4] for k in deg_of_freedom: linestyle = linestyles[k-1] ax.plot(x, stats.chi2.pdf(x, k), linestyle=linestyle, label=r\'$k=%i$\' % k) plt.xlim(0, 8) plt.ylim(0, 1.0) plt.legend() plt.show() 定義に従って各自由度ごとに2乗和を足し合わせてヒストグラムを作ってみる。 (ループのリストが自由度) cum = 0 for i in [1,2,3,4]: # 標準正規分布に従う乱数を1000個生成 x = np.random.normal(0, 1, 1000) # 2乗 x2 = x**2 cum += x2 plt.figure(figsize=(7,5)) plt.title(\"chi2 distribution.[k=1]\") plt.hist(cum, 80) 全部k=1ってなってしまった...。左上、右上、左下、右下の順に1,2,3,4。 頻度がこうなので、各階級の頻度の相対値を考えると上記の確率密度関数の形になりそう。 k=1,2が大きくことなるのがわかるし、3,4と進むにつれて変化が少なくなる。

default eye-catch image.

標本調査に必要なサンプル数の下限を与える2次関数

[mathjax] 2項分布に従う母集団の母平均を推測するために有意水準を設定して95%信頼区間を求めてみた。 母平均のあたりがついていない状況だとやりにくい。 [clink url=\"https://ikuty.com/2019/01/11/sampling/\"] (hat{p})がどんな値であっても下限は(hat{p})の関数で抑えられると思ったので、 気になって(hat{p})を変数のまま残すとどうなるかやってみた。 begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 \\ frac{1.96}{0.05}sqrt{hat{p}(1-hat{p})} le sqrt{n} \\ 39.2^2 hat{p}(1-hat{p}) le n end{eqnarray} 左辺を(f(hat{p}))と置くと (f(hat{p}))は下に凸の2次関数であって、 (frac{d}{dhat{p}}f(hat{p})=0)の時に最大となる。というか(hat{p}=0.5)。 (hat{p}=0.5)であるとすると、これはアンケートを取るときのサンプル数を求める式と同じで、 非常に有名な以下の定数が出てくる。 begin{eqnarray} 1537 * 0.5 (1-0.5) le n \\ 384 le n end{eqnarray} (hat{p})がどんな値であっても、サンプル数を400とれば、 有意水準=5%の95%信頼区間を得られる。 だから、アンケートの(n)数はだいたい400で、となる。 さらに、有意水準を10%にとれば、(n)の下限は100で抑えられる。 なるはやのアンケートなら100、ちゃんとやるには400、というやつがこれ。

default eye-catch image.

標本調査に必要なサンプル数を素人が求めてみる。

[mathjax] ちょっと不思議な計算をしてみる。 仮定に仮定を積み重ねた素人の統計。 成功か失敗かを応答する認証装置があったとする。 1回の試行における成功確率(p)は試行によらず一定でありベルヌーイ試行である。 (n)回の独立な試行を繰り返したとき、成功数(k)を確率変数とする離散確率変数に従う。 二項分布の確率密度関数は以下の通り。 begin{eqnarray} P(X=k)= {}_n C_k p^k (1-p)^{n-k} end{eqnarray} 期待値、分散は、 begin{eqnarray} E(X) &=& np \\ V(X) &=& np(1-p) end{eqnarray} (z)得点(偏差値,つまり平均からの誤差が標準偏差何個分か?)は、 begin{eqnarray} z &=& frac{X-E(X)}{sigma} \\ &=& frac{X-E(X)}{sqrt{V(X)}} \\ &=& frac{X-np}{sqrt{np(1-p)}} end{eqnarray} であり、(z)は標準正規分布に従う。 これを標本比率(hat{p}=frac{X}{n})を使うように式変形する。 begin{eqnarray} z &=& frac{frac{1}{n}}{frac{1}{n}} frac{X-np}{sqrt{np(1-p)}} \\ &=& frac{frac{X}{n}-p}{sqrt{frac{p(1-p)}{n}}} \\ &=& frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} end{eqnarray} (n)が十分に大きいとき、(z)は標準正規分布(N(0,1))に従う。 従って、(Z)の95%信頼区間は以下である。 begin{eqnarray} -1.96 le Z le 1.96 end{eqnarray} なので、 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{p(1-p)}{n}}} le 1.96 end{eqnarray} (hat{p})は(p)の一致推定量であるから、(n)が大なるとき(hat{p}=p)とすることができる。 begin{eqnarray} -1.96 le frac{hat{p}-p}{sqrt{frac{hat{p}(1-hat{p})}{n}}} le 1.96 \\ end{eqnarray} (p)について解くと(p)の95%信頼区間が求まる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} le p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} end{eqnarray} 上記のにおいて、標準誤差(1.96sqrt{frac{hat{p}(1-hat{p})}{n}})が小さければ小さいほど、 95%信頼区間の幅が狭くなる。この幅が5%以内であることを言うためには以下である必要がある。 (有意水準=5%) begin{eqnarray} 1.96sqrt{frac{hat{p}(1-hat{p})}{n}} le 0.05 end{eqnarray} 観測された(hat{p})が(0.9)であったとして(n)について解くと、 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.05 \\ frac{1.96}{0.05} sqrt{0.09} le sqrt{n} \\ 11.76 le sqrt{n} \\ 138.2 le n end{eqnarray} 139回試行すれば、100回中95回は(p)は以下の95%信頼区間に収まる。 つまり95%信頼区間は以下となる。 begin{eqnarray} hat{p}-1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} &le& p le hat{p}+1.96 sqrt{frac{hat{p}(1-hat{p})}{n}} \\ 0.9-1.96 frac{sqrt{0.09}}{sqrt{139}} &le& p le 0.9 + 1.96 frac{sqrt{0.09}}{sqrt{139}} \\ 0.9-1.96 frac{0.3}{11.78} &le& p le 0.9+1.96 frac{0.3}{11.78} \\ 0.85 &le& p le 0.95 end{eqnarray} (n)を下げたい場合は有意水準を下げれば良い。 統計的に有意水準=10%まで許容されることがある。 有意水準が10%であるとすると、(n)は35以上であれば良いことになる。 begin{eqnarray} 1.96sqrt{frac{0.9(1-0.9)}{n}} le 0.1 \\ frac{1.96}{0.1} sqrt{0.09} le sqrt{n} \\ 5.88 le sqrt{n} \\ 34.6 le n end{eqnarray} 信頼区間と有意水準の式において(p)を標本から取ってきたけど、 アンケートにおいてYes/Noを答える場合、(p)は標本における最大値(つまり0.5)を 設定して(n)を求める。 つまり、(p)として利用するのは標本比率ではないのかな?と。 このあたり、(hat{p})を変数として残すとどういうことがわかった。 [clink url=\"https://ikuty.com/2019/01/13/sampling_with2/\"]

default eye-catch image.

標本平均の母平均の推定

[mathjax] zozoスーツではないけれども、標本が正規分布に従うというのは、 真実の母平均に対して正規分布に従う計測誤差を含む分布を観測しているのと同じ。 母平均(mu)が未知である現象を計測誤差がわかっている計測手段で計測する話。 (n)回計測を行って得られた標本(X_1,X_2,cdots,X_n)は、母平均を中心として誤差分振れているはず。 つまり、(X_1=mu+e_1,X_2=mu+e_2,cdots,X_n=mu+e_N)。 誤差(e_i)は平均0、分散(sigma^2)の正規分布(N(0,sigma^2))に従うと考えると、 標本(X_i)は(mu)だけオフセットした正規分布(N(mu,sigma^2))に従うと考えられる。 標本平均は(bar{X}=frac{1}{n}(X_1+X_2+cdots+X_n))だから、 (n)に関係なく(E(bar{X})=mu)であって、(V(bar{X})=frac{sigma^2}{n})から(lim_{nrightarrow infty}V(bar{X})=0)。 中心極限定理により大なる(n)のとき(bar{X})の分布は正規分布(N(mu,frac{sigma^2}{n}))で近似できる。 標本(X_i)の標準偏差が(sigma)である一方、標本平均の標準偏差は(frac{sigma}{sqrt{N}})だから、 標本の分布より、標本平均の分布の方が裾が狭い。 正規分布(N(mu,frac{sigma^2}{n}))を標準化しておくと、標準正規分布の累積度数表を使って 平均(mu)、標準偏差(sigma)を評価できるようになる。z得点は以下の変換。 begin{eqnarray} Z=frac{bar{X}-mu}{frac{sigma}{sqrt{n}}} end{eqnarray} 分布(Z)は平均0、標準偏差1の標準正規分布になる。 見方としては、残差が標準偏差何個分か?の分布。全部足して1になる。 (bar{X},mu,sigma,n)として具体的な値を入れると数値(Z)が決まる。 ちなみに確率密度関数と累積度数は以下の通り。 begin{eqnarray} f(x) &=& frac{1}{sqrt{2pi}} exp left( -frac{x^2}{2} right) \\ int_{-infty}^{infty} f(x) dx &=& 1 end{eqnarray} (x=0)から(x=z)の面積(int_0^{z} frac{1}{sqrt{2pi}} left( -frac{x^2}{2} right) )を(Phi(z))とおき、 (Phi(z)=a)となる点を上側(a)パーセント点という名前が付いている。 (Phi(z))の積分は解析的に計算できないけれど、有用だし決まった数値なので、 ここみたいに表ができているからルックアップすれば良い様子。 (Z)得点が1.96であったとすると、標準正規分布表から(Phi(z=1.96)=0.475)であることがわかる。 これは上側確率が0.475という意味なので、両側確率は2をかけて0.975ということになる。 逆に言うと、(mu)だけが不明で、既知の母分散と標本平均から(mu)を推測することに、 この話を使うことができる。つまり、(-1.96 le z le +1.96)という式を立てると、 (mu)の信頼区間を作ることができる。つまり、(n)個の標本を取る操作を100回繰り返すと97.5回は 信頼区間が母平均を含まない区間になっている。 例 確率変数(X)が平均2、分散10の正規分布(N(2,10))に従うとする。 95%信頼区間は(-1.96 lt z lt 1.96)から、 (-1.96 sqrt{10} + 2 lt X lt 1.96 sqrt{10} + 2)。 (-4.2 lt X lt 8.2)。 100回試すと97.5回は母平均がこの区間にある。 (X)が負になる確率は、(Z=frac{X-2}{sqrt{10}})から、(sqrt{10}Z+2lt 0)、(Z lt -frac{2}{sqrt{10}})、(Z lt - 0.633)。 (P(X lt 0)=P(Z lt -0.633)=1-P(z lt 0.633))。

default eye-catch image.

たたみこみと正規分布の再生性

[mathjax] 正規母集団からの推定をやる前に、正規分布の再生性の理解が必要だったのでまとめてみる。 独立な確率変数(X_1)、(X_2)がそれぞれ確率分布(g(x))、(h(x))に従うとする。 各確率変数の和(X_1+X_2)が従う確率分布を(k(z))とする。 確率(P(X_1+X_2=z))を考えると、(X_1+X_2=z)となるのは、 (X_1=x, X_2=z-x)としたとき、両者を足して(z)になる全ての組み合わせ。 (X_1)は(g(x))、(X_2)は(h(z-x))に従うので、両者が同時に起こる確率は(g(x)h(z-x))。 これをまとめて書くと、 begin{eqnarray} k(z) = sum_x g(x)h(z-x) end{eqnarray} この形が「たたみこみ(convolution)」。  (k = g * x)と書く。 確率変数(X_1)、(X_2)が独立で、それぞれ平均(mu_1)、(mu_2)、分散(sigma_1^2)、(sigma_2^2)の正規分布に従うなら、 以下が成り立つ。 begin{eqnarray} N(mu_1,sigma_1^2) * N(mu_2,sigma_2^2) = N (mu_1+mu_2, sigma_1^2 + sigma_2^2) end{eqnarray} これ、モーメント母関数を使って証明できる様子。 ある分布のモーメント母関数があったとして、モーメント母関数を(n)回微分して変数を(0)と置くと、 分布の期待値、分散、歪度、突度など統計量を求められるやつ。 [clink url=\"https://ikuty.com/2018/09/22/moment_generating_funuction/\"] 正規分布の確率密度関数とモーメント母関数は以下の通り。 begin{eqnarray} f(x) &=& frac{1}{sqrt{2pisigma}} expleft( - frac{(x-mu)^2}{2sigma^2} right) \\ M(t) &=& exp left( mu t + frac{sigma^2 t^2}{2} right) end{eqnarray} もちろん、(N(mu_1,sigma_1^2))、(N(mu_2,sigma_2^2))のモーメント母関数は, begin{eqnarray} M_1(t) &=& exp left( mu_1 t + frac{sigma_1^2 t^2}{2} right) \\ M_2(t) &=& exp left( mu_2 t + frac{sigma_2^2 t^2}{2} right) \\ end{eqnarray} かけると、以下の通り(N(mu_1+mu_2,sigma_1^2+sigma_2^2))のモーメント母関数となる。 begin{eqnarray} M_1(t) M_2(t) &=& expleft( mu_1 t +frac{sigma_1^2 t^2}{2} right) expleft( mu_2 t + frac{sigma_2^2 t^2}{2} right) \\ &=& expleft( (mu_1+mu_2) t +frac{(sigma_1^2 + sigma_2^2) t^2}{2} right) end{eqnarray} たたみこみの操作は、独立な確率変数(X_1,X_2)について(X_1+X_2)の確率分布を求める操作だから、 この結果は独立な確率変数(X_1,X_2)が(N(mu_1,sigma_1^2))、(N(mu_2,sigma_2^2))に従うとき、 (X_1+X_2)が(N(mu_1+mu_2,sigma_1^2+sigma_2^2))に従うことを意味する。 ある確率分布のたたみ込みの結果が同じ確率分布になることを再生性(reproductive)というらしい。 正規分布の再生性を使った演算 正規分布には再生性があるので、以下みたいな演算ができる。 (X_1,X_2,cdots,X_n)が独立で、それぞれ正規分布(N(mu_1,sigma_1^2),N(mu_2,sigma_2^2),cdots,N(mu_N,sigma_N^2) )に 従うとき、(X_1+X_2+cdots+X_n)は(N(mu_1+mu2+cdots,mu_N,sigma_1^2+sigma_2^2+cdots+sigma_N^2))に従う。 (X_1,X_2,cdots,X_n)が全て同じ(N(mu,sigma^2))に従うなら、 (X_1+X_2+cdots+X_n)は、(N(nmu, nsigma^2))に従う。 (bar{X}=frac{X_1+X_2+cdots+X_n}{n})は(N(mu,frac{sigma^2}{n}))に従う。

default eye-catch image.

標本分散(sample variance)と不偏分散(unbiased variance)

[mathjax] 不偏分散は(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)。 分母から1を引く必要がある。なんでか調べてみたので書いてみる。 標本平均は(n)の大小によらず母平均の近傍にあって、母平均に確率収束する。 標本平均は(n)の大小に関係なく、その期待値と母平均が等しい(不偏)。 begin{eqnarray} E(bar{X}) &=& frac{1}{n}nmu = mu \\ lim_{n rightarrow infty} V(bar{X}) &=& 0 end{eqnarray} 100個のデータがあって、その中から5個取ったときの平均と、50個取ったときの平均に 母平均の推測という意味で違いがない。 では、分散はどうか。 定義通り標本の分散を(S^2 = frac{1}{n}{ (X_1-bar{X})^2 + (X_2-bar{X})^2 + cdots + (X_n-bar{X})^2 } )とすると、 (S^2)は母分散と等しくならない。不偏にならない。つまり、(E(S^2) ne sigma^2)。 その値が不偏であるか否かは、実際に期待値を式変形してみるとわかる。 結論を知っていないと出来ない変形ばかりだけども...。 begin{eqnarray} E(S^2) &=& Eleft[frac{1}{n} sum_{i=1}^n (x_i-bar{X})^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n left( (x_i-mu)-(bar{X}-mu) right)^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 -2sum_{i=1}^n(x_i-mu)(bar{X}-mu) + sum_{i=1}^n (bar{X}-mu)^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 -2n (bar{X}-mu) +n (bar{X}-mu)^2 right] \\ &=& frac{1}{n} Eleft[ sum_{i=1}^n (x_i-bar{X})^2 - n(bar{X}-mu)^2 right] \\ &=& frac{1}{n} sum_{i=1}^n Eleft[ (x_i-mu)^2 right] - Eleft[ (bar{X}-mu)^2 right] \\ &=& frac{1}{n} sum_{i=1}^n V(x_i) - V(bar{X}) \\ &=& sigma^2 - frac{1}{n} sigma^2 \\ &=& frac{n-1}{n} sigma^2 end{eqnarray} ということで、(E(S^2)ne sigma^2)。不偏でない。 では、どうすれば不偏な標本分散を得られるのか。 (E(S^2)=frac{n-1}{n} sigma^2)から、(frac{n}{n-1}E(S^2)=sigma^2)なので、(s^2=frac{n}{n-1}E(S^2))とすれば、 (s^2=sigma^2)ということになり、(s^2)は不偏となる。(s^2)を不偏分散という。 begin{eqnarray} s^2 = frac{n}{n-1} { (X_1-bar{X})^2 + (X_2-bar{X})^2 + cdots + (X_n-bar{X})^2 } end{eqnarray} 100個データがあって、10個データをとったときと、100個データをとったときの (E(S^2))の母分散とのズレは以下の通り。10個のとき(E(S^2))をそのまま計算してしまうと、 その値は母分散から10%もズレてしまう。100個にしても1%ずれる。 begin{eqnarray} E(S_{10}^2) &=& frac{9}{10}sigma^2 \\ E(S_{100}^2) &=& frac{99}{100}sigma^2 \\ end{eqnarray}