ブログ

Pythonを用いた推測統計のプログラミング 〜区間推定、仮説検定 etc〜

数式だけの解説ではわかりにくい場合もあると思われるので、統計学の手法や関連する概念をPythonのプログラミングで表現を行います。当記事では統計的推測(statistical inference)の理解にあたって、区間推定や仮説検定のPythonでの実装を取り扱いました。

・Pythonを用いた統計学のプログラミングまとめ
https://www.hello-statisticians.com/stat_program

基本事項の確認

区間推定

下記で詳しく取り扱った。
https://www.hello-statisticians.com/explain-terms-cat/flow_chart_stat1.html

仮説検定

推測統計でよく用いるSciPyのメソッド

SciPyの使用にあたっては、scipy.stats.xxx.oooのように使用するメソッドの指定を行うが、xxxに確率分布、oooに行う処理を指定する。

推測統計ではxxxに対して主に下記などの確率分布の指定を行う。

norm: 正規分布
t: $t$分布
chi2: $\chi^2$分布

同様に推測統計ではoooに対して、主に下記の処理の指定を行う。

cdf: 統計量から有意水準・$P$値の計算
ppf: 有意水準・$P$値から統計量の計算
pdf: 統計量から確率密度関数の計算

特にcdfppfは区間推定や検定で統計数値表の代わりに用いることができる。

以下、標準正規分布を例にcdfppfpdfの概要の確認を行う。

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

x = np.arange(-3.,3.01,0.01)
f_pdf = stats.norm.pdf(x)
F_cdf = stats.norm.cdf(x)

plt.plot(x,f_pdf,label="pdf")
plt.plot(x,F_cdf,label="cdf")

plt.legend()
plt.show()

・実行結果

上記のようにcdfpdfを用いることで累積分布関数と確率密度関数の描画を行うことができる。

次に、cdfppfの対応に関して確認を行う。

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

alpha = np.array([0.1, 0.05, 0.025, 0.01, 0.005])
x = np.arange(-3.,3.01,0.01)
F_cdf = stats.norm.cdf(x)

plt.plot(x,F_cdf,label="cdf")
for i in range(alpha.shape[0]):
    print("alpha: {}, statistic_upper: {:.2f}".format(alpha[i], stats.norm.ppf(1.-alpha[i])))
    x_upper = np.repeat(stats.norm.ppf(1.-alpha[i]), 100)
    y = np.linspace(0,1.-alpha[i],100)
    plt.plot(x_upper,y,label="alpha: {}".format(alpha[i]))

plt.legend(loc=2)
plt.show()

・実行結果

alpha: 0.1, statistic_upper: 1.28
alpha: 0.05, statistic_upper: 1.64
alpha: 0.025, statistic_upper: 1.96
alpha: 0.01, statistic_upper: 2.33
alpha: 0.005, statistic_upper: 2.58

上記よりcdfppfは逆関数であると考えることができる。それぞれの対応がわかりやすいように、下記では上図の編集を行った。

統計数値表とSciPyを用いた代用

scipy.statsscipy.stats.normscipy.stats.tscipy.stats.chi2を用いることで、統計的推測を行うことができる。特に区間推定や仮説検定でよく用いられる統計数値表が、累積分布関数に対応するscipy.stats.norm.cdfなどを活用することで作成できることは抑えておくと良い。

以下ではscipy.stats.norm.cdfscipy.stats.norm.ppfなどを用いることで統計数値表を作成し概要の理解を行うことに加えて、SciPyで代用する際の用い方について確認を行う。

標準正規分布

統計数値表の作成

上側確率が$100 \alpha’$%となる点を$z_{\alpha=\alpha’}$とおくとき、$z_{\alpha=\alpha’}=0$から$0.01$刻みで該当する$\alpha’$の値を表形式で表したものが統計数値表である。以下、$500$個の値を$50$行$10$列に並べることで標準正規分布の統計数値表の作成を行う。

import numpy as np
from scipy import stats

z = np.arange(0., 5., 0.01).reshape([50,10])
table_norm = 1. - stats.norm.cdf(z)

print(table_norm)

・実行結果

array([[  5.00000000e-01,   4.96010644e-01,   4.92021686e-01,
          4.88033527e-01,   4.84046563e-01,   4.80061194e-01,
          4.76077817e-01,   4.72096830e-01,   4.68118628e-01,
          4.64143607e-01],
       [  4.60172163e-01,   4.56204687e-01,   4.52241574e-01,
          4.48283213e-01,   4.44329995e-01,   4.40382308e-01,
          4.36440537e-01,   4.32505068e-01,   4.28576284e-01,
          4.24654565e-01],
.....
[  4.79183277e-07,   4.55381965e-07,   4.32721062e-07,
          4.11148084e-07,   3.90612854e-07,   3.71067408e-07,
          3.52465898e-07,   3.34764508e-07,   3.17921366e-07,
          3.01896462e-07]])

処理の理解にあたっては、上側確率であることから1. - stats.norm.cdf(z)のように$1$から引いたと考えると良い。

SciPyを用いた統計数値表の代用

・有意水準$\alpha$から標準化変量$Z$を導出

有意水準$\alpha$から標準化変量$Z$を導出する場合はscipy.stats.norm.ppfを用いればよい。

import numpy as np
from scipy import stats

alpha = np.array([0.1, 0.05, 0.025, 0.01, 0.005])
Z_Limit = stats.norm.ppf(1.-alpha)

print(Z_Limit)

・実行結果

> print(Z_Limit)
array([ 1.28155157,  1.64485363,  1.95996398,  2.32634787,  2.5758293 ])

標準化変量$Z$から有意水準$\alpha$を導出

標準化変量$Z$から有意水準$\alpha$を導出する場合は、統計数値表の作成と同様にscipy.stats.norm.cdfを用いればよい。

import numpy as np
from scipy import stats

Z = np.array([0., 0.5, 1., 1.5, 2., 2.5, 3.])
alpha = 1.-stats.norm.cdf(Z)

print(alpha)

・実行結果

> print(alpha)
[ 0.5         0.30853754  0.15865525  0.0668072   0.02275013  0.00620967
  0.0013499 ]

$t$分布

統計数値表の作成

自由度が$\nu$の$t$分布$t(\nu)$に関して、上側確率が$100 \alpha’$%となる点を$t_{\alpha=\alpha’}(\nu)$とおくとき、$\alpha’=\{0.25, 0.20, 0.15, 0.10, 0.05, 0.025, 0.01, 0.005, 0.001\}$に対して、$\nu$を変化させたときに$t_{\alpha=\alpha’}(\nu)$の値をまとめたものが$t$分布の統計数値表である。以下、それぞれの有意水準と自由度に対して$t$分布の統計数値表の作成を行う。

import numpy as np
from scipy import stats

alpha = np.array([0.25, 0.20, 0.15, 0.10, 0.05, 0.025, 0.01, 0.005, 0.0005])
nu = np.arange(1,56,1)
nu[50:] = np.array([60., 80., 120., 240., 100000])

table_t = np.zeros([nu.shape[0],alpha.shape[0]])
for i in range(nu.shape[0]):
    table_t[i,:] = stats.t.ppf(1.-alpha,nu[i])

print(table_t)

・実行結果

[[   1.            1.37638192    1.96261051    3.07768354    6.31375151
    12.70620474   31.82051595   63.65674116  636.61924943]
 [   0.81649658    1.06066017    1.38620656    1.88561808    2.91998558
     4.30265273    6.96455672    9.9248432    31.59905458]
 [   0.76489233    0.97847231    1.2497781     1.63774436    2.35336343
     3.18244631    4.54070286    5.8409093    12.92397864]
.....
 [   0.67551336    0.84312147    1.0386776     1.28508893    1.65122739
     1.96989764    2.34198547    2.59646918    3.33152484]
 [   0.6744922     0.84162483    1.03643876    1.28156003    1.64486886
     1.95998771    2.32638517    2.57587847    3.29062403]]

処理の理解にあたっては、上側確率であることから1. - alphaのように$1$から引いたと考えると良い。

SciPyを用いた統計数値表の代用

・有意水準$\alpha$からパーセント点$t_{\alpha=\alpha’}(\nu)$を導出

有意水準$\alpha$からパーセント点$t_{\alpha=\alpha’}(\nu)$を導出する場合は統計数値表の作成と同様にscipy.stats.t.ppfを用いればよい。

import numpy as np
from scipy import stats

alpha = np.array([0.1, 0.05, 0.025, 0.01, 0.005])
Z_Limit = stats.norm.ppf(1.-alpha)
t_Limit_20 = stats.t.ppf(1.-alpha, 20)
t_Limit_100 = stats.t.ppf(1.-alpha, 100)

print(Z_Limit)
print(t_Limit_20)
print(t_Limit_100)

・実行結果

> print(Z_Limit)
[ 1.28155157  1.64485363  1.95996398  2.32634787  2.5758293 ]
> print(t_Limit_20)
[ 1.32534071  1.72471824  2.08596345  2.527977    2.84533971]
> print(t_Limit_100)
[ 1.29007476  1.66023433  1.98397152  2.36421737  2.62589052]

上記より、$t$分布の棄却点は標準正規分布の棄却点より絶対値が大きくなることが確認できる。また、t_Limit_20t_Limit_100の比較により、自由度が大きくなるにつれて$t$分布が標準正規分布に近づくことも確認できる。

・パーセント点$t_{\alpha=\alpha’}(\nu)$から有意水準$\alpha$を導出

統計量$t$から有意水準$\alpha$を導出する場合は、scipy.stats.t.cdfを用いればよい。

import numpy as np
from scipy import stats

t = np.array([0., 0.5, 1., 1.5, 2., 2.5, 3.])
alpha_20 = 1.-stats.t.cdf(t,20)
alpha_100 = 1.-stats.t.cdf(t,100)

print(alpha_20)
print(alpha_100)

・実行結果

> print(alpha_20)
[ 0.5         0.31126592  0.16462829  0.07461789  0.02963277  0.01061677
  0.00353795]
> print(alpha_100)
[ 0.5         0.30908678  0.15986208  0.06838253  0.02410609  0.00702289
  0.00170396]

$\chi^2$分布

統計数値表の作成

自由度が$\nu$の$\chi^2$分布$\chi^2(\nu)$に関して、上側確率が$100 \alpha’$%となる点を$\chi^2_{\alpha=\alpha’}(\nu)$とおくとき、$\alpha’=\{0.995, 0.990, 0.975, 0.950, 0.900, 0.800, …, 0.100, 0.050, 0.025, 0.010, 0.005, 0.001\}$に対して、$\nu$を変化させたときに$t_{\alpha=\alpha’}(\nu)$の値をまとめたものが「基礎統計学Ⅰ」における$\chi^2$分布の統計数値表である。以下、それぞれの有意水準と自由度に対して$\chi^2$分布の統計数値表の作成を行う。

import numpy as np
from scipy import stats

alpha = np.array([0.995, 0.990, 0.975, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.025, 0.01, 0.005, 0.001])
nu = np.arange(1,62,1)
nu[50:] = np.array([60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 240.])

table_chi2 = np.zeros([nu.shape[0],alpha.shape[0]])
for i in range(nu.shape[0]):
    table_chi2[i,:] = stats.chi2.ppf(1.-alpha,nu[i])

print(table_chi2)

・実行結果

[[  3.92704222e-05   1.57087858e-04   9.82069117e-04 ...,   6.63489660e+00
    7.87943858e+00   1.08275662e+01]
 [  1.00250836e-02   2.01006717e-02   5.06356160e-02 ...,   9.21034037e+00
    1.05966347e+01   1.38155106e+01]
 [  7.17217746e-02   1.14831802e-01   2.15795283e-01 ...,   1.13448667e+01
    1.28381565e+01   1.62662362e+01]
.....
 [  1.52240992e+02   1.56431966e+02   1.62727983e+02 ...,   2.49445123e+02
    2.55264155e+02   2.67540528e+02]
 [  1.87324255e+02   1.91989896e+02   1.98983851e+02 ...,   2.93888101e+02
    3.00182237e+02   3.13436899e+02]]

SciPyを用いた統計数値表の代用

・有意水準$\alpha$からパーセント点$\chi^2_{\alpha=\alpha’}(\nu)$を導出

有意水準$\alpha$からパーセント点$\chi^2_{\alpha=\alpha’}(\nu)$を導出する場合は統計数値表の作成と同様にscipy.stats.chi2.ppfを用いればよい。

import numpy as np
from scipy import stats

alpha = np.array([0.975, 0.95, 0.1, 0.05, 0.025])
chi2_Limit_1 = stats.chi2.ppf(1.-alpha, 1)
chi2_Limit_10 = stats.chi2.ppf(1.-alpha, 10)
chi2_Limit_100 = stats.chi2.ppf(1.-alpha, 100)

print(chi2_Limit_1)
print(chi2_Limit_10)
print(chi2_Limit_100)

・実行結果

print(chi2_Limit_1)
[  9.82069117e-04   3.93214000e-03   2.70554345e+00   3.84145882e+00
   5.02388619e+00]
print(chi2_Limit_10)
[  3.24697278   3.94029914  15.98717917  18.30703805  20.48317735]
print(chi2_Limit_100)
[  74.22192747   77.92946517  118.49800381  124.3421134   129.56119719]

・パーセント点$\chi^2_{\alpha=\alpha’}(\nu)$から有意水準$\alpha$を導出

統計量$\chi^2$から有意水準$\alpha$を導出する場合は、scipy.stats.chi2.cdfを用いればよい。

import numpy as np
from scipy import stats

chi2 = np.array([2., 5., 10., 15., 20., 50., 100.])

alpha_1 = 1.-stats.chi2.cdf(chi2,1)
alpha_5 = 1.-stats.chi2.cdf(chi2,5)
alpha_10 = 1.-stats.chi2.cdf(chi2,10)
alpha_20 = 1.-stats.chi2.cdf(chi2,20)
alpha_100 = 1.-stats.chi2.cdf(chi2,100)

print(alpha_1)
print(alpha_5)
print(alpha_10)
print(alpha_20)
print(alpha_100)

・実行結果

> print(alpha_1)
[  1.57299207e-01   2.53473187e-02   1.56540226e-03   1.07511177e-04
   7.74421643e-06   1.53743684e-12   0.00000000e+00]
> print(alpha_5)
[  8.49145036e-01   4.15880187e-01   7.52352461e-02   1.03623379e-02
   1.24973056e-03   1.38579737e-09   0.00000000e+00]
> print(alpha_10)
[  9.96340153e-01   8.91178019e-01   4.40493285e-01   1.32061856e-01
   2.92526881e-02   2.66908343e-07   0.00000000e+00]
> print(alpha_20)
[  9.99999889e-01   9.99722648e-01   9.68171943e-01   7.76407613e-01
   4.57929714e-01   2.21476638e-04   1.25965904e-12]
> print(alpha_100)
[ 1.          1.          1.          1.          1.          0.99999305
  0.48119168]

区間推定・仮説検定

母分散既知の場合の母平均

区間推定

$$
\large
\begin{align}
\bar{x}-1.96\frac{\sigma}{\sqrt{n}} \leq &\mu \leq \bar{x}+1.96\frac{\sigma}{\sqrt{n}}
\end{align}
$$
母平均の$95$%区間推定は上記のように表される。以下、上記を元に具体的な例に対して計算を行う。

import numpy as np
from scipy import stats

sigma = 1.
x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
ave_x = np.mean(x)

Lower_mu = ave_x + stats.norm.ppf(0.025)*sigma/np.sqrt(x.shape[0])
Upper_mu = ave_x + stats.norm.ppf(0.975)*sigma/np.sqrt(x.shape[0])

print("Estimated mu-interval: [{:.3f}, {:.3f}]".format(Lower_mu, Upper_mu))

・実行結果

> print("Estimated mu-interval: [{:.3f}, {:.3f}]".format(Lower_mu, Upper_mu))
Estimated mu-interval: [0.203, 1.957]

以下、$\sigma$の値を変化させた際の区間推定の結果の確認を行う。

import numpy as np
from scipy import stats

sigma = np.array([0.1, 0.25, 0.5, 1., 2.])
x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
ave_x = np.mean(x)

Lower_mu = ave_x + stats.norm.ppf(0.025)*sigma/np.sqrt(x.shape[0])
Upper_mu = ave_x + stats.norm.ppf(0.975)*sigma/np.sqrt(x.shape[0])

for i in range(sigma.shape[0]):
    print("Estimated mu-interval: [{:.3f}, {:.3f}]".format(Lower_mu[i], Upper_mu[i]))

・実行結果

Estimated mu-interval: [0.992, 1.168]
Estimated mu-interval: [0.861, 1.299]
Estimated mu-interval: [0.642, 1.518]
Estimated mu-interval: [0.203, 1.957]
Estimated mu-interval: [-0.673, 2.833]

Lower_muUpper_muの計算にあたって、それぞれ0.0250.975を引数に設定したが、これは上側確率を$\alpha’$とおく際の$1-\alpha’$に対応し、このように考えると少々わかりにくい。そのため、標準正規分布の累積分布関数の値から変数の値を計算すると考える方が良い。

これらの違いは元々が「統計数値表」に基づいて考えていた一方で、SciPyなどで実装を行う際には「数値表」が必要ないことからより関数的な表し方に変わったのではないかと考えられる。

SciPyでの表し方がより自然に思われるので、表よりもSciPyなどを用いた形式で理解する方がわかりやすいかもしれない。

仮説検定① 〜norm.ppfと区間推定の活用〜

前項のxave_xに対し、帰無仮説を$H_0: \mu=0$と設定し、上側確率$\alpha’=\{0.05, 0.025, 0.005\}$を用いた仮説検定を行う。
$$
\large
\begin{align}
z_{\alpha=1-\alpha’} \leq &\frac{\bar{x}-\mu}{\sigma/\sqrt{n}} \leq z_{\alpha=\alpha’} \\
-z_{\alpha=\alpha’} \leq &\frac{\sqrt{n}}{\sigma}\bar{x} \leq z_{\alpha=\alpha’}
\end{align}
$$
ここでの検定統計量の$\displaystyle Z = \frac{\sqrt{n}}{\sigma}\bar{x}$に関して上記が成立しなければ$H_0$を棄却し、対立仮説の$H_1: \mu \neq 0$を採択する。

import numpy as np
from scipy import stats

alpha = np.array([0.05, 0.025, 0.005])
sigma = 1.
x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
Z_observed = np.sqrt(x.shape[0])*np.mean(x)/sigma

Lower_Z = -stats.norm.ppf(1.-alpha)*sigma
Upper_Z = stats.norm.ppf(1.-alpha)*sigma

for i in range(alpha.shape[0]):
    if Z_observed > Upper_Z[i] or Z_observed < Lower_Z[i]:
        print("Z: {:.3f}. If alpha: {} -> Upper_Z: {:.2f}, reject H_0: mu=0. ".format(Z_observed, alpha[i], Upper_Z[i]))
    else:
        print("Z: {:.3f}. If alpha: {} -> Upper_Z: {:.2f}, accept H_0: mu=0. ".format(Z_observed, alpha[i], Upper_Z[i]))

・実行結果

Z: 2.415. If alpha: 0.05 -> Upper_Z: 1.64, reject H_0: mu=0. 
Z: 2.415. If alpha: 0.025 -> Upper_Z: 1.96, reject H_0: mu=0. 
Z: 2.415. If alpha: 0.005 -> Upper_Z: 2.58, accept H_0: mu=0.

上記より、上側確率を小さくするにつれて棄却されにくいことが確認できる。

仮説検定② 〜norm.cdfと$P$値の活用〜

前項の「仮説検定①」では区間推定の結果と$Z$の値を元に検定を行なったが、$Z$の値に対してnorm.cdfを用いることで$P$値を計算し、検定を行うこともできる。

import numpy as np
from scipy import stats

alpha = np.array([0.05, 0.025, 0.005])
sigma = 1.
x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
Z_observed = np.sqrt(x.shape[0])*np.mean(x)/sigma

P_value = 1.-stats.norm.cdf(Z_observed)

for i in range(alpha.shape[0]):
    if P_value > 1.-alpha[i] or P_value < alpha[i]:
        print("Z: {:.3f}, P_value: {:.3f}. If alpha: {}, reject H_0: mu=0. ".format(Z_observed, P_value, alpha[i]))
    else:
        print("Z: {:.3f}, P_value: {:.3f}. If alpha: {}, accept H_0: mu=0. ".format(Z_observed, P_value, alpha[i]))

・実行結果

Z: 2.415, P_value: 0.008. If alpha: 0.05, reject H_0: mu=0. 
Z: 2.415, P_value: 0.008. If alpha: 0.025, reject H_0: mu=0. 
Z: 2.415, P_value: 0.008. If alpha: 0.005, accept H_0: mu=0. 

仮説検定は②で取り扱った「$P$値の活用」を用いる方が「ノンパラメトリック法」と同様に取り扱うことができるなど、わかりやすい。また、この際にppfcdfの図を元に下記のように考えることもできる。

母分散未知の場合の母平均

区間推定

$$
\large
\begin{align}
\bar{x}-t_{\alpha=0.025}(n-1) \frac{s}{\sqrt{n}} \leq &\mu \leq \bar{x}+t_{\alpha=0.025}(n-1) \frac{s}{\sqrt{n}} \\
s = \frac{1}{n-1} & \sum_{i=1}^{n} (x_i-\bar{x})^2
\end{align}
$$
母分散が未知の際の母平均の$95$%区間推定は上記のように表される。以下、上記を元に具体的な例に対して計算を行う。

import numpy as np
from scipy import stats

x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
ave_x = np.mean(x)
s = np.sqrt(np.sum((x-ave_x)**2/(x.shape[0]-1)))

Lower_mu = ave_x - stats.t.ppf(1-0.025, x.shape[0]-1)*s/np.sqrt(x.shape[0])
Upper_mu = ave_x + stats.t.ppf(1-0.025, x.shape[0]-1)*s/np.sqrt(x.shape[0])

print("s: {:.2f}".format(s))
print("Estimated mu-interval: [{:.3f}, {:.3f}]".format(Lower_mu, Upper_mu))

・実行結果

> print("s: {:.2f}".format(s))
s: 0.30
> print("Estimated mu-interval: [{:.3f}, {:.3f}]".format(Lower_mu, Upper_mu))
Estimated mu-interval: [0.703, 1.457]

母分散既知の場合と同様の事例を用いたが、$s=0.3$の周辺の$\sigma=0.25$の際の区間が[0.861, 1.299]、$\sigma=0.5$の際の区間が[0.642, 1.518]であったので、概ね同様の結果が得られていることが確認できる。

仮説検定

前項のxave_xに対し、帰無仮説を$H_0: \mu_0=0.95$と設定し、上側確率$\alpha’=\{0.05, 0.025, 0.005\}$を用いた仮説検定を行う。
$$
\large
\begin{align}
t_{\alpha=1-\alpha’}(n-1) \leq &\frac{\bar{x}-\mu_0}{s/\sqrt{n}} \leq t_{\alpha=\alpha’}(n-1) \\
-t_{\alpha=\alpha’}(n-1) \leq &\frac{\sqrt{n}}{s} (\bar{x}-0.95) \leq t_{\alpha=\alpha’}(n-1)
\end{align}
$$
ここでの検定統計量の$\displaystyle t = \frac{\sqrt{n}}{s} (\bar{x}-\mu_0)$に関して上記が成立しなければ$H_0$を棄却し、対立仮説の$H_1: \mu_0 \neq 0.95$を採択する。

import numpy as np
from scipy import stats

mu_0 = 0.95
alpha = np.array([0.05, 0.025, 0.005]) 

x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
ave_x = np.mean(x)
s = np.sqrt(np.sum((x-ave_x)**2/(x.shape[0]-1)))
t_observed = np.sqrt(x.shape[0])*(ave_x-mu_0)/s

Lower_t = -stats.t.ppf(1.-alpha, x.shape[0])*s
Upper_t = stats.t.ppf(1.-alpha, x.shape[0])*s

for i in range(alpha.shape[0]):
    if t_observed > Upper_t[i] or t_observed < Lower_t[i]:
        print("t: {:.3f}. If alpha: {}, reject H_0: mu=0. ".format(t_observed, alpha[i]))
    else:
        print("t: {:.3f}. If alpha: {}, accept H_0: mu=0. ".format(t_observed, alpha[i]))

print(Lower_t)
print(Upper_t)

・実行結果

t: 0.958. If alpha: 0.05, reject H_0: mu=0. 
t: 0.958. If alpha: 0.025, reject H_0: mu=0. 
t: 0.958. If alpha: 0.005, accept H_0: mu=0. 

> print(Lower_t)
[-0.61119443 -0.77969608 -1.22300952]
> print(Upper_t)
[ 0.61119443  0.77969608  1.22300952]

上記より、上側確率を小さくするにつれて棄却されにくいことが確認できる。

母分散

区間推定

$$
\large
\begin{align}
\frac{(n-1)s^2}{\chi_{\alpha=0.025}^2(n-1)} \leq & \sigma^2 \leq \frac{(n-1)s^2}{\chi_{\alpha=0.975}^2(n-1)} \\
s = \frac{1}{n-1} & \sum_{i=1}^{n} (x_i-\bar{x})^2
\end{align}
$$
母分散の$95$%区間推定は上記のように表される。以下、上記を元に具体的な例に対して計算を行う。

import numpy as np
from scipy import stats

x = np.array([1.2, 1.5, 0.9, 0.7, 1.1])
ave_x = np.mean(x)
s = np.sqrt(np.sum((x-ave_x)**2/(x.shape[0]-1)))

Lower_sigma2 = (x.shape[0]-1)*s**2/stats.chi2.ppf(1.-0.025, x.shape[0]-1)
Upper_sigma2 = (x.shape[0]-1)*s**2/stats.chi2.ppf(1.-0.975, x.shape[0]-1)

print("s^2: {:.2f}".format(s**2))
print("Estimated sigma2s-interval: [{:.3f}, {:.3f}]".format(Lower_sigma2, Upper_sigma2))

・実行結果

> print("s^2: {:.2f}".format(s**2))
s^2: 0.09
> print("Estimated sigma2s-interval: [{:.3f}, {:.3f}]".format(Lower_sigma2, Upper_sigma2))
Estimated sigma2s-interval: [0.033, 0.760]

以下、サンプルを増やした際に$95$%区間はどのように変化するかに関して確認を行う。

import numpy as np
from scipy import stats

x = np.array([1.2, 1.5, 0.9, 0.7, 1.1, 0.5, 1.0, 1.3, 0.9, 1.2])
ave_x = np.mean(x)
s = np.sqrt(np.sum((x-ave_x)**2/(x.shape[0]-1)))

Lower_sigma2 = (x.shape[0]-1)*s**2/stats.chi2.ppf(1.-0.025, x.shape[0]-1)
Upper_sigma2 = (x.shape[0]-1)*s**2/stats.chi2.ppf(1.-0.975, x.shape[0]-1)

print("s^2: {:.2f}".format(s**2))
print("Estimated sigma2s-interval: [{:.3f}, {:.3f}]".format(Lower_sigma2, Upper_sigma2))

・実行結果

> print("s^2: {:.2f}".format(s**2))
s^2: 0.09
> print("Estimated sigma2s-interval: [{:.3f}, {:.3f}]".format(Lower_sigma2, Upper_sigma2))
Estimated sigma2s-interval: [0.041, 0.289]

上記よりサンプルを増やすと$\sigma^2$の$95$%区間が小さくなることが確認できる。

仮説検定

前項のxに対し、帰無仮説を$H_0: \sigma_0^2=0.3$と設定し、上側確率$\alpha’=\{0.05, 0.025, 0.005\}$を用いた仮説検定を行う。
$$
\large
\begin{align}
\chi_{\alpha=0.975}(n-1) \leq & \frac{(n-1)s^2}{\sigma_0^2} \leq \chi_{\alpha=0.025}(n-1) \\
s = \frac{1}{n-1} & \sum_{i=1}^{n} (x_i-\bar{x})^2
\end{align}
$$
ここでの検定統計量の$\displaystyle \chi^2 = \frac{(n-1)s^2}{\sigma_0^2}$に関して上記が成立しなければ$H_0$を棄却し、対立仮説の$H_1: \sigma^2 \neq 0.3$を採択する。

import numpy as np
from scipy import stats

sigma2_0 = 0.3
alpha = np.array([0.05, 0.025, 0.005]) 

x = np.array([1.2, 1.5, 0.9, 0.7, 1.1, 0.5, 1.0, 1.3, 0.9, 1.2])
ave_x = np.mean(x)
s = np.sqrt(np.sum((x-ave_x)**2/(x.shape[0]-1)))

chi2_observed = (x.shape[0]-1)*s**2/sigma2_0

Lower_chi2 = stats.chi2.ppf(alpha, x.shape[0]-1)
Upper_chi2 = stats.chi2.ppf(1.-alpha, x.shape[0]-1)

for i in range(alpha.shape[0]):
    if chi2_observed > Upper_chi2[i] or chi2_observed < Lower_chi2[i]:
        print("chi2: {:.3f}. If alpha: {}, reject H_0: sigma^2=0.3. ".format(chi2_observed, alpha[i]))
    else:
        print("chi2: {:.3f}. If alpha: {}, accept H_0: sigma^2=0.3. ".format(chi2_observed, alpha[i]))

print(Lower_chi2)
print(Upper_chi2)

・実行結果

chi2: 2.603. If alpha: 0.05, reject H_0: sigma^2=0.3. 
chi2: 2.603. If alpha: 0.025, reject H_0: sigma^2=0.3. 
chi2: 2.603. If alpha: 0.005, accept H_0: sigma^2=0.3. 

> print(Lower_t)
[ 3.32511284  2.7003895   1.7349329 ]
> print(Upper_t)
[ 16.9189776   19.0227678   23.58935078]

上記より、上側確率を小さくするにつれて棄却されにくいことが確認できる。

母平均の差

区間推定

仮説検定

母分散の比

区間推定

仮説検定

事前分布・ベイズとMAP推定・予測分布|問題演習で理解する統計学【20】

最尤法は観測値のみを元に推定を行うが、トピックによっては事前知識がわかっている場合や、サンプル数が少ない場合に最尤法が極端な結果を示す場合の補正にあたって、事前分布に基づくベイズ法は役に立つ。事前分布やMAP推定、予測分布が理解できるように演習を取り扱った。
・現代数理統計学 Ch.$14$ 「ベイズ法」の章末演習の解答例
https://www.hello-statisticians.com/explain-books-cat/math_stat_practice_ch14.html

・標準演習$100$選
https://www.hello-statisticians.com/practice_100

基本問題

ベイズの定理と事前分布・事後分布

・問題
パラメータ推定におけるベイズの定理の導入は「ベイズ統計」や「ベイズ推定」など様々な表現があるが、数式に基づく基本的な流れを抑えておけば基本的に十分である。この問題では以下、最尤推定に対し事前分布を元にベイズの定理を考える一連の流れについて取り扱う。下記の問いにそれぞれ答えよ。

i) 事象$A,B$がある際に、$A \cap B$は条件付き確率を元に下記のように表せる。
$$
\begin{align}
P(A \cap B) = P(A|B)P(B) = P(B|A)P(A) \quad (1)
\end{align}
$$

上記の式が妥当であることを直感的に解釈せよ。

ⅱ) $(1)$式より、$P(A|B)$を$P(B), P(B|A), P(A)$を用いて表せ。

ⅲ) ⅱ)の式がベイズの定理の式に一致する。以下ではベイズの定理を最尤法に適用する流れに関して二項分布を例に確認を行う。$X \sim \mathrm{Bin}(n,p)$の確率関数$p(x)=P(x|p)$は下記のように表せる。
$$
\begin{align}
p(x) = {}_{n} C_{x} p^{x} (1-p)^{n-x} \quad (2)
\end{align}
$$

ここで$n$回のベルヌーイ試行に基づいて観測値$X=x$が得られた際に$p$の推定を行うことを考える。最尤推定は尤度を$L(p)=P(x|p)$とおき、$L(p)$が最大になる$p$を導出する手法である。
$(2)$式を元に対数を考えることで、$L(p)$が最大値を取る$p$を$x, n$を用いて表せ。

iv) $n=2$のとき、ⅲ)の式を用いて$x=0,1,2$に関して$p$の推定を行え。

v) $n=2$のようにサンプルが少ない場合、iv)の推定結果は$x$の値次第で大きく変動する。この変動の緩和にあたって下記のように事前分布$P(p)$を導入することを考える。
$$
\begin{align}
P(p) = 6p(1-p), \quad 0 \leq p \leq 1
\end{align}
$$

ここで$P(x)=c(x)$のように表すとき、ベイズの定理に基づいて事後分布$P(p|x)$を表せ。

vi) v)で得られた事後分布$P(p|x)$が最大になる$p$を$n, x$を用いて表せ。
vⅱ) $n=2$のときvi)の式を用いて$x=0,1,2$に関しての$p$の推定を行い、事前分布や事後分布を導入した意義に関して論じよ。

・解答
i)
$(1)$式は「$2$つの事象が同時に発生する確率」と「どちらか一方が発生し、その後にもう一方が発生する確率」のどちらで考えても式が同様であると直感的に解釈できる。

ⅱ)
下記のように表せる。
$$
\large
\begin{align}
P(A|B) = \frac{P(B|A)P(A)}{P(B)}
\end{align}
$$

ⅲ)
対数尤度$\log{L(p)}=\log{P(x|p)}$は下記のように考えることができる。
$$
\large
\begin{align}
\log{L(p)} &= \log{P(x|p)} \\
&= \log{({}_{n} C_{x} p^{x} (1-p)^{n-x})} \\
&= x \log{p} + (n-x) \log{(1-p)} + \log{{}_{n} C_{x}}
\end{align}
$$

iv)
ⅲ)の結果に$n=2$と$x=0,1,2$を代入することで$\hat{p}=0, \hat{p}=0.5, \hat{p}=1$がそれぞれ得られる。

v)
ベイズの定理に基づいて事後分布$P(p|x)$は下記のように表せる。
$$
\large
\begin{align}
P(p|x) &= \frac{P(x|p)P(p)}{P(x)} \\
&= \frac{6 {}_{n} C_{x}}{c(x)} p^{x+1} (1-p)^{n-x+1}
\end{align}
$$

vi)
事後分布の対数$\log{P(p|x)}$は下記のように表せる。
$$
\large
\begin{align}
\log{P(p|x)} &= \log{ \left[ \frac{6 {}_{n} C_{x}}{c(x)} p^{x+1} (1-p)^{n-x+1} \right] } \\
&= (x+1) \log{p} + (n-x+1) \log{(1-p)} + \log{ \left[ \frac{6 {}_{n} C_{x}}{c(x)} \right] }
\end{align}
$$

上記を$p$に関して微分すると下記が得られる。
$$
\large
\begin{align}
\frac{\partial \log{P(p|x)}}{\partial p} &= \frac{x+1}{p} – \frac{n-x+1}{1-p} \\
&= \frac{(x+1)(1-p) – (n-x+1)p}{p(1-p)} \\
&= \frac{x+1-(n+2)p}{p(1-p)}
\end{align}
$$

上記は$p$に関して単調増加であるので、$\log{L(p)}$や$L(p)$は$\displaystyle p = \frac{x+1}{n+2}$のとき最大値を取る。

vⅱ)
vi)の結果に$n=2$と$x=0,1,2$を代入することで$\hat{p}=0.25, \hat{p}=0.5, \hat{p}=0.75$がそれぞれ得られる。ここでiv)では$\hat{p}=0, \hat{p}=0.5, \hat{p}=1$が得られたことを鑑みると、事前分布を導入することで$0.5$に近い推定結果が得られたことが確認できる。
このように事前分布を導入することでサンプル数が少ない場合のサンプルのばらつきに対して頑健である結果が得られることが多い。

・解説
この問題では「二項分布の共役事前分布」を題材にベイズの定理と事前分布・事後分布に関して確認を行いました。詳しい内容は「二項分布の共役事前分布と事後分布の解釈」で取り扱いました。
vⅱ)では事前分布を導入した意義に関して論じましたが、$p=0.5$に近い値の場合は妥当である一方で$p=0.01$や$p=0.99$の場合には妥当でない結果になり得ることに注意が必要です。手法選択の際には一見ベイズ手法が万能であるように見えがちですが、「ベイズ推定は最尤推定よりも必ず良い」ではなく、「妥当な事前分布が設定できる場合はサンプルが少ない際に有用である」のように理解しておくのが良いように思います。

二項分布の共役事前分布と事後分布の解釈

・問題
共役事前分布(conjugate prior distribution)を用いることで、事後分布の導出をシンプルに行うことができる。共役事前分布はフォーマルな表記だと難しい印象があると思われるので、以下では二項分布を例に共役事前分布について確認し、事後分布の解釈を行う。

ベイズ法ではパラメータの推定にあたって、下記のように事後分布の$P(\theta|x)$が事前分布の$P(\theta)$と尤度の$p(x|\theta)$の積に比例することを利用する。
$$
\large
\begin{align}
P(\theta|x) \propto P(\theta)p(x|\theta) \quad (1)
\end{align}
$$

上記の式は多くの書籍などで出てくるが、抽象的な表記より具体的な確率分布を例に考える方がわかりやすい。以下では二項分布に関して$(1)$式を表すことを考える。

ここまでの内容を元に下記の問いに答えよ。
i) 二項分布におけるパラメータは事象が起こる確率であるので、以下では$\theta=p$のように、パラメータを$\theta$ではなく$p$を用いて表す。確率$p$の試行を$n$回繰り返すとき、事象が起こる回数を$x$、この確率関数を$p(x|p,n)$のように表すとき、確率関数$p(x|p,n)$の式を表せ。
ⅱ) ベータ分布$Be(a,b)$の確率密度関数を答えよ。
ⅲ) i)で取り扱った確率関数$p(x|p,n)$は尤度と考えることができる。ここで$p(x|p,n)$を$x$ではなく$p$の関数と見るときベータ分布と同様の形を持つことをⅱ)の結果と見比べることで示せ。
iv) $p$の事前分布の$P(p)$がベータ分布$Be(a,b)$であるとき、事前分布$P(p)$の関数を記載せよ。
v) $p(x|p,n)P(p)$を計算せよ。
vi) 事後分布$P(p|x)$が$P(p|x) \propto p(x|p,n)P(p)$のように表されることを元に、事後分布$P(p|x)$がどのような分布になるか答えよ。
vⅱ) v)、vi)の結果を確認することでベータ分布$Be(a,b)$の$a,b$のこの事例での意味を考察せよ。

・解答
i)
確率関数$p(x|p,n)$は下記のように表すことができる。
$$
\large
\begin{align}
p(x|p,n) = {}_{n} C_{x} p^{x} (1-p)^{n-x}
\end{align}
$$

ⅱ)
ベータ分布$Be(a,b)$の確率密度関数を$f(x|a,b)$とおくとき、$f(x|a,b)$は下記のように表される。
$$
\large
\begin{align}
f(x|a,b) = \frac{1}{B(a,b)} x^{a-1} (1-x)^{b-1}
\end{align}
$$
ここで上記における$B(a,b)$はベータ関数を表す。

ⅲ)
i)とⅱ)の結果を見比べることで、どちらも変数を$x$と見た際に、$x^{a-1} (1-x)^{b-1}$の形状の関数であることを読み取ることができる。

iv)
事前分布$P(p)$はⅱ)の結果の$x$を$p$で置き換えることで得られるので、下記のように表すことができる。
$$
\large
\begin{align}
P(p) = \frac{1}{B(a,b)} p^{a-1} (1-p)^{b-1}
\end{align}
$$

v)
$p(x|p,n)P(p)$は下記のように計算できる。
$$
\large
\begin{align}
p(x|p,n)P(p) &= {}_{n} C_{x} p^{x} (1-p)^{n-x} \times \frac{1}{B(a,b)} p^{a-1} (1-p)^{b-1} \\
&= \frac{{}_{n} C_{x}}{B(a,b)} p^{a+x-1} (1-p)^{b+n-x-1}
\end{align}
$$

vi)
v)の結果より、下記が成立する。
$$
\large
\begin{align}
P(p|x) & \propto p(x|p,n)P(p) \\
&= {}_{n} C_{x} p^{x} (1-p)^{n-x} \times \frac{1}{B(a,b)} p^{a-1} (1-p)^{b-1} \\
& \propto p^{a+x-1} (1-p)^{b+n-x-1}
\end{align}
$$
よって事後分布$P(p|x)$はベータ分布$Be(a+x,b+n-x)$に一致する。

vⅱ)
二項分布における$x, n-x$はコイン投げにおける表と裏の回数に一致する。事前分布における$(a,b)$が事後分布において$(a+x,b+n-x)$になることは、実際の実験を行う前に事前に表が$a$回、裏が$b$回出たと考えることに一致する。

・解説
vⅱ)より、事前分布のベータ分布におけるパラメータの$a, b$に事前知識を組み込むことが可能になります。このようにすることで、試行回数が$5$回や$10$回などの少ない際に事前知識による補正ができると抑えておくと良いと思います。

ポアソン分布の共役事前分布と事後分布の解釈

・問題
ポアソン分布の共役事前分布にはガンマ分布を考えることができるので、ガンマ分布の解釈と合わせて理解すれば良い。以下、導出の流れとガンマ分布の解釈に関して確認を行う。

ポアソン分布の確率関数を$f(x|\lambda)$とおくと、$f(x|\lambda)$は下記のように表される。
$$
\begin{align}
f(x|\lambda) = \frac{\lambda^{x} e^{-\lambda}}{x!} \\
\end{align}
$$

ここで観測値$x_1, x_2, …, x_n$が観測される同時確率を$P(x_1,…,x_n|\lambda)$とおくと、$P(x_1,…,x_n|\lambda)$は下記のように表すことができる。
$$
\begin{align}
P(x_1,…,x_n|\lambda) &= \prod_{i=1}^{n} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} \\
&= \lambda^{\sum_{i=1}^{n} x_i} e^{-n \lambda} \prod_{i=1}^{n} (x_i!)^{-1} \\
&= \lambda^{n \bar{x}} e^{-n \lambda} \prod_{i=1}^{n} (x_i!)^{-1}
\end{align}
$$

上記を$\lambda$の関数と見ると、ガンマ分布の形状と一致するので、$\lambda$の事前分布$P(\lambda)$に$Ga(a,1)$を考える。このとき$P(\lambda)$は下記のように表すことができる。
$$
\begin{align}
P(\lambda) = \frac{1}{\Gamma(a)} \lambda^{a-1} e^{-\lambda}
\end{align}
$$

ここまでの内容を元に以下の問題に答えよ。
i) 事後分布$P(\lambda|x_1,x_2,…,x_n) \propto P(x_1,…,x_n|\lambda)P(\lambda)$を計算せよ。ただし、$\propto$を用いて$\lambda$に関連しない項は省略して良い。
ⅱ)
$$
\begin{align}
f(x|a,b) = \frac{b^a}{\Gamma(a)} \lambda^{a-1} e^{-b \lambda}
\end{align}
$$
$Ga(a,b)$の確率密度関数$f(x|a,b)$が上記のように表されるとき、i)の結果より、$\lambda$の事後分布$P(\lambda|x_1,x_2,…,x_n)$がどのようなパラメータのガンマ分布に従うか答えよ。
ⅲ) ガンマ分布$Ga(a,b)$のモーメント母関数を$m(t)$のようにおくと、$\lambda > t$のとき$m(t)$は下記のように表すことができる。
$$
\begin{align}
m(t) = \left( \frac{b}{b-t} \right)^{a}
\end{align}
$$
上記に対し、$\displaystyle \frac{d}{dt}m(t), \frac{d^2}{dt^2}m(t)$を導出せよ。
iv) ⅲ)の結果を元にガンマ分布$Ga(a,b)$の期待値$E[\lambda]$と分散$V[\lambda]$を計算せよ。
v) ⅱ)とiv)の結果を元に、観測値が観測されるにつれて事後分布は事前分布に比べてどのように変化するかを考察せよ。

・解答
i)
下記のように事後分布$P(\lambda|x_1,x_2,…,x_n)$に関して計算を行うことができる。
$$
\large
\begin{align}
P(\lambda|x_1,x_2,…,x_n) & \propto P(x_1,…,x_n|\lambda)P(\lambda) \\
&= \lambda^{n \bar{x}} e^{-n \lambda} \prod_{i=1}^{n} (x_i!)^{-1} \times \frac{1}{\Gamma(a)} \lambda^{a-1} e^{-\lambda} \\
& \propto \lambda^{a + n \bar{x} – 1} e^{-(n+1) \lambda}
\end{align}
$$

ⅱ)
事後分布$P(\lambda|x_1,x_2,…,x_n)$はガンマ分布$Ga(a + n \bar{x}, n+1)$のガンマ分布に従う。

ⅲ)
$\displaystyle \frac{d}{dt}m(t), \frac{d^2}{dt^2}m(t)$は下記のように計算できる。
$$
\large
\begin{align}
\frac{d}{dt}m(t) &= \frac{d}{dt} \left( \frac{b}{b-t} \right)^{a} \\
&= b^a \frac{d}{dt} (b-t)^{-a} \\
&= b^a \times -a(b-t)^{-a-1} \times (-1) = \frac{ab^{a}}{(b-t)^{a+1}} \\
\frac{d^2}{dt^2}m(t) &= \frac{d}{dt} \left( \frac{ab^{a}}{(b-t)^{a+1}} \right) \\
&= ab^{a} \frac{d}{dt} (b-t)^{-(a+1)} \\
&= ab^{a} \times -(a+1) (b-t)^{-(a+2)} \times (-1) = \frac{a(a+1)b^{a}}{(b-t)^{a+2}}
\end{align}
$$

iv)
期待値$E[\lambda]$と分散$V[\lambda]$は下記のように計算することができる。
$$
\large
\begin{align}
E[\lambda] &= \frac{d}{dt}m(t) \Bigr|_{t=0} \\
&= \frac{ab^{a}}{(b-0)^{a+1}} \\
&= \frac{a}{b} \\
V[\lambda] &= \frac{d^2}{dt^2}m(t) \Bigr|_{t=0} – \left( \frac{d}{dt}m(t) \Bigr|_{t=0} \right)^{2} \\
&= \frac{a(a+1)b^{a}}{(b-0)^{a+2}} – \left( \frac{a}{b} \right)^{2} \\
&= \frac{a(a+1)}{b^2} – \frac{a^2}{b^2} \\
&= \frac{a^2+a-a^2}{b^2} \\
&= \frac{a}{b^2}
\end{align}
$$

v)
iv)の結果を元にⅱ)で確認を行なった事前分布$Ga(a, 1)$と事後分布$Ga(a + n \bar{x}, n+1)$を確認する。事前分布の期待値・分散をそれぞれ$E[\lambda], V[\lambda]$、事後分布の期待値・分散をそれぞれ$E[\lambda|x_1,x_2,…,x_n], V[\lambda|x_1,x_2,…,x_n]$とおく。$E[\lambda], E[\lambda|x_1,x_2,…,x_n], V[\lambda], V[\lambda|x_1,x_2,…,x_n]$はそれぞれ下記のように表すことができる。
$$
\large
\begin{align}
E[\lambda] &= \frac{a}{1} = a \\
E[\lambda|x_1,x_2,…,x_n] &= \frac{a + n \bar{x}}{n+1} \\
V[\lambda] &= \frac{a}{1^2} = a \\
V[\lambda|x_1,x_2,…,x_n] &= \frac{a + n \bar{x}}{(n+1)^2}
\end{align}
$$

上記より、$n$が大きくなるにつれて$E[\lambda], E[\lambda|x_1,x_2,…,x_n]$が標本平均の$\bar{x}$に近づき、$V[\lambda], E[\lambda|x_1,x_2,…,x_n]$が$0$に近づくことがわかる。

これは$n$が大きくなるにつれて最尤推定の結果に近づくことを意味し、ベイズの定理を用いた推定が最尤法の拡張であると考えることができることを表している。

・解説
主に下記の内容に基づいて作成を行いました。
https://www.hello-statisticians.com/explain-terms-cat/conjugate_dist1.html

ⅲ)で取り扱ったガンマ分布のモーメント母関数に関しては下記で詳しい導出を行いましたので、こちらも参考になると思います。
https://www.hello-statisticians.com/explain-books-cat/toukeigaku-aohon/ch1_blue_practice.html#19

v)の結果の解釈は最尤法とベイズを比較する際によく出てくるので、事後分布を用いた推定結果の解釈に関してはなるべく行うようにすると良いと思います。

正規分布の共役事前分布の設定と事後分布

・問題
観測値が正規分布に従うと仮定した際の尤度のパラメータに対して正規分布の共役事前分布を考えることで、ベイズの定理を用いて事後分布を計算することができる。

この計算は重要なトピックである一方で平方完成に関連する計算が複雑になり、計算ミスなどが起こりやすい。この問題では以下、関連の計算について可能な限りシンプルに計算ができるように導出の流れの確認を行う。

以下の問題に答えよ。
i) 正規分布$\mathcal{N}(\mu, \sigma^2)$の確率密度関数を$f(x|\mu,\sigma^2)$のようにおくとき、$f(x|\mu,\sigma^2)$を$x, \mu, \sigma^2$の式で表せ。
ⅱ) 確率変数列${X_1,…,X_n}$に対して$X_1,…,X_n \sim \mathcal{N}(\mu, \sigma^2) \quad i.i.d.,$が成立するとき、確率変数列$\{X_1,…,X_n\}$に対応する観測値の列$\{x_1,…,x_n\}$が観測されたときの$\mu$に関する尤度$L(\mu|x_1,…,x_n)$を表せ。
ⅲ) $\displaystyle \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i$のようにおくとき、Fisher-Neymanの分解定理を用いて$\bar{x}$が$\mu$の十分統計量であることを示せ。
iv) $\mu$の事前分布が正規分布$\mathcal{N}(\lambda,\tau^2)$に従うとき、確率密度関数$P(\mu)$を$\lambda,\tau^2$を用いて表せ。
v) $a>0$のとき$g(x) = ax^2 + bx + c$の平方完成を行え。また、$g(x)$が最小となる$x$の値を答えよ。
vi) $x_1,x_2,…,x_n$が観測された際の$\mu$の事後分布の確率密度関数$P(\mu|x_1,…,x_n)$をベイズの定理に基づいて計算し、事後分布が$\displaystyle \mathcal{N} \left( \frac{n \tau^2 \bar{x} + \sigma^2 \lambda}{n \tau^2 + \sigma^2} , \frac{\sigma^2 \tau^2}{n \tau^2 + \sigma^2} \right)$に一致することを確認せよ。
vⅱ) vi)の結果を事後分布の平均に着目して解釈せよ。

・解答
i)
確率密度関数$f(x|\mu,\sigma^2)$は下記のように表すことができる。
$$
\large
\begin{align}
f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ -\frac{(x-\mu)^2}{2 \sigma^2} \right]
\end{align}
$$

ⅱ)
尤度$L(\mu|x_1,…,x_n)$は同時確率に一致するので、$i.i.d.,$であることから下記のように導出することができる。
$$
\large
\begin{align}
L(\mu|x_1,…,x_n) &= f(x_1,…,x_n|\mu,\sigma^2) = \prod_{i=1}^{n} f(x_i|\mu,\sigma^2) \\
&= \prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ -\frac{(x_i-\mu)^2}{2 \sigma^2} \right] \\
&= \left( \frac{1}{\sqrt{2 \pi} \sigma} \right)^{n} \exp \left[ -\sum_{i=1}^{n} \frac{(x_i-\mu)^2}{2 \sigma^2} \right]
\end{align}
$$

ⅲ)
$(x_i-\mu) = (x_i-\bar{x}) + (\bar{x}-\mu)$より、 $\displaystyle \sum_{i=1}^{n} (x_i-\mu)^2$は下記のように変形できる。
$$
\large
\begin{align}
\sum_{i=1}^{n} (x_i-\mu)^2 &= \sum_{i=1}^{n} ((x_i-\bar{x}) + (\bar{x}-\mu))^2 \\
&= \sum_{i=1}^{n} ((x_i-\bar{x})^2 + 2(\bar{x}-\mu)^2 + (x_i-\bar{x})(\bar{x}-\mu)) \\
&= \sum_{i=1}^{n} (x_i-\bar{x})^2 + \sum_{i=1}^{n} (\bar{x}-\mu)^2 + 2(\bar{x}-\mu) \sum_{i=1}^{n}(x_i-\bar{x}) \\
&= \sum_{i=1}^{n} (x_i-\bar{x})^2 + \sum_{i=1}^{n} (\bar{x}-\mu)^2 \\
&= n(x_i-\bar{x})^2 + \sum_{i=1}^{n} (\bar{x}-\mu)^2
\end{align}
$$

よって同時確率密度関数$f(x_1,…,x_n|\mu,\sigma^2)$は下記のように変形できる。
$$
\large
\begin{align}
f(x_1,…,x_n|\mu,\sigma^2) &= \prod_{i=1}^{n} f(x_i|\mu,\sigma^2) \\
&= \left( \frac{1}{\sqrt{2 \pi} \sigma} \right)^{n} \exp \left[ -\sum_{i=1}^{n} \frac{(x_i-\mu)^2}{2 \sigma^2} \right] \\
&= \exp \left[ -\frac{n(\bar{x}-\mu)^2}{2 \sigma^2} \right] \times \left( \frac{1}{\sqrt{2 \pi} \sigma} \right)^{n} \exp \left[ -\sum_{i=1}^{n} \frac{(x_i-\bar{x})^2}{2 \sigma^2} \right]
\end{align}
$$

上記のように同時確率密度関数を$\mu, \bar{x}$の関数と$x_i,\bar{x}$の関数に分解することができたのでFisher-Neymanの定理より$\bar{x}$が$\mu$の十分統計量であることが示せる。

iv)
確率密度関数$P(\mu)$は下記のように表せる。
$$
\large
\begin{align}
f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2 \pi} \tau} \exp \left[ -\frac{(\mu-\lambda)^2}{2 \tau^2} \right]
\end{align}
$$

v)
下記のように$g(x)$の平方完成を行うことができる。
$$
\large
\begin{align}
g(x) &= ax^2+bx+c \\
&= a \left( x^2 + \frac{b}{a}x \right) + c \\
&= a \left( x + \frac{b}{2a} \right)^2 + c – a \left( \frac{b}{2a} \right)^2 \\
&= a \left( x + \frac{b}{2a} \right)^2 + c – \frac{ab^2}{4a^2} \\
&= a \left( x + \frac{b}{2a} \right)^2 + c – \frac{b^2}{4a}
\end{align}
$$

$g(x)$は$a>0$から下に凸の関数であるので、上記より、$\displaystyle x = -\frac{b}{2a}$のとき$g(x)$は最小値を取ることがわかる。

vi)
事後分布の確率密度関数$P(\mu|x_1,…,x_n)$に関しては下記のように考えることができる。
$$
\large
\begin{align}
P(\mu|x_1,…,x_n) & \propto L(\mu|x_1,…,x_n) P(\mu) \\
&= \exp \left[ -\frac{n(\bar{x}-\mu)^2}{2 \sigma^2} \right] \times \left( \frac{1}{\sqrt{2 \pi} \sigma} \right)^{n} \exp \left[ -\sum_{i=1}^{n} \frac{(x_i-\bar{x})^2}{2 \sigma^2} \right] \times \frac{1}{\sqrt{2 \pi} \tau} \exp \left[ -\frac{(\mu-\lambda)^2}{2 \tau^2} \right] \\
& \propto \exp \left[ -\frac{n(\bar{x}-\mu)^2}{2 \sigma^2} \right] \times \exp \left[ -\frac{(\mu-\lambda)^2}{2 \tau^2} \right] \\
&= \exp \left[ -\frac{\sigma^2(\mu-\lambda)^2 + n \tau^2(\mu-\bar{x})^2}{2 \sigma^2 \tau^2} \right] \\
&= \exp \left[ -\frac{(n \tau^2 + \sigma^2)\mu^2 – 2(n \tau^2 \bar{x} + \sigma^2 \lambda)\mu + …}{2 \sigma^2 \tau^2} \right] \\
& \propto \exp \left[ -\frac{n \tau^2 + \sigma^2}{2 \sigma^2 \tau^2} \left(\mu – \frac{n \tau^2 \bar{x} + \sigma^2 \lambda}{n \tau^2 + \sigma^2} \right)^2 \right]
\end{align}
$$

vⅱ)
$\displaystyle \frac{n \tau^2 \bar{x} + \sigma^2 \lambda}{n \tau^2 + \sigma^2}$を解釈するにあたっては、標本平均$\bar{x}$と事前分布の平均の$\lambda$がそれぞれどのような配分で事後分布の平均に作用するかを確認すればよい。

ここで事後分布の平均は$\bar{x}$と$\lambda$を$n \tau^2:\sigma^2$に内分する点であることが確認でき、$n$が大きい時は$\bar{x}$に限りなく近づくことも同時にわかる。

・解説
下記の内容に基づいて作成を行いました。
https://www.hello-statisticians.com/explain-terms-cat/conjugate_dist1.html#i-6

ⅲ)のようにFisher-Neymanの定理の式の形で同時確率密度関数、尤度を表すことで$\mu$に関連する部分の表記をシンプルに行えます。このような表記を用いることで、vi)での計算で$\displaystyle \sum$を考えなくてよく、これにより式が見やすくなります。

vi)の計算にあたってはv)のように全ての項を考えるのではなく、v)の結果を活用しさらに$\propto$を用いて定数項は考えないことで計算をシンプルに行うことができます。

EAP推定量とMAP推定量

・問題
ベイズの定理を用いてパラメータの事後分布を得た場合、予測分布のように事後分布をそのまま用いる場合もあるが、事後分布を$1$つの値に要約して用いる場合も多い。この要約にあたっては事後分布の期待値を考えるEAP推定量と事後分布が最大となる点に着目するMAP推定量がよく用いられる。EAP推定量は期待事後推定量(Expected a Posterior Estimator)の略で、MAP推定量は最大事後確率推定量(Maximum a Posterior Estimator)の略であることも合わせて抑えておくと良い。

この問題では以下、EAP推定量とMAP推定量に関して具体的な事後分布に基づいて確認するにあたって、「ベータ分布」、「ガンマ分布」、「正規分布」に対してそれぞれEAP推定量とMAP推定量の導出を行う。「ベータ分布$\mathrm{Be}(\alpha,\beta)$」、「ガンマ分布$\mathrm{Ga}(\alpha,\beta)$」、「正規分布$\mathcal{N}(\mu,\sigma^2)$」の確率密度関数を$f_1(x),f_2(x),f_3(x)$とおくとそれぞれ下記のように表すことができる。
$$
\large
\begin{align}
f_1(x) &= \frac{1}{B(\alpha,\beta)} x^{\alpha-1}(1-x)^{\beta-1} = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1} \\
f_2(x) &= \frac{1}{\beta^{\alpha} \Gamma(\alpha)} x^{\alpha-1} \exp{\left( -\frac{x}{\beta} \right)} \\
f_3(x) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{\left( -\frac{(x-\mu)^2}{2 \sigma^2} \right)}
\end{align}
$$

ここまでの内容に基づいて下記の問いにそれぞれ答えよ。
i) EAP推定量とMAP推定量は記述統計における「代表値」と対応して理解することができるが、EAP推定量とMAP推定量はそれぞれどのような代表値に対応するか答えよ。
ⅱ) 「二項分布の共役事前分布と事後分布の解釈」では、二項分布$\mathrm{Bin}(n,p)$の確率パラメータ$p$の事後分布にベータ分布$\mathrm{Be}(\alpha+x,\beta+n-x)$が得られた。このとき$p$のEAP推定量$\hat{p}_{EAP}$とMAP推定量$\hat{p}_{EAP}$を$n,x,\alpha,\beta$を用いて表せ。
ⅲ) ⅱ)と同様に事後分布にベータ分布$\mathrm{Be}(\alpha+x,\beta+n-x)$が得られたとき、$p$のMAP推定量$\hat{p}_{MAP}$を$n,x,\alpha,\beta$を用いて表せ。
iv) ⅱ)、ⅲ)で得られた$\hat{p}_{EAP}$と$\hat{p}_{MAP}$を最尤推定量$\displaystyle \hat{p}_{MLE}=\frac{x}{n}$に基づいて解釈せよ。
v) $\lambda \sim \mathrm{Ga}(\alpha,\beta)$のとき、$\hat{\lambda}_{EAP}$と$\hat{\lambda}_{MAP}$をそれぞれ計算せよ。
vi) $\mu \sim \mathcal{N}(\mu_N,\sigma_N^2)$のとき、$\hat{\mu}_{EAP}$と$\hat{\mu}_{MAP}$をそれぞれ答えよ。

・解答
i)
EAP推定量は平均に対応し、MAP推定量はモードに対応すると考えることができる。

ⅱ)
$p$のEAP推定量$\hat{p}_{EAP}$は下記のように得られる。
$$
\large
\begin{align}
\hat{p}_{EAP} &= \int_{0}^{1} p \times f_1(p) dp \\
&= \frac{1}{B(a+x,\beta+n-x)} \int_{0}^{1} p \times p^{\alpha+x-1} (1-p)^{\alpha+n-x-1} dp \\
&= \frac{1}{B(a+x,\beta+n-x)} \int_{0}^{1} p^{\alpha+x+1-1} (1-p)^{\alpha+n-x-1} dp \\
&= \frac{B(a+x+1,\beta+n-x)}{B(a+x,\beta+n-x)} \\
&= \frac{\Gamma(\alpha+x+1)\Gamma(\beta+n-x)}{\Gamma(n+\alpha+\beta+1)} \times \frac{\Gamma(n+\alpha+\beta)}{\Gamma(\alpha+x)\Gamma(\beta+n-x)} \\
&= \frac{(\alpha+x)!(\beta+n-x-1)!}{(n+\alpha+\beta)!} \times \frac{(n+\alpha+\beta-1)!}{(\alpha+x-1)!(\beta+n-x-1)!} \\
&= \frac{\alpha+x}{n+\alpha+\beta}
\end{align}
$$

ⅲ)
$p$のMAP推定量$\hat{p}_{MAP}$は下記のように得られる。
$$
\large
\begin{align}
\frac{\partial \log{(P(p|x))}}{\partial p} &= \frac{\partial}{\partial p} \log{\left[ \frac{1}{B(\alpha+x-1,\beta+n-x)} p^{\alpha+x-1}(1-p)^{\beta+n-x-1} \right]} \\
&= \frac{\partial}{\partial p} \left[ (\alpha+x-1)\log{p} + (\beta+n-x-1)\log{(1-p)} – \log{B(\alpha+x-1,\beta+n-x)} \right] \\
&= \frac{\alpha+x-1}{p} – \frac{\beta+n-x-1}{1-p} \\
&= \frac{(\alpha+x-1)(1-p) – (\beta+n-x-1)p}{p(1-p)} \\
&= \frac{\alpha+x-1 – (n+\alpha+\beta-2)p}{p(1-p)} \\
\hat{p}_{MAP} &= \frac{\alpha+x-1}{n+\alpha+\beta-2}
\end{align}
$$

iv)
最尤推定量$\displaystyle \hat{p}_{MLE}=\frac{x}{n}$は$n$回の試行の中で事象が$x$回観測されたのでパラメータを$\displaystyle \hat{p}_{MLE}=\frac{x}{n}$のように推定したと考えられる。
事後分布に基づく推定は予めサンプルがある上で観測された割合を計算したと見ることもでき、EAP推定量は「事前に$\alpha+\beta$回中$\alpha$回計測された場合」、MAP推定量は「事前に$\alpha+\beta-2$回中$\alpha-1$回計測された場合」とそれぞれ解釈できる。

v)
$\hat{\lambda}_{EAP}$は下記のように計算できる。
$$
\large
\begin{align}
\hat{\lambda}_{EAP} &= \int_{0}^{\infty} \lambda \times f_2(\lambda) d \lambda \\
&= \frac{1}{\beta^{\alpha} \Gamma(\alpha)} \int_{0}^{\infty} \lambda \times \lambda^{\alpha-1} \exp{\left( -\frac{\lambda}{\beta} \right)}d \lambda \\
&= \frac{1}{\beta^{\alpha} \Gamma(\alpha)} \int_{0}^{\infty} \lambda^{(\alpha+1)-1} \exp{\left( -\frac{\lambda}{\beta} \right)}d \lambda \\
&= \frac{1}{\beta^{\alpha} \Gamma(\alpha)} \times \beta{\alpha+1} \Gamma(\alpha+1) \\
&= \alpha \beta
\end{align}
$$

$\hat{\lambda}_{MAP}$の計算にあたっては$\log{f_2(\lambda)}$を考え、$\lambda$で微分を行えばよい。$\log{f_2(\lambda)}$は下記のように表すことができる。
$$
\large
\begin{align}
\log{f_2(\lambda)} &= \log{ \left[ \frac{1}{\beta^{\alpha} \Gamma(\alpha)} \lambda^{\alpha-1} \exp{\left( -\frac{\lambda}{\beta} \right)} \right] } \\
&= (\alpha-1) \log{\lambda} – \frac{\lambda}{\beta} – \log{\left[ \beta^{\alpha} \Gamma(\alpha) \right]}
\end{align}
$$

上記を$\lambda$で偏微分を行うことで下記が得られる。
$$
\large
\begin{align}
\frac{\partial \log{f_2(\lambda)}}{\partial \lambda} &= \frac{\alpha-1}{\lambda} – \frac{1}{\beta} \\
&= \frac{\alpha \beta – \beta – \lambda}{\lambda \beta} \\
\hat{\lambda}_{MAP} &= (\alpha-1) \beta
\end{align}
$$

vi)
正規分布の確率密度関数の形状より、$\hat{\mu}_{EAP}=\mu_N, \hat{\mu}_{MAP}=\mu_N$である。

・解説
ⅱ)〜iv)では二項分布の確率パラメータ$p$の事後分布に基づいてEAP推定量、MAP推定量の計算を行いました。v)、vi)の事後分布のパラメータはそれぞれ「ポアソン分布の共役事前分布と事後分布の解釈」と「正規分布の共役事前分布の設定と事後分布」で取り扱いました。ⅱ)〜iv)では事後分布のパラメータを元に推定量の計算を行いましたが、このように計算すると式が複雑になるので、v)、vi)のように「①それぞれの確率分布に関する推定量の計算」と「②事後分布のパラメータの計算」を分けて考えるとシンプルに取り扱えて良いのではないかと思います。v)、vi)では簡略化にあたって②は省略し、①のみ取り扱いました。
また、EAP推定量は積分、MAP推定量は微分を用いることでそれぞれ計算できることは抑えておくと良いと思います。
「統計検定準$1$級のワークブック」ではEAP推定量を「ベイズ推定量」のように表現されますが、ベイズ推定が何を指すのかがわかりにくくなるので、この問題では「数理統計学 統計的推論の基礎」に基づいて「EAP推定量」と表記しました。

発展問題

共役事前分布と指数型分布族

予測分布

・問題
パラメータ$\theta$に関する事後分布は観測された標本$x_1,\cdots,x_n$を元に$P(\theta|x_1,\cdots,x_n)$を考える。ここで得られた事後分布を元に新規の標本に関して予測を行うことを考える。

予測を行うにあたっては、「EAP推定量とMAP推定量」で考えたように事後分布を元にパラメータの推定量を計算し、推定量を元に予測を行うのがシンプルである。が、新規のサンプル$x_{n+1}$に関して下記を考えることで予測を行うことも可能である。
$$
\large
\begin{align}
P(x_{n+1}|x_1,\cdots,x_n) = \int P(x_{n+1}|\theta) P(\theta|x_1,\cdots,x_n) d \theta
\end{align}
$$

上記に基づいて定められる$P(x_{n+1}|x_1,\cdots,x_n)$を予測分布(Predictive distribution)という。以下、標本がポアソン分布に基づいて観測された際の予測分布の導出に関して演習形式で確認を行う。下記の問いにそれぞれ答えよ。
i) ポアソン分布$\mathrm{Po}(\lambda)$の確率関数を$p(x)$とおくとき、$p(x)$を$x, \lambda$の式で表せ。
ⅱ) $x_1=2, x_2=3, x_3=5, x_4=2, x_5=3$が観測されたとき、i)の式を元に尤度$L(\lambda)$を表せ。
ⅲ) パラメータ$\lambda$の事前分布にガンマ分布$\displaystyle \mathrm{Ga} \left( 3,\frac{1}{2} \right)$を考える。この確率密度関数を$P(\lambda)$とおくとき$P(\lambda)$を表せ。
iv) ⅱ)の尤度$L(\lambda)=P(x_1,\cdots,x_5|\lambda)$とⅲ)の事前確率$P(\lambda)$を用いて$\lambda$の事後分布を求めよ。
v) iv)で求めた事後分布の確率密度関数$P(\lambda|x_1,\cdots,x_5)$を表せ。
vi) v)の結果などを用いて$(1)$式を元に$P(x_{n+1}|x_1,\cdots,x_5)$を計算せよ。
vⅱ) 負の二項分布$\mathrm{NB}(r,p)$の確率関数を$p(x)$とおくと、$p(x)$は下記のように表せる。
$$
\begin{align}
p(x) &= {}_r H_{x} p^{r} (1-p)^{x} \\
&= {}_{r+x-1} C_{x} p^{r} (1-p)^{x} \\
&= \frac{(r+x-1)!}{x! (r-1)!} p^{r} (1-p)^{x} \quad (2)
\end{align}
$$
$(2)$式を元にvi)の導出結果の$x_6$に対応する確率変数が従う負の二項分布のパラメータを答えよ。

・解答
i)
確率関数$p(x)$は下記のように表すことができる。
$$
\large
\begin{align}
p(x) = \frac{\lambda^{x} e^{-\lambda}}{x!}
\end{align}
$$

ⅱ)
尤度$L(\lambda)$は標本が観測される同時確率に一致するので下記のように考えることができる。
$$
\large
\begin{align}
L(\lambda) &= \prod_{i=1}^{5} p(x_i) = \prod_{i=1}^{5} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} \\
&= \lambda^{\sum_{i=1}^{5} x_i} e^{-5 \lambda} \left( \prod_{i=1}^{5} x_i! \right)^{-1} \\
&= \lambda^{15} e^{-5 \lambda} \left( \prod_{i=1}^{5} x_i! \right)^{-1}
\end{align}
$$

ⅲ)
パラメータ$\lambda$の事前分布$\displaystyle \mathrm{Ga} \left( 3,\frac{1}{2} \right)$の確率密度関数$P(\lambda)$は下記のように表すことができる。
$$
\large
\begin{align}
P(\lambda) = \frac{2^{3}}{\Gamma(3)} \lambda^{3-1} e^{-2 \lambda}
\end{align}
$$

iv)
事後分布の確率密度関数を$P(\lambda|x_1,\cdots,x_5)$とおくと、$P(\lambda|x_1,\cdots,x_5)$は下記のように考えることができる。
$$
\large
\begin{align}
P(\lambda|x_1,\cdots,x_n) & \propto P(x_1,\cdots,x_5|\lambda)P(\lambda) \\
& \propto \lambda^{15} e^{-5 \lambda} \times \lambda^{3-1} e^{-2 \lambda} \\
&= \lambda^{18-1} e^{-7 \lambda}
\end{align}
$$

上記より、$\lambda$の事後分布はガンマ分布$\displaystyle \mathrm{Ga} \left( 18,\frac{1}{7} \right)$であると考えることができる。

v)
$P(\lambda|x_1,\cdots,x_n)$はガンマ分布$\displaystyle \mathrm{Ga} \left( 18,\frac{1}{7} \right)$の確率密度関数であるので下記のように表せる。
$$
\large
\begin{align}
P(\lambda|x_1,\cdots,x_n) = \frac{7^{18}}{\Gamma(18)} \lambda^{18-1} e^{-7 \lambda}
\end{align}
$$

vi)
$(1)$式を元に$P(x_{6}|x_1,\cdots,x_5)$は下記のように計算できる。
$$
\large
\begin{align}
P(x_{6}|x_1,\cdots,x_5) &= \int_{0}^{\infty} P(x_{6}|\lambda) P(\lambda|x_1,\cdots,x_5) d \lambda \quad (1) \\
&= \int_{0}^{\infty} \frac{\lambda^{x_6} e^{-\lambda}}{x_6!} \times \frac{7^{18}}{\Gamma(18)} \lambda^{18-1} e^{-7 \lambda} d \lambda \\
&= \frac{7^{18}}{x_6! \Gamma(18)} \int_{0}^{\infty} \lambda^{18+x_6-1} e^{-8 \lambda} d \lambda \\
&= \frac{7^{18}}{x_6! \Gamma(18)} \times \frac{\Gamma(18+x_6)}{8^{18+x_6}} \\
&= \frac{(18+x_6-1)!}{x_6! (18-1)!} \cdot \left( \frac{7}{8} \right)^{18} \cdot \left( \frac{1}{8} \right)^{x_6}
\end{align}
$$

vⅱ)
$(2)$式に対し、$\displaystyle r=18, x=x_6, p=\frac{7}{8}$で置き換えることでvi)の$P(x_{6}|x_1,\cdots,x_5)$の式が得られる。よって、$x_6$に対応する確率変数は負の二項分布$\displaystyle \mathrm{NB} \left( 18,\frac{7}{8} \right)$に従うと考えられる。

・解説
問題の作成にあたっては下記の問題の設定を主に用いました。

参考

統計検定2級問題解説 ~2021年6月実施~ (問13~問22)

過去問題

過去問題は統計検定公式問題集が問題と解答例を公開しています。こちらを参照してください。


問13 解答

(推定量)

$\boxed{ \ \mathsf{17}\ }$ ①

Ⅰ.母集団の特徴づける定数のことを母数というが、母集団から測定された標本データをもとに、確率分布の(現実には測定できない)母数を推定した数量を推定量という。推定量は標本データの関数として表されるが、確率変数の関数は確率変数なので、推定量は確率変数となる。
Ⅱ.推定量$\hat\theta_n$(推定量の標本分布が標本数$n$によって変化するものとする)が母数$\theta$に確率収束している場合、$\hat\theta$を一致推定量という。
Ⅲ.推定量$\hat\theta$の期待値が常に母数$\theta$に等しくなる場合、$\hat\theta$を不偏推定量という。例えば母分散の推定量である標本分散は一致推定量であるが不偏推定量ではない。


問14 解答

(推定量の期待値、分散)

[1]

$\boxed{ \ \mathsf{18}\ }$ ⑤

$E[\hat\theta]=\theta$となる場合、$\hat\theta$は不偏推定量という。
$$\begin{align}
E[\hat\mu_1]&=E\left[\frac1n\sum_{i=1}^nX_i\right]=\frac1n\sum_{i=1}^nE[X_i]=\frac1n\times n\mu=\mu\\
E[\hat\mu_2]&=E\left[\frac12(X_1+X_2)\right]=\frac12(E[X_1]+E[X_2])=\frac12(\mu+\mu)=\mu\\
E[\hat\mu_3]&=E[X_1]=\mu\\
E[\hat\mu_4]&=E\left[\frac2{n(n+1)}\sum_{i=1}^niX_i\right]=\frac2{n(n+1)}\sum_{i=1}^niE[X_i]=\frac2{n(n+1)}\sum_{i=1}^ni\mu\\&=\frac2{n(n+1)}\times\frac{n(n+1)}{2}\times\mu=\mu
\end{align}$$

[2]

$\boxed{ \ \mathsf{19}\ }$ ①

$$\begin{align}
V[\hat\mu_1]&=V\left[\frac1n\sum_{i=1}^nX_i\right]=\frac1{n^2}\sum_{i=1}^nV[X_i]=\frac1{n^2}\times n\sigma^2=\frac1n\sigma^2\\
V[\hat\mu_2]&=V\left[\frac12(X_1+X_2)\right]=\frac14(V[X_1]+V[X_2])=\frac14(\sigma^2+\sigma^2)=\frac12\sigma^2\\
V[\hat\mu_3]&=V[X_1]=\sigma^2\\
V[\hat\mu_4]&=V\left[\frac2{n(n+1)}\sum_{i=1}^niX_i\right]=\frac4{n^2(n+1)^2}\sum_{i=1}^ni^2V[X_i]=\frac4{n^2(n+1)^2}\sum_{i=1}^ni^2\sigma^2\\&=\frac4{n^2(n+1)^2}\times\frac{n(n+1)(2n+1)}{6}\times\sigma^2=\frac{2(2n+1)}{3n(n+1)}\sigma^2
\end{align}$$
ここで、$n\gt3$なので、$\displaystyle \frac1n\lt\frac12\lt1$
また、$\displaystyle \frac1n=\frac{3(n+1)}{3n(n+1)}<\frac{2(2n+1)}{3n(n+1)}$となることから$V[\hat\mu_1]$が最小となる。


問15 解答

(正規母集団の区間推定とサンプルサイズ)

[1]

$\boxed{ \ \mathsf{20}\ }$ ④

正規母集団から抽出した標本の標本平均は$\bar{X}\sim N(\mu,\sigma^2/n)$なので,$\begin{align}\frac{\bar{X}-\mu}{\sqrt{\sigma^2/n}}\end{align}$は標準正規分布に従う。よって「標準正規分布の上側確率」の表から
$$P\left(|\bar{X}-\mu|\le1.96{\sqrt{\sigma^2/n}}\right)=0.95$$
したがって、真の$\mu$が含まれる確率が$95\%$となる区間($\mu$の$95\%$信頼区間)は以下の通りとなる。
$$\begin{align}
\bar{X}-1.96\frac{\sigma}{\sqrt{n}}\le&\mu\le\bar{X}+1.96\frac{\sigma}{\sqrt{n}}\\
5.25-1.96\times\frac{12}{\sqrt{100}}\le&\mu\le5.25+1.96\times\frac{12}{\sqrt{100}}\\
2.90\le&\mu\le7.60
\end{align}$$

[2]

$\boxed{ \ \mathsf{21}\ }$ ③

[1]から$\mu$の$95\%$信頼区間の幅は$\displaystyle\ 2\times1.96\frac{\sigma}{\sqrt{n}}$。これを$4$以下にしたいので、
$$2\times1.96\times\frac{12}{\sqrt{n}}\le4\ \Rightarrow\ n\ge\left(\frac{2\times1.96\times12}{4}\right)^2=138.3$$


問16 解答

(単回帰モデル、最小二乗法)

[1]

$\boxed{ \ \mathsf{22}\ }$ ①

最小二乗法は実際の値$y_i$と回帰式によって予測された値$\hat{y}_i$との差(残差)の二乗和を最小にするように回帰係数を求める手法である。残差の二乗和(残差平方和)は
$$S=\sum_{i=1}^n(y_i-\hat{y}_i)^2=\sum_{i=1}^n(y_i-\hat\beta x_i)^2=\sum_{i=1}^n(y_i^2-2\hat\beta x_iy_i+\hat\beta^2x_i^2)$$
$S$を最小とする$\hat\beta$を求めるために、$\hat\beta$で偏微分し$0$に等しいとすると、
$$\begin{eqnarray}
\frac{\partial S}{\partial\hat\beta}=\sum_{i=1}^n(-2x_iy_i+2\hat\beta x_i^2)=0\\
\hat\beta\sum_{i=1}^nx_i^2=\sum_{i=1}^nx_iy_i\\
\therefore\ \hat\beta=\frac{\sum_{i=1}^nx_iy_i}{\sum_{i=1}^nx_i^2}
\end{eqnarray}$$

[2]

$\boxed{ \ \mathsf{23}\ }$ ②

Ⅰ.[1]の結果から、一般的に
$$\hat\beta=\frac{\sum_{i=1}^nx_iy_i}{\sum_{i=1}^nx_i^2}\neq\frac{\sum_{i=1}^ny_i}{\sum_{i=1}^nx_i}$$
であるから、
$$\begin{eqnarray}
\hat\beta\sum_{i=1}^nx_i&\neq&\sum_{i=1}^ny_i\\
\sum_{i=1}^ny_i-\hat\beta\sum_{i=1}^nx_i&\neq&0\\
\sum_{i=1}^n\hat u_i&\neq&0
\end{eqnarray}$$
Ⅱ.[1]から
$$\begin{eqnarray}
\hat\beta\sum_{i=1}^nx_i^2=\sum_{i=1}^nx_iy_i\\
\sum_{i=1}^nx_i(y_i-\hat\beta x_i)=0\\
\sum_{i=1}^nx_i(y_i-\hat y_i)=0\\
\sum_{i=1}^nx_i\hat u_i=0
\end{eqnarray}$$
Ⅲ.[1]の結果から、一般的に
$$\hat\beta=\frac{\sum_{i=1}^nx_iy_i}{\sum_{i=1}^nx_i^2}\neq\frac{\sum_{i=1}^ny_i}{\sum_{i=1}^nx_i}$$
であるから、
$$\begin{eqnarray}
\frac1n\hat\beta\sum_{i=1}^nx_i&\neq&\frac1n\sum_{i=1}^ny_i\\
\frac1n\sum_{i=1}^n\hat\beta x_i&\neq&\bar{y}\\
\frac1n\sum_{i=1}^n\hat y_i&\neq&\bar{y}
\end{eqnarray}$$
Ⅳ.Ⅲ.と同じく
$$\begin{eqnarray}
\frac1n\hat\beta\sum_{i=1}^nx_i&\neq&\frac1n\sum_{i=1}^ny_i\\
\hat\beta\bar{x}&\neq&\bar{y}
\end{eqnarray}$$

※定数項を含む単回帰モデル$$y_i=\alpha+\beta x_i+u_i$$の場合、残差平方和は
$$S=\sum_{i=1}^n(y_i-\hat{y}_i)^2=\sum_{i=1}^n(y_i-\hat\alpha-\hat\beta x_i)^2$$
$S$を最小とする$\hat\beta$を求めるために、$\hat\alpha,\hat\beta$で偏微分し$0$に等しいとすると、
$$
\sum_{i=1}^n(y_i^2-2\hat\alpha y_i-2\hat\beta x_iy_i+\hat\alpha^2+2\hat\alpha\hat\beta x_i+\hat\beta^2 x_i^2)
$$
$$\begin{eqnarray}
\frac{\partial S}{\partial\hat\alpha}=\sum_{i=1}^n(-2y_i+2\hat\alpha+2\hat\beta x_i)&=&0\\
n\hat\alpha+\hat\beta\sum_{i=1}^nx_i&=&\sum_{i=1}^ny_i&\cdots(A)\\
\frac{\partial S}{\partial\hat\beta}=\sum_{i=1}^n(-2x_iy_i+2\hat\alpha x_i+2\hat\beta x_i^2)&=&0\\
\hat\alpha\sum_{i=1}^nx_i+\hat\beta\sum_{i=1}^nx_i^2&=&\sum_{i=1}^nx_iy_i&\cdots(B)
\end{eqnarray}$$
$(A)$から
$$\hat\alpha=\frac1n\sum_{i=1}^ny_i-\frac1n\hat\beta\sum_{i=1}^nx_i=\bar{y}-\hat\beta\bar{x}$$
$(B)$から
$$\begin{eqnarray}
(\bar{y}-\hat\beta\bar{x})n\bar{x}+\hat\beta\sum_{i=1}^nx_i^2=\sum_{i=1}^nx_iy_i\\
\hat\beta(\sum_{i=1}^nx_i^2-n\bar{x}^2)=\sum_{i=1}^nx_iy_i-n\bar{x}\bar{y}\\
\hat\beta=\frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})^2}\\
\end{eqnarray}$$
さらに、(A)から
$$\begin{eqnarray}
\sum_{i=1}^ny_i-n\hat\alpha-\hat\beta\sum_{i=1}^nx_i&=&0\\
\sum_{i=1}^n(y_i-\hat\alpha-\hat\beta x_i)&=&0\\
\sum_{i=1}^n(y_i-\hat y_i)&=&0\\
\sum_{i=1}^n\hat u_i&=&0
\end{eqnarray}$$
また
$$\begin{eqnarray}
\bar y=\frac1n\sum_{i=1}^ny_i=\frac1n\sum_{i=1}^n(\hat\alpha+\hat\beta x_i)=\frac1n\sum_{i=1}^n\hat y\\
\bar y=\frac1n\sum_{i=1}^ny_i=\hat\alpha+\frac1n\hat\beta\sum_{i=1}^n x_i=\hat\alpha+\hat\beta\bar x
\end{eqnarray}$$
$(B)$から
$$\begin{eqnarray}
\sum_{i=1}^nx_iy_i-\hat\alpha\sum_{i=1}^nx_i-\hat\beta\sum_{i=1}^nx_i^2&=&0\\
\sum_{i=1}^nx_i(y_i-\hat\alpha-\hat\beta x_i)&=&0\\
\sum_{i=1}^nx_i(y_i-\hat y_i)&=&0\\
\sum_{i=1}^nx_i\hat u_i&=&0
\end{eqnarray}$$
[2]のⅠ.Ⅲ.Ⅳ.に相当する関係はいずれも$(A)$の式から導き出されるもので、定数項を含まないモデルでは$(A)$に相当する条件がなく、Ⅰ.Ⅲ.Ⅳ.の関係は成り立たない。


問17 解答

(母比率の区間推定)

[1]

$\boxed{ \ \mathsf{24}\ }$ ②

成功確率$p$の試行を$n$回行うときに成功する回数$X$は二項分布$B(n,p)$に従う。
  $\therefore\ \ E(X)=np,\ V(X)=np(1-p)$
このとき,$n$がある程度大きいときは,中心極限定理によって,$B(n,p)$は正規分布$N(np,np(1-p))$に近似できる。よって,$X$を標準化すると標準正規分布$N(0,1)$に従う。$$Z=\frac{X-np}{\sqrt{np(1-p)}}=\frac{X/n-p}{\sqrt{\frac{p(1-p)}n}}\sim N(0,1)$$ここで,標本平均 $\hat p=x/n$は$p$の一致推定量なので,$n$が十分大きいとき$p$は$\hat p$に置き換えられる。
したがって,母比率の$100(1-\alpha)\%$信頼区間は,標準正規分布の上側 $100\alpha/2\%$ 点を $z_{\alpha/2}$とすると,$$P\left(\hat p-z_{\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}n}\le p\le\hat p+z_{\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}n}\right)=1-\alpha$$
$500$ 回画びょうを投げて $284$ 回表が出たので,比率の推定値は $\hat p=284/500=0.568$
これから,表が出る確率の$95\%$信頼区間は,$n=500$,$\hat p=0.568$,$\alpha=0.05$として$$\begin{align}\hat p\pm z_{\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}n}=&0.568\pm1.96\times\sqrt{\frac{0.568\times(1-0.568)}{500}}\\=&0.568\pm0.043\\=&[0.525,0.611]\end{align}$$

[2]

$\boxed{ \ \mathsf{25}\ }$ ①

・確率 $p=1/2$ の試行を $n=8$ 回行って成功する回数の分布なので,二項分布 $B(n,p)=B(8,1/2)$ に従う。$$\begin{align}P(X=4)=&{}_8\mathrm{C}_4\times(1/2)^4\times(1-1/2)^{8-4}\\=&\frac{8\times7\times6\times5}{4\times3\times2\times1}\times(1/2)^4\times(1/2)^4=0.273\end{align}$$
・帰無仮説$H_0:p=1/2$、対立仮説$H_1:p\gt1/2$として、$X\ge c_1$のとき$H_0$を棄却する検定は、片側検定となるので、$X=7$のときの$P_-$値は
$$P(X\ge7|H_0)=P(X=7)+P(X=8)=0.031+0.004=0.035$$
・帰無仮説$H_0:p=1/2$、対立仮説$H_1:p\neq1/2$として、$|X-4|\ge c_2$のとき$H_0$を棄却する検定は、両側検定となり確率分布が対称であるので、実現値が$X=7$のときの$P_-$値は
$$\begin{align}P(|X-4|\ge3|H_0)=&P(X=0)+P(X=1)+P(X=7)+P(X=8)\\=&0.004+0.0031+0.031+0.004=0.070\end{align}$$


問18 解答

(母平均の差の検定(分散未知であるが等分散))

$\boxed{ \ \mathsf{26}\ }$ ②

(ア)正規母集団$N(\mu_1,\sigma^2)$から抽出した標本平均$\bar{X}=\frac1m\sum_{i=1}^mX_i$は$N(\mu_1,\sigma^2/m)$に従い、正規母集団$N(\mu_2,\sigma^2)$から抽出した標本平均$\bar{Y}=\frac1n\sum_{i=1}^nY_i$は$N(\mu_2,\sigma^2/n)$に従う。
したがって、正規分布の再生性から、標本平均の差$\bar{X}-\bar{Y}$は$N(\mu_1-\mu_2,\sigma^2/m+\sigma^2/n)$に従う。その結果、$\bar{X}-\bar{Y}$を標準化して$$A=\frac{\bar{X}-\bar{Y}-(\mu_1-\mu_2)}{\sqrt{\frac{\sigma^2}{m}+\frac{\sigma^2}{n}}}\sim N(0,1)$$となる。
(イ)群$1,2$について$$\begin{eqnarray}
\sum_{i=1}^m\frac{(X_i-\bar{X})^2}{\sigma^2}=\frac{(m-1)U_X^2}{\sigma^2}\sim \chi^2(m-1)\\
\sum_{i=1}^n\frac{(Y_i-\bar{Y})^2}{\sigma^2}=\frac{(n-1)U_Y^2}{\sigma^2}\sim \chi^2(n-1)
\end{eqnarray}$$したがって、$\chi^2$分布の再生性から、$$B=\frac{(m-1)U_X^2+(n-1)U_Y^2}{\sigma^2}\sim \chi^2(m+n-2)$$となる。
(ウ)独立な$2$つの確率変数$Z\sim N(0,1)$と$W\sim\chi^2(m)$があるとき、$$\frac{Z}{\sqrt{W/m}}$$は自由度$m$の$t$分布に従う。よって、(ア)と(イ)から$$T=\frac{A}{\sqrt{\frac{B}{m+n-2}}}$$は自由度$m+n-2$の$t$分布に従う。


問19 解答

(独立性の検定)

独立性の検定は,2つの属性$A,B$が独立かどうかの検定。
属性$A$のカテゴリが$A_i$,属性$B$のカテゴリが$B_j$の観測度数を$f_{ij}=O_{ij}$とし,
 $f_{i\cdot}=\sum_jf{ij}$を$i$行の度数合計(行和),
 $f_{\cdot j}=\sum_if{ij}$を$j$列の度数合計(列和),
 $f_{\cdot\cdot}=\sum_i\sum_jf_{ij}=\sum_if_{i\cdot}=\sum_jf_{\cdot j}=n$を全度数合計という。
属性$A,B$が独立という帰無仮説は,$H_0:P(A\cap B)=P(A)P(B)$が成り立つことである。
ここで,カテゴリ$A_i,B_j$の出現確率はそれぞれ$f_{i\cdot}/n,f_{\cdot j}/n$であるので,$H_0$のもとで,属性$A$のカテゴリが$A_i$,属性$B$のカテゴリが$B_j$の期待度数は$$E_{ij}=n(f_{i\cdot}/n)(f_{\cdot j}/n)=f_{i\cdot}f_{\cdot j}/n$$となる。
帰無仮説$H_0$の下で,次検定の統計量$\chi^2$は度数が大きいときに近似的に$\chi^2$分布に従う。行和と列和が固定されていることから自由度は$(r$(行の数)$-1)\times(c$(列の数)$-1)$となる。$$\chi^2=\sum_{i=1}^r\sum_{j=1}^c\frac{(O_{ij}-E_{ij})^2}{E_{ij}}\sim\chi^2((r-1)(c-1))$$有意水準$100\alpha\%$で帰無仮説が棄却されるには,$\chi^2$分布の上側$\alpha$点より上で求めた$\chi^2$統計量が大きくなればよい。

$\boxed{ \ \mathsf{27}\ }$ ③

期待度数は
 喫煙歴あり・心筋梗塞あり $10\times15/20=7.5$
 喫煙歴あり・心筋梗塞なし $10\times15/20=7.5$
 喫煙歴なし・心筋梗塞あり $10\times5/20=2.5$
 喫煙歴なし・心筋梗塞なし $10\times5/20=2.5$
よって、$\chi^2$統計量の実現値は
$$\chi^2=\frac{(9-7.5)^2}{7.5}+\frac{(6-7.5)^2}{7.5}+\frac{(1-2.5)^2}{2.5}+\frac{(4-2.5)^2}{2.5}=2.40$$
$\chi^2$統計量は帰無仮説の下で近似的に自由度$(2-1)\times(2-1)=1$の$\chi^2$分布に従う。
ここで、確率変数$W$が自由度$1$の$\chi^2$分布に従うとき、標準正規分布に従う確率変数$Z$を用いて$W=Z^2$と表わされる。したがって$P_-$値は「正規標準分布の上側確率」の表を用いて、
$$P(W\gt2.40)=P(|Z|\gt\sqrt{2.40})=P(|Z|\gt1.55)=2\times0.0606=0.1212$$


問20 解答

(第一種の過誤)

真実
帰無仮説が正しい対立仮説が正しい
検定の結果帰無仮説を棄却しない
(対立仮説が正しいとは言えない)
正しい第二種の過誤(β)
帰無仮説を棄却する
(対立仮説が正しい)
第一種の過誤(α)
有意水準
正しい
検出力(1-β)

[1]

$\boxed{ \ \mathsf{28}\ }$ ③

$X_j\sim N(\mu_j,1),\ X_k\sim N(\mu_k,1)$ であるから、$X_j-X_k\sim N(\mu_j-\mu_k,2)$
よって、帰無仮説$H_0:\mu_j=\mu_k$の下で$$Z=\frac{X_j-X_k-(\mu_j-\mu_k)}{\sqrt{2}}=\frac{X_j-X_k}{\sqrt{2}}\sim N(0,1)$$が成り立つ。したがって、第1種過誤の確率$\alpha_{12}(1.96\sqrt{2})$の値は$$\alpha_{12}(1.96\sqrt{2})=P(|X_j-X_k|\gt1.96\sqrt{2})=P(|Z|\gt1.96)=2\times0.025=0.050$$

[2]

$\boxed{ \ \mathsf{29}\ }$ ④

$\alpha_{12}(z)$が$(5/3)\%$となるような$z$を定める。
$$\begin{eqnarray}\alpha_{12}(z)=P(|X_j-X_k|\gt z)=P(|Z|\gt z/\sqrt{2})&=&0.05/3\\P(Z\gt z/\sqrt{2})&=&0.05/6=0.0083\end{eqnarray}$$
「正規標準分布の上側確率」の表から、
$P(Z>2.39)=0.0084,\ P(Z>2.40)=0.0082$
$\therefore\ z/\sqrt{2}=2.395\ \Rightarrow\ z=2.395\times\sqrt{2}=3.387$


問21 解答

(一元配置分散分析)

[1]

$\boxed{ \ \mathsf{30}\ }$ ③

対象とするパソコン、調査する対策、計測の順序をランダムに決めているので③が最も適切である。
(①②はパソコンの購入時期の影響を受ける。④⑤は固有のパソコンの性能と対策の順番の影響を受ける。)

[2]

$\boxed{ \ \mathsf{31}\ }$ ②

対策の平方和(水準間平方和)の自由度は対策(水準)の数$-1$なので、$3-1=2$。
誤差の平方和の自由度は総データ数$-$水準の数なので、$12-3=9$。

[3]

$\boxed{ \ \mathsf{32}\ }$ ⑤

水準数$a$、総観測値数$n$の一元配置分散分析において、水準$j$の標本平均及び観測値数を$y_{j\cdot}, n_j$、残差平方和を$S_e$、残差の自由度を$\phi_e$、残差の平均平方を$V_e$とする。
水準$j$の母平均の$100(1-\alpha)\%$信頼区間は$t_{0.05/2}(12-3)=t_{0.025}(9)=2.262$
$$y_{j\cdot}\pm t_{\alpha/2}(\phi_e)\sqrt{\frac{V_e}{n_j}}=y_{j\cdot}\pm t_{\alpha/2}(n-a)\sqrt{\frac{S_e}{(n-a)n_j}}$$
対策$3$の効果の点推定値が$-49.9$なのでこの効果の$95\%$信頼区間
$$\mu-49.9\pm 2.262\times\sqrt{\frac{1890.1}{(12-3)\times4}}=\mu-49.9\pm16.39$$
効果の信頼区間は$[-66.29,-33.51]$となる。


問22 解答

(重回帰モデル,統計ソフトウェアの活用)

※重回帰モデルの統計ソフトウェアによる出力結果の主な項目
$\mathtt{Estimate}$:回帰係数の推定値
$\mathtt{Std.Error}$:回帰係数の推定値の標準誤差
$\mathtt{t\ value}$:$t$値,$\mathtt{Pr(\gt|t|)}$:$P_-$値・・・回帰係数の検定で使う
$\mathtt{Rasidual\ standard\ error}$:誤差項の標準偏差の推定値
$\mathtt{degrees\ of\ freedom}$:自由度
$\mathtt{Multiple\ R-squared}$:決定係数($R^2$)
$\mathtt{Adjusted\ R-squared}$:自由度調整済み決定係数($R^{*2}$)
$\mathtt{F-statistic}$:$F$検定統計量,$\mathtt{p-value}$:$P_-$値・・・回帰の有意性の検定で使う

[1]

$\boxed{ \ \mathsf{33}\ }$ ⑤

$\mathtt{t\ value}$はある説明変数$x_j$は被説明変数$y$の予測に役立たない$iff\ H_0:\beta_j=0$($beta_j$は説明変数$x_j$の回帰係数)とする帰無仮説のもとで、$\hat\beta_j$に基づく$t$統計量の実現値である。回帰係数の推定値の標準誤差を$se(\hat\beta_j)$とすると、$t$統計量の実現値は、$$t=\hat\beta_j/se(\hat\beta_j)\sim t(n-p-1)$$である($n-p-1$は残差の自由度)。
したがって、(ア)の値は$-9.614/3.575=-2.689$

[2]

$\boxed{ \ \mathsf{34}\ }$ ②

① $P_-$値が最も小さい説明変数は$\log($人口密度$)$である。誤り。
② 政令指定都市ダミーの回帰係数が$-0.198$であるので、政令指定都市であれば$1$人あたり社会体育施設数は$e^-0.198=0.82$倍($2$割減)となる。正しい。
③ $15$歳未満人口の割合の$P_-$値は$0.333\ge10\%$なので、帰無仮説は棄却できない。誤り。
④ $log(1$人当たり所得$)$の回帰係数が正なので、$1$人当たり所得が低ければ、$1$人あたり社会体育施設数は少なくなる傾向にある。誤り。
⑤ 統計的に優位性が確認されたとしても、それは説明変数と被説明変数の間に相関関係が見られるということであって、説明変数の被説明変数への因果関係が存在するとは必ずしも言えない。誤り。

[3]

$\boxed{ \ \mathsf{35}\ }$ ④

Ⅰ. $log(1$人当たり所得$)$は、モデルAの結果において有意水準$5\%$で有意でない($P_-$値$\gt0.05$)。誤り。
Ⅱ. 2つのモデルの自由度調整済み決定係数($\mathtt{Adjusted\ R-squared}$)の値を比較すると、モデルBのほうが値が高いので、モデルBのほうがより良いモデルである。正しい。
Ⅲ. 2つのモデルの$F$検定の$P_-$値を見ると、それぞれ$4.494\times10^{-16}, 2.2\times10^{-16}$未満であり、極めて小さい値であることから、説明変数にかかるすべての係数がゼロであるという帰無仮説は棄却される。正しい。



統計検定2級問題解説 ~2021年6月実施~ (問1~問12)

過去問題

過去問題は統計検定公式問題集が問題と解答例を公開しています。こちらを参照してください。


問1 解答

(歪度)

$\boxed{ \ \mathsf{1}\ }$ ①

歪度は,右に裾が長い分布では正の値になり,左に裾が長い分布では負の値になる。
問題のグラフは全体的に右に裾の長い分布となっている。


問2 解答

(年次変化率,幾何平均)

$\boxed{ \ \mathsf{2}\ }$ ②

時点(年,月など) $t$の観測値を$y_t$としたとき,
年次変化率 $(y_{t+1}-y_t)/y_t$ または $y_{t+1}/y_t-1$
$t$年$(1950\le{t}\le1954)$の観測値を$y_t$,年次変化率を$r$としたとき,変化率の平均は幾何平均となるので、
$$
\begin{align} 
r&=(y_{1951}/y_{1950}\times y_{1952}/y_{1951}\times y_{1953}/y_{1952}\times y_{1954}/y_{1953})^{1/4}-1\\
&=(y_{1954}/y_{1950})^{1/4}-1\fallingdotseq 0.154=15.4\%
\end{align}
$$
上式から明らかなことは,各年の変化率の幾何平均は
最初年の値を$y_0$,最後年の値を$y_t$,期間を$t$とすると $\left({y_t}/{y_0}\right)^{1/t}$ で求められる。


問3 解答

(価格指数)

$\boxed{ \ \mathsf{3}\ }$ ②

主な物価指数(デフレータ
・ラスパイレス物価指数…基準年の購入量や取引量等を重みとして算出した価格指数。
品目$i$の基準年価格$=p_{oi}$,基準年数量$=q_{oi}$,比較年価格$=p_{ti}$として$$P_L=\frac{\sum_ip_{ti}q_{0i}}{\sum_ip_{0i}q_{0i}}$$
・パーシェ物価指数…比較年の購入量や取引量等を重みとして算出した価格指数。
品目$i$の基準年価格$=p_{oi}$,比較年数量$=q_{ti}$,比較年価格$=p_{ti}$として$$P_P=\frac{\sum_ip_{ti}q_{ti}}{\sum_ip_{0i}q_{ti}}$$
・フィッシャー物価指数…ラスパイレス指数とパーシェ指数の幾何平均で求められる価格指数。$$P_F=\sqrt{P_L\times P_P}$$

問題は比較年のパーシェ指数を求めるので(基準年を$100$とする)、$$P_P=\frac{80\times80+90\times70}{78\times80+84\times70}\times100=104.8$$


問4 解答

(時系列データの指数化)

[1]

$\boxed{ \ \mathsf{4}\ }$ ③

東京都の$2019$年の新聞発行部数の指数はグラフから$57$
したがって、$1990$年から$2019$年にかけての新聞発行部数は、$1990$年を$100$として$100-57=43$減少したことになる。よって、減少部数は
 $5,190$万部$\times 13\%\times 43/100\fallingdotseq 290$万部

[2]

$\boxed{ \ \mathsf{5}\ }$ ①

問題にある前年比増加率の定義から、前年から減少している場合、増加率は負の値となる。
発行部数指数のグラフを見ると、$2005$年以降は東京都、鳥取県ともに毎年減少している。このことから、前年比増加率のグラフも$2005$年以降は両都県ともに負の値となっている必要がある。これを満たすグラフは①のみである。
(②は鳥取県の$2013$年、④は東京都の$2008$年、⑤は東京都の$2006$年でそれぞれ正の値をとっている。)


問5 解答

(散布図,相関係数,共分散)

[1]

$\boxed{ \ \mathsf{6}\ }$ ⑤

Ⅰ.散布図から、飛型点のとりうる範囲に比べ、飛距離点の取りうる範囲のほうが大きいく、ばらつきが大きいので、飛距離点のほうが分散が大きくなる。正しい。
Ⅱ.2つの散布図を比較すると、飛距離点と飛型点の散布図のほうがプロットされた点が直線状に分布しているので相関が強く、相関係数も高くなる。正しい。
Ⅲ.飛距離点と飛型点の散布図に回帰直線をあてはめると、明らかに$y$切片は正の値をとる。正しい。

[2]

$\boxed{ \ \mathsf{7}\ }$ ④

飛距離点の定義から、飛距離$x$は飛距離点$y$の$1/2$に比例する($x=y/2+68$)。
このとき、飛距離の分散は飛距離点の分散の$1/4$,標準偏差は$1/2$になる。
共分散は,2つのデータの平均からの偏差の積の和なので,片方のデータの平均からの偏差が$1/2$倍となれば,$1/2$倍となる。
相関係数は、共分散をそれぞれの標準偏差で割ったものだから、飛距離と飛型点の相関係数は飛距離点と飛型点の相関係数と等しくなる。


問6 解答

(相関係数,共分散)

$\boxed{ \ \mathsf{8}\ }$ ③

共分散$$\begin{align}s_{xy}&=\frac1{10}\sum_{i=1}^{10}(x_i-\bar{x})(y_i-\bar{y})\\&=\frac1{10}\sum_{i=1}^{10}x_iy_i-\frac1{10}\sum_{i=1}^{10}x_i\bar{y}-\frac1{10}\sum_{i=1}^{10}\bar{x}y_i+\frac1{10}\sum_{i=1}^{10}\bar{x}\bar{y}\\&=\frac1{10}\sum_{i=1}^{10}x_iy_i-\bar{x}\bar{y}\\&=\frac1{10}\sum_{i=1}^{10}x_iy_i-\frac1{10}\sum_{i=1}^{10}x_i\times\frac1{10}\sum_{i=1}^{10}y_i\\&= \frac1{10} \times4548.7-\frac1{10}\times346.3\times\frac1{10}\times121.8=33.08\end{align}$$
標準偏差$$\sigma_x=\sqrt{\frac{10-1}{10}s_x^2}=\sqrt{\frac{9}{10}\times167.4}=12.27$$$$\sigma_y=\sqrt{\frac{10-1}{10}s_y^2}=\sqrt{\frac{9}{10}\times11.6}=3.23$$
よって、相関係数は$$r=\frac{s_{xy}}{\sigma_x\sigma_y}=\frac{33.08}{12.27\times3.23}=0.83$$


問7 解答

(標本抽出法)

$\boxed{ \ \mathsf{9}\ }$ ④

多段抽出法・・・母集団をいくつかのグループ(第1段抽出単位)に分け,そこから無作為抽出でいくつかグループを選び,さらにその中から無作為抽出でいくつかのグループ(第2段抽出単位)を選び・・を何段か繰り返してそこから標本を無作為に抽出する。
層化抽出法・・・母集団をあらかじめいくつかの層(グループ)に分けておき、各層の中から必要な数の調査対象を無作為に抽出する方法。
集落(クラスター)抽出法・・・母集団を小集団であるクラスター(集落)に分け,その中からいくつかのクラスターを無作為に抽出し,それぞれのクラスターにおいて全数調査を行う。


問8 解答

(同時確率関数と相関係数)

$\boxed{ \ \mathsf{10}\ }$ ⑤

まず、$X, Y$の周辺分布を求める。
$\begin{align}
P(X=-1)=P(X=1)&=0+1/4+0=1/4\\
P(X=0)&=1/4+0+1/4=1/2\\
P(Y=-1)=P(Y=1)&=0+1/4+0=1/4\\
P(Y=0)&=1/4+0+1/4=1/2
\end{align}$
これらから$X^2,Y^2$の期待値$E[X^2],E[Y^2]$、分散$V[X^2],V[Y^2]$及び$X^2$と$Y^2$の共分散$Cov(X^2,Y^2)$、相関係数$r$を求める。
$$\begin{align}
E[X^2]&= (-1)^2\times P(X=-1)+0^2\times P(X=0)+1^2\times P(X=1)\\
&=1/4+0+1/4=1/2=\mu_{X^2}\\
E[Y^2]&= (-1)^2\times P(Y=-1)+0^2\times P(Y=0)+1^2\times P(Y=1)\\
&=1/4+0+1/4=1/2=\mu_{Y^2}\\
V[X^2]&=E[(X^2-\mu_{X^2})^2]\\
&=\{(-1)^2-1/2\}^2\times P(X=-1)+(0^2-1/2)^2\times P(X=0)+(1^2-1/2)^2\times P(X=1)\\
&=1/16+1/8+1/16=1/4\\
V[Y^2]&=E[(Y^2-\mu_{Y^2})^2]\\
&={(-1)^2-1/2}^2\times P(Y=-1)+(0^2-1/2)^2\times P(Y=0)+(1^2-1/2)^2\times P(Y=1)\\
&=1/16+1/8+1/16=1/4\\
Cov[X^2,Y^2]&=E[(X^2-\mu_{X^2})(Y^2-\mu_{Y^2})\\
&={(-1)^2-1/2}\times{(-1)^2-1/2}\times P(X=-1,Y=-1)\\
&+{(-1)^2-1/2}\times(0^2-1/2)\times P(X=-1,Y=0)\\
&+{(-1)^2-1/2}\times(1^2-1/2)\times P(X=-1,Y=1)\\
&+(0^2-1/2)\times{(-1)^2-1/2}\times P(X=0,Y=-1)\\
&+(0^2-1/2)\times(0^2-1/2)\times P(X=0,Y=0)\\
&+(0^2-1/2)\times(1^2-1/2)\times P(X=0,Y=1)\\
&+(1^2-1/2)\times{(-1)^2-1/2}\times P(X=1,Y=-1)\\
&+(1^2-1/2)\times(0^2-1/2)\times P(X=1,Y=0)\\
&+(1^2-1/2)\times(1^2-1/2)\times P(X=1,Y=1)\\
&=0-1/16+0-1/16+0-1/16+0-1/16+0=-1/4\\
\therefore r&=\frac{Cov[X^2,Y^2]}{\sqrt{V[X^2]V[Y^2]}}=\frac{-1/4}{\sqrt{1/4\times1/4}}=-1
\end{align}$$
また
$P(X^2=0,Y^2=0)=P(X=0,Y=0)=0$
$P(X^2=0)=P(X=0)=1/2, P(Y^2=0)=P(Y=0)=1/2$
であることから
$P(X^2=0,Y^2=0)\neq P(X^2=0)P(Y^2=0)$
となり、$P(X^2\cap Y^2)=P(X^2)P(Y^2)$が成り立たないため、$X^2$と$Y^2$は互いに独立ではない。


問9 解答

(非復元抽出の確率)

$\boxed{ \ \mathsf{11}\ }$ ⑤

無作為に集められた$25$人の中に同じ誕生日の人が存在する確率を求めるためには、同じ誕生日の人が全くいない確率を求めて$1$から引けばよい。
同じ誕生日がないということなので、$365$日から重複を許さずに$25$日を抽出する確率を求める。
$$\underbrace{\frac{365}{365}\times\frac{364}{365}\times\frac{363}{365}\times\cdots\times\frac{341}{365}}_{25}=\frac1{365^{25}}\times\frac{365!}{340!}$$
よって、同じ誕生日の人が存在する確率は
$$1-\frac{365!}{365^{25}\times340!}$$


問10 解答

(正規分布、標準正規分布)

$\boxed{ \ \mathsf{12}\ }$ ④

確率変数$X$が正規分布$N(60,9^2)$に従うとき、$$Z=\frac{X-60}{9}$$は標準正規分布$N(0,1)$に従う。そこで、
$$P(X\le c)=0.011\iff P(Z\le \frac{c-60}{9})=0.011$$
「標準正規分布の上側確率」の表から$P(X\ge2.29)=0.011$なので、$P(X\le-2.29)=0.011$、$$\begin{align}\therefore \frac{c-60}{9}&=-2.29\\c&=39.39\end{align}$$


問11 解答

(連続型確率変数)

[1]

$\boxed{ \ \mathsf{13}\ }$ ①

$P(X\gt 1)=1-P(X\le 1)=1-F(1)=1-1=0$

[2]

$\boxed{ \ \mathsf{14}\ }$ ③

確率密度関数$f(x)$は累積分布関数$F(x)$を微分して求める。
$$\begin{eqnarray}
f(x)=\frac{d}{dx}F(x)=
\begin{cases}
1&(0\ge x\lt 1)\\
0&(x\lt 0,\ 1\ge x)
\end{cases}
\end{eqnarray}$$
期待値$E(X)$は、
$$E(X)=\int_{-\infty}^\infty xf(x)dx=\int_0^1 x\cdot1dx=\left[\frac{1}{2}x^2\right]_0^1=\frac{1}{2}$$


問12 解答

(幾何分布、チェビシェフの不等式)

[1]

$\boxed{ \ \mathsf{15}\ }$ ③

幾何分布$P(X=x)=p(1-p)^{x-1}$の期待値(平均)は$1/p$
$$\therefore\ P(X)=\frac13\left(\frac23\right)^{n-1}\ \Rightarrow\ E(X)=\frac1{1/3}=3$$

※成功か失敗しかない試行をベルヌーイ試行という。成功確率は $p$。
このベルヌーイ試行を独立に何回も行うとき,初めて成功するまでに“試行”した回数を $X$ とすると,$X$ の確率関数は$$P(X=x)=p(1-p)^{x-1}$$となり,この確率分布をパラメータ $p$ の幾何分布という。(本によっては,初めて成功するまでに“失敗”した回数を $X$ とする定義の仕方もある。)
ここで,等比級数の和$$\displaystyle \sum_{x=0}^\infty a^x=\frac1{1-a}\ \ \ (|a|<1)$$の両辺を $a$ で微分すると$$\displaystyle \sum_{x=0}^\infty xa^{x-1}=\frac1{(1-a)^2}$$さらに,この式の両辺を $a$ で微分すると,$$\displaystyle \sum_{x=0}^\infty x(x-1)a^{x-2}=\frac2{(1-a)^3}$$となる。これを利用して,幾何分布の期待値と分散を求める。$$\begin{align}E[X]=&\sum_{x=0}^\infty xp(1-p)^{x-1}=p\sum_{x=0}^\infty x(1-p)^{x-1}\\=&\frac{p}{\{1-(1-p)\}^2}=\frac1p\\V[X]=&E[X(X-1)]+E[X]-E[X]^2\\=&\sum_{x=0}^\infty x(x-1)p(1-p)^{x-1}+\frac1p-\frac1{p^2}\\=&p(1-p)\sum_{x=0}^\infty x(x-1)(1-p)^{x-2}+\frac1p-\frac1{p^2}\\=&\frac{2p(1-p)}{\{1-(1-p)\}^3}+\frac1p-\frac1{p^2}\\=&\frac{2-2p}{p^2}+\frac{p}{p^2}-\frac1{p^2}=\frac{1-p}{p^2}\\\end{align}$$

[2]

$\boxed{ \ \mathsf{16}\ }$ ④

期待値$E[X]$、分散$V[X]$を持つ確率分布に従う確率変数$X$について、任意の$\epsilon\gt 0$に対して、チェビシェフの不等式$$P(|X-E[X]|\ge\epsilon)\le V[X]/\epsilon^2$$が成り立つ。
また、母平均$\mu=3$、母分散$\sigma^2=6$の母集団から抽出した標本$X_1,\cdots,X_n$の標本平均$\displaystyle\bar{X}=\frac1n\sum_{i=1}^nX_i$の期待値と分散は
$$\begin{eqnarray}E[\bar{X}]=E\left[\frac1n\sum_{i=1}^nX_i\right]=\frac1n\sum_{i=1}^nE[X_i]=\frac{n\mu}n=\mu=3\\V[\bar{X}]=V\left[\frac1n\sum_{i=1}^nX_i\right]=\frac1{n^2}\sum_{i=1}^nV[X_i]=\frac{n\sigma^2}{n^2}=\frac{\sigma^2}{n}=\frac{6}{n}\end{eqnarray}$$
これをチェビシェフの不等式に当てはめると
$$P(|X-3|\ge\epsilon)\le\frac{6/n}{\epsilon^2}$$


推測統計フローチャート(推定、検定を考えるにあたっての解法の整理)

中心極限定理などに基づいて母集団の確率分布のパラメータの点推定・区間推定や、パラメータに関する仮説の検定を行う推測統計は、基本的な考え方は一貫している一方で推定の対象や分散の既知・未知などに置ける場合分けなど、関連する概念が多くわかりにくい。
そのため当稿では解法の整理の補助となるように、推測統計に関連するトピックをフローチャートの形式にまとめる。作成にあたっては、「基礎統計学Ⅰ 統計学入門(東京大学出版会)」の$9$章〜$12$章を主に参考にした。

大枠の整理

推測統計を考える際の前提

推測統計を考える際に前提となるのが母集団(population)と標本(sample)である。記述統計学(descriptive statistics)では得られた標本についてのみ考えるが、得られた標本の裏側の母集団についても考察を行うのが推測統計である。

推測統計では母集団の持つ分布である、母集団分布(population distribution)について知ることが目的となる。母集団分布を考えるにあたっては、事象に対して基本的には理論的・経験的に正規分布やポアソン分布などの確率分布をあてはめることができるとされることが多い。正規分布やポアソン分布を母集団分布に仮定できる場合、正規分布の$\mu$や$\sigma^2$、二項分布の$p$、ポアソン分布の$\lambda$のようにいくつかの確率分布のパラメータさえわかれば母集団分布について全て知ることができる。統計的推測(statistical inference)ではこのパラメータを母数(parameter)と呼び、母数の推測が推測統計の目的となる。

よって、推測統計では母集団分布が正規分布の際の$\mu$や$\sigma^2$などが推測の対象となる。

フローチャート

推測統計の手法を用いる際のフローチャートは概ね上記のようになる。書籍などで確認するとトピックが多く大変だが、このように図で整理することでパターンの把握が容易になるのではと思われる。
どの問題もまず最初に推測の対象を考えるとよく、基本的には「母平均」、「母分散」、「その他」で把握しておくと良い。

以下、フローチャートに対応させつつそれぞれのパターンについて確認する。区間推定、検定のどちらも考え方自体はそれほど変わらないため、ここでは区間推定を元に考えることとする。また、全てを同時に抑えるのは大変なので、発展項目と思われるものについては*をつけた。

上側確率について

区間推定や検定を行う際に、確率密度関数の変数$x$の位置とその位置における積分値の対応を考える際に、$\alpha$を導入して考えることが多い。たとえば区間推定を行う際に、母集団分布の母数$\theta$が$L \leq \theta \leq U$のように$L$(lower confidence limit)、$U$(upper confidence limit)を用いて表せると考えると、確率分布と$\alpha$、$L$、$U$の関係は下記のようになる。
$$
\begin{align}
P(L \leq \theta \leq U) = 1 – \alpha
\end{align}
$$
上記が基本的な$\alpha$の考え方であり、標準正規分布表や$t$分布表、$\chi^2$分布の表などと見比べて区間推定や検定を行う。一方で、この際に、表の読み取りの際に$\alpha/2$などが出てきてわかりにくくなるケースがある。

このような状況を防ぐにあたって、「$XX$分布において上側確率が$100\alpha$%となるパーセント点に対応する$XX$の値を$XX_{\alpha}$とする」と表記し、$5$%区間を考えるにあたっては上側確率が2.5%の際は$XX_{\alpha=0.025}$のように、$\alpha/2$ではなく$\alpha=0.025$のように数値を変更する形式で表す方が取り扱いやすいと思われる。また、$x=0$を中心とする左右対称の確率分布の場合は$XX_{\alpha=0.975}=-XX_{\alpha=0.025}<0$であることも抑えておくと良い。これは標準正規分布と$t$分布に当てはまる。

母平均の区間推定(母分散既知)

標本数を$n$、標本平均を$\bar{x}$、母平均を$\mu$、母分散を$\sigma^2$とする。
$$
\begin{align}
z = \frac{\bar{x}-\mu}{\sigma/\sqrt{n}}
\end{align}
$$
上記のように$z$を計算したとき、$z$は平均を引いたのちに標準偏差で割っているので標準化されたと考えることができる。よって$z$は標準正規分布$N(0,1)$に従うと考えることができる。

これにより、$z$値を標準正規分布と見比べることで区間推定を行うことができる。母平均$\mu$の$95$%区間を考えるにあたって、対応する標準正規分布の区間が$z_{\alpha=0.975} \leq z \leq z_{\alpha=0.025}$のように表せるとする。ここで正規分布表より、$z_{\alpha=0.975}=-1.96$、$z_{\alpha=0.025}=1.96$を満たすことが読み取れる。このとき、下記の式変形によって$\mu$の$95$%区間を求めることができる。
$$
\begin{align}
z_{\alpha=0.975} \leq &z \leq z_{\alpha=0.025} \\
z_{\alpha=0.975} \leq &\frac{\bar{x}-\mu}{\sigma/\sqrt{n}} \leq z_{\alpha=0.025} \\
z_{\alpha=0.975}\frac{\sigma}{\sqrt{n}} \leq &\bar{x}-\mu \leq z_{\alpha=0.025}\frac{\sigma}{\sqrt{n}} \\
-z_{\alpha=0.025}\frac{\sigma}{\sqrt{n}} \leq &\mu-\bar{x} \leq -z_{\alpha=0.975}\frac{\sigma}{\sqrt{n}} \\
\bar{x}-z_{\alpha=0.025}\frac{\sigma}{\sqrt{n}} \leq &\mu \leq \bar{x}-z_{\alpha=0.975}\frac{\sigma}{\sqrt{n}} \\
\bar{x}-1.96\frac{\sigma}{\sqrt{n}} \leq &\mu \leq \bar{x}+1.96\frac{\sigma}{\sqrt{n}}
\end{align}
$$

母平均の区間推定(母分散未知)

標本数を$n$、標本平均を$\bar{x}$、母平均を$\mu$、標本分散を$s^2$とする。
$$
\begin{align}
t = \frac{\bar{x}-\mu}{s/\sqrt{n}}
\end{align}
$$
上記のように$t$を計算したとき、$t$はt分布$t(n-1)$に従うと考えることができる。

これにより、$t$値をt分布$t(n-1)$と見比べることで区間推定を行うことができる。母平均$\mu$の95%区間を考えるにあたって、対応する$t$分布の区間が$t_{\alpha=0.975}(n-1) \leq t \leq t_{\alpha=0.025}(n-1)$のように表せるとする。ここで$n=10$の際はt分布の表より、$t_{\alpha=0.025}(10-1)=t_{\alpha=0.025}(9)=2.262$、$t_{\alpha=0.975}(10-1)=-t_{\alpha=0.025}(9)=-2.262$を満たすことが読み取れる。このとき、下記の式変形によって$\mu$の95%区間を求めることができる。
$$
\begin{align}
t_{\alpha=0.975}(9) \leq &t \leq t_{\alpha=0.025}(9) \\
t_{\alpha=0.975}(9) \leq &\frac{\bar{x}-\mu}{s/\sqrt{n}} \leq t_{\alpha=0.025}(9) \\
t_{\alpha=0.975}(9)\frac{s}{\sqrt{n}} \leq &\bar{x}-\mu \leq t_{\alpha=0.025}(9)\frac{s}{\sqrt{n}} \\
-t_{\alpha=0.025}(9)\frac{s}{\sqrt{n}} \leq &\mu-\bar{x} \leq -t_{\alpha=0.975}(9)\frac{s}{\sqrt{n}} \\
\bar{x}-t_{\alpha=0.025}(9)\frac{s}{\sqrt{n}} \leq &\mu \leq \bar{x}+t_{\alpha=0.025}(9)\frac{s}{\sqrt{n}} \\
\bar{x}-2.262\frac{s}{\sqrt{n}} \leq &\mu \leq \bar{x}+2.262\frac{s}{\sqrt{n}}
\end{align}
$$
サンプル数が多くなるにつれて$t$分布は正規分布に近づく。一方で、サンプル数が少ない際は母分散既知の場合に比較して$95$%区間は大きくなる。このことは「母分散がわからない方が母集団の不確実性が大きい」と直感的に解釈しておくと良いと思われる。

母分散の区間推定

標本数を$n$、母分散を$\sigma^2$、標本分散を$s^2$とする。
$$
\begin{align}
\chi^2 = \frac{(n-1)s^2}{\sigma^2}
\end{align}
$$
上記のように$\chi^2$を計算したとき、$t$は$\chi^2$分布$\chi^2(n-1)$に従うと考えることができる。これにより、$\chi^2$値を$\chi^2$分布$\chi^2(n-1)$と見比べることで区間推定を行うことができる。

ここで$\chi^2$分布において上側確率が$100\alpha$%となるパーセント点に対応する$\chi^2$の値を$\chi^2_{\alpha}$とする。このとき母分散$\sigma^2$の$95$%区間は、$\chi^2_{\alpha=0.975}(n-1) \leq \chi^2 \leq \chi^2_{\alpha=0.025}(n-1)$のように表せる。ここで$n=10$の際は$\chi^2$分布の表より、$\chi^2_{\alpha=0.975}(10-1)=\chi^2_{\alpha=0.975}(9)=2.70039$、$\chi^2_{\alpha=0.025}(10-1)=\chi^2_{\alpha=0.025}(9)=19.0228$を満たすことが読み取れる。このとき、下記の式変形によって$\sigma^2$の$95$%区間を求めることができる。
$$
\begin{align}
\chi^2_{\alpha=0.975}(9) \leq &\chi^2 \leq \chi^2_{\alpha=0.025}(9) \\
\chi^2_{\alpha=0.975}(9) \leq &\frac{(n-1)s^2}{\sigma^2} \leq \chi^2_{\alpha=0.025}(9) \\
\frac{1}{\chi^2_{\alpha=0.025}(9)} \leq &\frac{\sigma^2}{(n-1)s^2} \leq \frac{1}{\chi^2_{\alpha=0.975}(9)} \\
\frac{(n-1)s^2}{\chi^2_{\alpha=0.025}(9)} \leq &\sigma^2 \leq \frac{(n-1)s^2}{\chi^2_{\alpha=0.975}(9)} \\
\frac{9s^2}{19.0228} \leq &\sigma^2 \leq \frac{9s^2}{2.70039}
\end{align}
$$

母平均の差の区間推定(母分散が未知だが等しいと仮定)

母集団分布を$N(\mu_1, \sigma^2)$と$N(\mu_2, \sigma^2)$で表すことのできる母分散の等しい二つの正規母集団から個別に標本$X_1, X_2, …, X_m$と$Y_1, Y_2, …, Y_n$を抽出した際の母平均の差$\mu_1-\mu_2$の区間推定について考える。このとき不偏分散を$s^2$とすると$s^2$は下記のように計算できる。
$$
\begin{align}
s^2 = \frac{1}{m+n-2} \left( \sum_{j=1}^{m}(X_i-\bar{X})^2 + \sum_{j=1}^{n}(Y_i-\bar{Y})^2 \right)
\end{align}
$$

ここで下記のように$t$値を計算する。
$$
\begin{align}
t = \frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{s \sqrt{1/m+1/n}}
\end{align}
$$

このとき上記は自由度$m+n-2$の$t$分布$t(m+n-2)$に従うため、$t(m+n-2)$見比べることで区間推定を行うことができる。母平均の差$\mu_1-\mu_2$の$95$%区間を考えるにあたって、対応する$t$分布$t(m+n-2)$の区間が$t_{\alpha=0.975}(m+n-2) \leq t \leq t_{\alpha=0.025}(m+n-2)$のように表せるとする。ここで$m=10$、$n=10$の際は$t$分布の表より、$t_{\alpha=0.975}(10+10-2)=-t_{\alpha=0.025}(18)=-2.101$、$t_{\alpha=0.025}(10+10-2)=t_{\alpha=0.025}(18)=2.101$を満たすことが読み取れる。このとき、下記の式変形によって$\mu_1-\mu_2$の$95$%区間を求めることができる。
$$
\begin{align}
t_{\alpha=0.975}(18) \leq &t \leq t_{\alpha=0.025}(18) \\
t_{\alpha=0.975}(18) \leq &\frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{s \sqrt{1/m+1/n}} \leq t_{\alpha=0.025}(18) \\
-t_{\alpha=0.025}(18) \leq &\frac{(\mu_1-\mu_2)-(\bar{X}-\bar{Y})}{s \sqrt{1/m+1/n}} \leq -t_{\alpha=0.975}(18) \\
-t_{\alpha=0.025}(18)s \sqrt{1/m+1/n} \leq &(\mu_1-\mu_2)-(\bar{X}-\bar{Y}) \leq -t_{\alpha=0.975}(18)s \sqrt{1/m+1/n} \\
(\bar{X}-\bar{Y})-t_{\alpha=0.025}(18)s \sqrt{1/m+1/n} \leq &(\mu_1-\mu_2) \leq (\bar{X}-\bar{Y})+t_{\alpha=0.025}(18)s \sqrt{1/m+1/n} \\
(\bar{X}-\bar{Y})-2.101s \sqrt{1/10+1/10} \leq &(\mu_1-\mu_2) \leq (\bar{X}-\bar{Y})+2.101s \sqrt{1/10+1/10} \\
(\bar{X}-\bar{Y})-\frac{2.101s}{\sqrt{5}} \leq &(\mu_1-\mu_2) \leq (\bar{X}-\bar{Y})+\frac{2.101s}{\sqrt{5}}
\end{align}
$$

母平均の差の区間推定(母分散が未知であるかつ等しいと仮定できない)*

二つの母分散が等しいと仮定できない場合の母平均の差の場合は、下記のように集団ごとに不偏分散を計算する。
$$
\begin{align}
s_1^2 &= \frac{1}{m-1} \sum_{j=1}^{m}(X_i-\bar{X})^2 \\
s_2^2 &= \frac{1}{n-1} \sum_{j=1}^{n}(Y_i-\bar{Y})^2
\end{align}
$$
上記のようにそれぞれの集団の不偏分散を計算したのちにウェルチの近似法を用いて$t$値を計算する。
$$
\begin{align}
t = \frac{(\bar{X}-\bar{Y}) – (\mu_1-\mu_2)}{\sqrt{s_1^2/m + s_2^2/n}}
\end{align}
$$

このとき下記のように$\nu$を計算する。
$$
\begin{align}
\nu = \frac{(s_1^2/m + s_2^2/n)^2}{\frac{(s_1^2/m)^2}{m-1} + \frac{(s_2^2/n)^2}{n-1}}
\end{align}
$$
ここで$\nu$に一番近い整数を$\nu’$とすると前述のように計算した$t$値は自由度$\nu’$の$t$分布$t(\nu’)$に近似的に従う。

このとき$\mu_1-\mu_2$の$95$%区間はこれまでと同様に求めることができる。
$$
\begin{align}
t_{\alpha=0.975}(\nu’) \leq &t \leq t_{\alpha=0.025}(\nu’) \\
(\bar{X}-\bar{Y})-t_{\alpha=0.025}(\nu’)\sqrt{s_1^2/m + s_2^2/n} \leq &(\mu_1-\mu_2) \leq (\bar{X}-\bar{Y})+t_{\alpha=0.025}(\nu’)\sqrt{s_1^2/m + s_2^2/n}
\end{align}
$$

他と比較しても式が複雑なので途中計算は省略した。

母分散の比の区間推定 *

二つの集団のサンプル数を$m$、$n$、母分散を$\sigma_1^2$、$\sigma_2^2$とし、不偏標本分散を$s_1^2$、$s_2^2$とする。
$$
\begin{align}
F = \frac{\sigma_2^2 s_1^2}{\sigma_1^2 s_2^2}
\end{align}
$$

このとき上記のように$F$値を定義すると$F$値は自由度$(m-1,n-1)$の$F$分布$F(m-1,n-1)$に従う。この導出にあたっての詳細は省略したが「基礎統計学Ⅰ 統計学入門(赤本)」の$10.5.2$節の記載が参考になる。

さて、このとき$\displaystyle \frac{\sigma_2^2}{\sigma_1^2}$の$95$%区間を求める。
$$
\begin{align}
F_{\alpha=0.975}(m-1,n-1) \leq &F \leq F_{\alpha=0.025}(m-1,n-1) \\
F_{\alpha=0.975}(m-1,n-1) \leq &\frac{\sigma_2^2 s_1^2}{\sigma_1^2 s_2^2} \leq F_{\alpha=0.025}(m-1,n-1) \\
F_{\alpha=0.975}(m-1,n-1)\frac{s_2^2}{s_1^2} \leq &\frac{\sigma_2^2}{\sigma_1^2} \leq F_{\alpha=0.025}(m-1,n-1)\frac{s_2^2}{s_1^2}
\end{align}
$$

負の値や逆数を取るなどの不等号が反転する演算を行っていないことに注意しておくとよい。不等号の反転がある場合の計算がややこしい場合は不等号を一つずつ計算する方がミスが減らせるので、わからなくなったら一つずつ計算する方が良いと思われる。

母比率の区間推定(二項分布)

サンプル数$n$、母比率$p$の二項分布$Binom(n,p)$に従う確率変数$X$の母平均は$np$、母分散は$np(1-p)$と表すことができる。
https://www.hello-statisticians.com/explain-books-cat/toukeigakunyuumon-akahon/ch6_practice.html#61
$$
\begin{align}
z = \frac{X-np}{\sqrt{np(1-p)}}
\end{align}
$$

このとき中心極限定理に基づいて、上記が従う分布は標準正規分布で近似できる。よって、$z$値を標準正規分布と見比べることで区間推定を行うことができる。母比率$p$の$95$%区間を考えるにあたって、対応する標準正規分布の区間が$z_{\alpha=0.975} \leq z \leq z_{\alpha=0.025}$のように表せるとする。ここで正規分布表より、$z_{\alpha=0.975}=-1.96$、$z_{\alpha=0.025}=1.96$を満たすことが読み取れる。このとき、下記の式変形によって$p$の$95$%区間を求めることができる。
$$
\begin{align}
z_{\alpha=0.975} \leq &z \leq z_{\alpha=0.025} \\
z_{\alpha=0.975} \leq &\frac{X-np}{\sqrt{np(1-p)}} \leq z_{\alpha=0.025} \\
-z_{\alpha=0.025} \leq &\frac{p-X/n}{\sqrt{p(1-p)/n}} \leq -z_{\alpha=0.975} \\
-z_{\alpha=0.025} \leq &\frac{p-\hat{p}}{\sqrt{p(1-p)/n}} \leq -z_{\alpha=0.975} \\
-z_{\alpha=0.025}\sqrt{p(1-p)/n} \leq &p-\hat{p} \leq -z_{\alpha=0.975}\sqrt{p(1-p)/n} \\
\hat{p}-z_{\alpha=0.025}\sqrt{p(1-p)/n} \leq &p \leq \hat{p}+z_{\alpha=0.025}\sqrt{p(1-p)/n} \\
\hat{p}-1.96\sqrt{p(1-p)/n} \leq &p \leq \hat{p}+1.96\sqrt{p(1-p)/n}
\end{align}
$$
「基礎統計学Ⅰ 統計学入門」の$11.5.3$の議論により、「$n$が大きい場合大数の法則に基づいて、不等号の一番左と右の式における$p$は$\hat{p}$で近似できる」と考えることができる。よって、$95$%区間は下記のように近似することができる。
$$
\begin{align}
\hat{p}-1.96\sqrt{\hat{p}(1-\hat{p})/n} \leq p \leq \hat{p}+1.96\sqrt{\hat{p}(1-\hat{p})/n}
\end{align}
$$

母比率の区間推定(ポアソン分布)

解法の整理

区間推定

点推定

検定

まとめ

問題演習(基礎統計学Ⅰ)

https://www.hello-statisticians.com/explain-books-cat/toukeigakunyuumon-akahon/ch10_practice.html
https://www.hello-statisticians.com/explain-books-cat/toukeigakunyuumon-akahon/ch11_practice.html
https://www.hello-statisticians.com/explain-books-cat/toukeigakunyuumon-akahon/ch12_practice.html

クラメールの公式(Cramer’s rule)を用いた連立一次方程式の解法

クラメールの公式(Cramer’s rule)を用いることで連立一次方程式を行列式(determinant)の計算に基づいて解くことが可能になります。当記事ではクラメールの公式の式と問題演習を通した計算例の確認について取りまとめを行いました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$4$章「行列式」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

クラメールの公式

$n$個の未知数についての$n$個方程式からなる連立$1$次方程式$A \mathbf{x} = \mathbf{b}$を下記のように定義する。
$$
\large
\begin{align}
A \mathbf{x} &= \mathbf{b} \\
A &= \left( \begin{array}{ccc} \mathbf{a}_{1} & \cdots & \mathbf{a}_{n} \end{array} \right) \\
\mathbf{a}_{i} &= \left( \begin{array}{c} a_{1i} \\ \vdots \\ \mathbf{a}_{ni} \end{array} \right), \, \mathbf{x} = \left( \begin{array}{c} x_{1} \\ \vdots \\ x_{n} \end{array} \right), \, \mathbf{b} = \left( \begin{array}{c} b_{1} \\ \vdots \\ b_{n} \end{array} \right)
\end{align}
$$

このとき、連立$1$次方程式$A \mathbf{x} = \mathbf{b}$の解$\mathbf{x}$の$i$番目の成分$x_{i}$は下記に基づいて計算することができる。
$$
\large
\begin{align}
x_{i} = \frac{\det{\left( \begin{array}{ccccccc} \mathbf{a}_{1} & \cdots & \mathbf{a}_{i-1} & \mathbf{b} & \mathbf{a}_{i+1} & \cdots & \mathbf{a}_{n} \end{array} \right)}}{\det{(A)}}
\end{align}
$$

上記をクラメールの公式という。

クラメールの公式の使用例

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

基本例題$073$

・$[1]$
$$
\large
\begin{cases}
2x – 3y = 1 \\
3x + 4y = 2
\end{cases}
$$

上記の連立方程式は下記のような行列表記で表すことができる。
$$
\large
\begin{align}
\left( \begin{array}{cc} 2 & -3 \\ 3 & 4 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} 1 \\ 2 \end{array} \right)
\end{align}
$$

上記に対してクラメールの公式を適用することで、$x, y$は下記のように得られる。
$$
\large
\begin{align}
x &= \frac{\left| \begin{array}{cc} 1 & -3 \\ 2 & 4 \end{array} \right|}{\left| \begin{array}{cc} 2 & -3 \\ 3 & 4 \end{array} \right|} \\
&= \frac{4+6}{8+9} \\
&= \frac{10}{17} \\
y &= \frac{\left| \begin{array}{cc} 2 & 1 \\ 3 & 2 \end{array} \right|}{\left| \begin{array}{cc} 2 & -3 \\ 3 & 4 \end{array} \right|} \\
&= \frac{4-3}{8+9} \\
&= \frac{1}{17}
\end{align}
$$

対称行列における正定値行列・非負定値行列・負定値行列の性質とその導出

正定値行列(positive definite matrix)は行列に対する任意のベクトルの$2$次形式が常に正である場合の行列に対応します。当記事では正定値行列、非負定値行列、負定値行列について成立する式や性質とその導出について取り扱いました。
「統計学のための数学入門$30$講」の$23.2$節の「正定値行列・非負定値行列・負定値行列の性質」や「Matrix Computations」Section$4.2$の「Positive Definite Systems」を参考に作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

正定値行列・非負定値行列・負定値行列の定義と性質

正方行列$A$の$2$次形式

$p$次正方行列$A \in \mathbb{R}^{p \times p}$のベクトル$\mathbf{x} \in \mathbb{R}^{p}$についての$2$次形式を$q(\mathbf{x})$とおくと、$q(\mathbf{x})$は下記のように定義される。
$$
\large
\begin{align}
q(\mathbf{x}) = \mathbf{x}^{\mathrm{T}} A \mathbf{x}
\end{align}
$$

$2$次形式について詳しくは下記で取り扱った。

正定値行列・非負定値行列・負定値行列の定義

$p$次正方行列$A \in \mathbb{R}^{p \times p}$が正定値行列であることは任意のベクトル$\mathbf{x} \in \mathbb{R}^{p}$について下記が成立することと同値である。
$$
\large
\begin{align}
\mathbf{x}^{\mathrm{T}} A \mathbf{x} > 0
\end{align}
$$

正方行列$A$が正定値行列であることを$A > O$のように表す場合もある。同様に$A$が非負定値行列・負定値行列である場合は下記のように表される。

・非負定値行列
$$
\large
\begin{align}
A \geq O \iff \mathbf{x}^{\mathrm{T}} A \mathbf{x} \geq 0, \, {}^{\forall} \mathbf{x} \in \mathbb{R}^{p}
\end{align}
$$

・負定値行列
$$
\large
\begin{align}
A < O \iff \mathbf{x}^{\mathrm{T}} A \mathbf{x} < 0, \, {}^{\forall} \mathbf{x} \in \mathbb{R}^{p}
\end{align}
$$

対称行列における正定値行列・非負定値行列・負定値行列の性質

$A$が$p$次対称行列である場合、下記がそれぞれ成立する。
$(1) \,$ $A > O$ $\, \iff \,$ 行列$A$の固有値が全て正
$(2) \,$ $A \geq O$ $\, \iff \,$ 行列$A$の固有値が全て$0$以上
$(3) \,$ $A < O$ $\, \iff \,$ 行列$A$の固有値が全て負
$(4) \,$ $A > O \, \implies \, \det{(A)} > 0$
$(5) \,$ $A \geq O \, \implies \, \det{(A)} \geq 0$
$(6) \,$ $A > O$なら逆行列$A^{-1} > O$が存在し、$A^{-1} > O$である
$(7) \,$ $A \geq B > O$なら$|A| \geq |B|$である
$(8) \,$ $A \geq B > O$なら$B^{-1} \geq A^{-1}$である
$(9) \,$ $B \in \mathbb{R}^{p \times q}$のとき、$B^{\mathrm{T}}B$と$BB^{\mathrm{T}}$の正の固有値の個数と固有値の値は一致する
$(10) \,$ $A > O$のとき$-A < O$である

導出

$(1)$〜$(3)$の導出

$(1)$の導出

・「$A > O$ $\, \implies \,$ $\lambda > 0$」の導出
行列$A \in \mathbb{R}^{p \times p}$の固有値を$\lambda$、$\lambda$に対応する長さ$1$の固有ベクトルを$\mathbf{x} \in \mathbb{R}^{p}$とおくと下記が成立する。
$$
\large
\begin{align}
A \mathbf{x} = \lambda \mathbf{x}
\end{align}
$$

ここで上記の式に$\mathbf{x}^{\mathrm{T}}$を左からかけるとき、$\lambda$がスカラーであるので下記のような式変形が得られる。
$$
\large
\begin{align}
\mathbf{x}^{\mathrm{T}} A \mathbf{x} &= \mathbf{x}^{\mathrm{T}} (\lambda \mathbf{x}) \\
\mathbf{x}^{\mathrm{T}} A \mathbf{x} &= \lambda \mathbf{x}^{\mathrm{T}} \mathbf{x} \\
&= \lambda
\end{align}
$$

$A > O$であるとき${}^{\forall} \mathbf{x} \in \mathbb{R}^{p}$について$\mathbf{x}^{\mathrm{T}} A \mathbf{x} > 0$であるので「$A > O$ $\, \implies \,$ $\lambda > 0$」が成立する。

・「$\lambda > 0$ $\, \implies \,$ $A > O$」の導出
$A$の固有値を$\lambda_{i}, \, i=1, \cdots , p$とおき、$\lambda_{i} > 0$を仮定する。このとき$\lambda_{i}$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{i}$とおくと、スペクトル分解の式に基づいて下記が成立する。
$$
\large
\begin{align}
A = \sum_{i=1}^{p} \lambda_{i} \mathbf{u}_{i} \mathbf{u}_{i}^{\mathrm{T}}
\end{align}
$$

ここで$\mathbf{x}^{\mathrm{T}} A \mathbf{x}$に上記を代入することで下記が得られる。
$$
\large
\begin{align}
\mathbf{x}^{\mathrm{T}} A \mathbf{x} &= \mathbf{x}^{\mathrm{T}} \left( \sum_{i=1}^{p} \lambda_{i} \mathbf{u}_{i} \mathbf{u}_{i}^{\mathrm{T}} \right) \mathbf{x} \\
&= \sum_{i=1}^{p} \lambda_{i} \mathbf{x}^{\mathrm{T}} \mathbf{u}_{i} \mathbf{u}_{i}^{\mathrm{T}} \mathbf{x} \\
&= \sum_{i=1}^{p} \lambda_{i} (\mathbf{x}^{\mathrm{T}} \mathbf{u}_{i})^{2}
\end{align}
$$

「$\lambda_{i} > 0$を仮定した」かつ「$\mathbf{x}$が固有ベクトル$\mathbf{u}_{1}, \cdots , \mathbf{u}_{p}$全てと直交することはない1」ので、上記より$\mathbf{x}^{\mathrm{T}} A \mathbf{x} > 0$が成立する。よって$A > O$である。

$(2), (3)$の導出

$(1)$式の$>$を$\geq$や$<$に置き換えることで同様に示すことができる。

$(4), (5)$の導出

固有値分解の式より、直交行列$U$を用いて$A = U^{\mathrm{T}} \Lambda U$が成立する。このとき$|U|=1$であるので$\det{(A)} = |A|$は下記のように計算できる。
$$
\large
\begin{align}
\det{(A)} &= |A| = |U^{\mathrm{T}} \Lambda U| \\
&= |U^{\mathrm{T}}| |\Lambda| |U| \\
&= |\Lambda| = \prod_{i=1}^{p} \lambda_{i}
\end{align}
$$

ここで「$A > O \, \iff \, \lambda_{i} > 0$」、「$A \geq O \, \iff \, \lambda_{i} \geq 0$」なので、$(4)$と$(5)$がそれぞれ成立する。

$(6)$の導出

正定値行列は固有値が全て正であるので、行列式が正であり逆行列を持つ。ここで$A = U^{\mathrm{T}} \Lambda U$のように固有値分解の式を仮定するとき、$A^{-1}$について下記が成立する。
$$
\large
\begin{align}
A^{-1} &= (U^{\mathrm{T}} \Lambda U)^{-1} \\
&= U^{-1} \Lambda^{-1} (U^{\mathrm{T}})^{-1} \\
&= U^{\mathrm{T}} \Lambda^{-1} U
\end{align}
$$

ここで$\Lambda$の対角成分を$\lambda_{1}, \cdots , \lambda_{p}$とおくと、$\Lambda^{-1}$の対角成分は$\displaystyle \frac{1}{\lambda_{1}} , \cdots , \frac{1}{\lambda_{p}}$で表される。

$\lambda_{i}>0, \, {}^{\forall} i \in \{ 1, \cdots , p \}$のとき$\displaystyle \frac{1}{\lambda_{i}}>0, \, {}^{\forall} i \in \{ 1, \cdots , p \}$が成立するので$A > O$のとき$A^{-1} > O$である。

$(7), (8)$の導出

$(9)$の導出

$\lambda$を$B^{\mathrm{T}} B$の固有値、$\mathbf{x}$を対応する固有ベクトルとする。このとき下記が成立する。
$$
\large
\begin{align}
B^{\mathrm{T}} B \mathbf{x} = \lambda \mathbf{x}
\end{align}
$$

ここで上記に$B$を左からかけると下記が成立する。
$$
\large
\begin{align}
B B^{\mathrm{T}} B \mathbf{x} &= B \lambda \mathbf{x} \\
(B B^{\mathrm{T}}) B \mathbf{x} &= \lambda B \mathbf{x}
\end{align}
$$

上記より、$\lambda$が$B B^{\mathrm{T}}$の固有値、$B \mathbf{x}$が対応する固有ベクトルであることが確認できる。このことは全ての$\lambda$について成立するので「$B^{\mathrm{T}}B$と$BB^{\mathrm{T}}$の正の固有値の個数と固有値の値は一致する」ことが確認できる。

$(10)$の導出

$(1)$の導出と同様に示すことができる。

  1. $\mathbf{u}_{1}, \cdots \mathbf{u}_{p}$はそれぞれ$1$次独立であるから、$\mathbf{x}$が全ての固有ベクトル$\mathbf{u}_{1}, \cdots \mathbf{u}_{p}$とは直交できない。 ↩︎

行列の対角化③ ユニタリ行列を用いたエルミート行列の対角化(diagonalization)

行列の各固有値(eigen value)に対応する固有ベクトル(eigen vector)を用いることで行列の対角化(diagonalization)を行うことが可能です。当記事ではユニタリ行列を用いたエルミート行列の対角化(diagonalization)について取りまとめました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$8$章「固有値問題と行列の対角化」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

ユニタリ行列を用いたエルミート行列の対角化

標準内積と複素計量ベクトル空間

ユニタリ行列

行列$A$がユニタリ行列であるとき、ユニタリ行列の随伴行列$A^{*}$と単位行列$I$に関して下記が成立する。
$$
\large
\begin{align}
A A^{*} = A^{*} A = I
\end{align}
$$

上記より$A^{*}=A^{-1}$も同時に成立する。

エルミート行列

$n$次正方行列$A=(a_{ij})$を考え、$a_{ij}$の複素共役を$\overline{a_{ij}}$のように定める。このとき$a_{ij}=\overline{a_{ji}}$が成立する場合$A$はエルミート行列であるといい、$A^{*}=A$のように表す。

対角成分に関しては$a_{ii}=\overline{a_{ii}}$が成立するので、エルミート行列の対角成分が全て実数であることも抑えておくと良い。

テプリッツの定理

$A$が$n$次正方行列であるとき、「行列$A$がユニタリ行列で対角化できる」$\iff$「$A$が正規行列であり$A A^{*} = A^{*} A$」が成立する。

エルミート行列と正規行列

エルミート行列の定義より$A^{*}=A$が成立する。よって$A A^{*} = A^{2} = A^{*} A$が成立するのでエルミート行列は正規行列であり、テプリッツの定理が適用できる。

計算例

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

重要例題$083$

$$
\large
\begin{align}
A = \left( \begin{array}{ccc} 3 & i & -1 \\ -i & 5 & i \\ -1 & -i & 3 \end{array} \right)
\end{align}
$$

行列$A$の固有多項式を$F_{A}(t)=\det{(tI_{3} \, – \, A)}$とおくと、$F_{A}(t)$は下記のように計算できる。
$$
\large
\begin{align}
F_{A}(t) &= \det{(tI_{3} \, – \, A)} = \left| \begin{array}{ccc} t \, – \, 3 & -i & 1 \\ i & t \, – \, 5 & -i \\ 1 & i & t \, – \, 3 \end{array} \right| \\
&= \left| \begin{array}{ccc} 0 & -i(t \, – 2) & -(t \, – \, 2)(t \, – \, 4) \\ 0 & t \, – \, 4 & -i(t \, – \, 2) \\ 1 & i & t \, – \, 3 \end{array} \right| \\
&= 1 \cdot (-1)^{1+3} \left| \begin{array}{cc} -i(t \, – 2) & -(t \, – \, 2)(t \, – \, 4) \\ t \, – \, 4 & -i(t \, – \, 2) \end{array} \right| \\
&= i^{2}(t \, – \, 2)^{2} + (t \, – \, 2)(t \, – \, 4)^{2} \\
&= -(t \, – \, 2)^{2} + (t \, – \, 2)(t^{2} \, – \, 8t + 16) \\
&= (t \, – \, 2)(t^{2} \, – \, 8t + 16 \, – \, (t \, – \, 2)) \\
&= (t \, – \, 2)(t^{2} \, – \, 9t + 18) \\
&= (t \, – \, 2)(t \, – \, 3)(t \, – \, 6)
\end{align}
$$

上記より、行列$A$の固有値は$2, \, 3, \, 6$である。以下、それぞれの固有値について長さ$1$の固有ベクトルを求める。

・固有値が$2$のとき
行列$A \, – 2I_{3}$は下記のように行基本変形を行うことができる。
$$
\large
\begin{align}
A \, – 2I_{3} &= \left( \begin{array}{ccc} 1 & i & -1 \\ -i & 3 & i \\ -1 & -i & 1 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & -1 \\ -i & 3 & i \\ 0 & 0 & 0 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & -1 \\ 0 & 2 & 0 \\ 0 & 0 & 0 \end{array} \right)\longrightarrow \left( \begin{array}{ccc} 1 & i & -1 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{array} \right)
\end{align}
$$

上記より固有値$2$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{1}$とおくと、$\mathbf{u}_{1}$は下記のように表せる。
$$
\large
\begin{align}
\mathbf{u}_{1} = \frac{1}{\sqrt{2}} \left( \begin{array}{c} 1 \\ 0 \\ 1 \end{array} \right)
\end{align}
$$

・固有値が$3$のとき
行列$A \, – 3I_{3}$は下記のように行基本変形を行うことができる。
$$
\large
\begin{align}
A \, – 3I_{3} &= \left( \begin{array}{ccc} 0 & i & -1 \\ -i & 2 & i \\ -1 & -i & 0 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & 0 \\ -i & 2 & i \\ 0 & i & -1 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & 0 \\ 0 & 1 & i \\ 0 & i & -1 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & i \\ 0 & i & -1 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & i \\ 0 & 0 & 0 \end{array} \right)
\end{align}
$$

上記より固有値$3$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{3}$とおくと、$\mathbf{u}_{3}$は下記のように表せる。
$$
\large
\begin{align}
\mathbf{u}_{2} = \frac{1}{\sqrt{3}} \left( \begin{array}{c} -1 \\ -i \\ 1 \end{array} \right)
\end{align}
$$

・固有値が$6$のとき
行列$A \, – 6I_{3}$は下記のように行基本変形を行うことができる。
$$
\large
\begin{align}
A \, – 3I_{3} &= \left( \begin{array}{ccc} -3 & i & -1 \\ -i & -1 & i \\ -1 & -i & -3 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & 3 \\ i & 1 & -i \\ 3 & -i & -3 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & 3 \\ 0 & 2 & -4i \\ 0 & -4i & -8 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & i & 3 \\ 0 & 1 & -2i \\ 0 & i & 2 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & -2i \\ 0 & i & 2 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & -2i \\ 0 & 0 & 0 \end{array} \right)
\end{align}
$$

上記より固有値$6$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{3}$とおくと、$\mathbf{u}_{3}$は下記のように表せる。
$$
\large
\begin{align}
\mathbf{u}_{3} = \frac{1}{\sqrt{6}} \left( \begin{array}{c} -1 \\ 2i \\ 1 \end{array} \right)
\end{align}
$$

ここで$\displaystyle U=\left( \begin{array}{ccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} \end{array} \right)$とおくと、$U$は下記のように表すことができる。
$$
\large
\begin{align}
U = \left( \begin{array}{ccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} \end{array} \right) = \frac{1}{\sqrt{6}} \left( \begin{array}{ccc} \sqrt{3} & -\sqrt{2} & -1 \\ 0 & -\sqrt{2}i & 2i \\ \sqrt{3} & \sqrt{2} & 1 \end{array} \right)
\end{align}
$$

上記の$U$はユニタリ行列であり、下記のように$A$の対角化を行うことができる。
$$
\large
\begin{align}
U^{*} A U &= \frac{1}{6} \left( \begin{array}{ccc} \sqrt{3} & -\sqrt{2} & -1 \\ 0 & \sqrt{2}i & -2i \\ \sqrt{3} & \sqrt{2} & 1 \end{array} \right)^{\mathrm{T}} \left( \begin{array}{ccc} 3 & i & -1 \\ -i & 5 & i \\ -1 & -i & 3 \end{array} \right) \left( \begin{array}{ccc} \sqrt{3} & -\sqrt{2} & -1 \\ 0 & -\sqrt{2}i & 2i \\ \sqrt{3} & \sqrt{2} & 1 \end{array} \right) \\
&= \frac{1}{6} \left( \begin{array}{ccc} \sqrt{3} & 0 & \sqrt{3} \\ -\sqrt{2} & \sqrt{2}i & \sqrt{2} \\ -1 & -2i & 1 \end{array} \right) \left( \begin{array}{ccc} 3 & i & -1 \\ -i & 5 & i \\ -1 & -i & 3 \end{array} \right) \left( \begin{array}{ccc} \sqrt{3} & -\sqrt{2} & -1 \\ 0 & -\sqrt{2}i & 2i \\ \sqrt{3} & \sqrt{2} & 1 \end{array} \right) \\
&= \frac{1}{6} \left( \begin{array}{ccc} 2\sqrt{3} & 0 & 2\sqrt{3} \\ -3\sqrt{2} & 3\sqrt{2}i & 3\sqrt{2} \\ -6 & -12i & 6 \end{array} \right) \left( \begin{array}{ccc} \sqrt{3} & -\sqrt{2} & -1 \\ 0 & -\sqrt{2}i & 2i \\ \sqrt{3} & \sqrt{2} & 1 \end{array} \right) \\
&= \frac{1}{6} \left( \begin{array}{ccc} 12 & 0 & 0 \\ 0 & 18 & 0 \\ 0 & 0 & 36 \end{array} \right) \\
&= \left( \begin{array}{ccc} 2 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 6 \end{array} \right)
\end{align}
$$

基本変形・基本行列(elementary matrix)と行列式(determinant)

行基本変形は基本行列(elementary matrix)の積による操作によって表すことができるなど、基本行列は
よく出てくるので抑えておくと良いです。当記事では基本変形・基本行列(elementary matrix)の行列式(determinant)の計算について取りまとめを行いました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$3$節「行列の構造」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

基本行列の積と行列式

行列の積の行列式

行列$A$と$B$の積$AB$の行列式$\det{(AB)}$について下記が成立する。
$$
\large
\begin{align}
\det{(AB)} = \det{(A)} \det{(B)} \quad (1.1)
\end{align}
$$

基本行列の行列式

基本行列$P_{ij}, P_{i}(c), P_{ij}(a)$の行列式はそれぞれ下記のように得られる。
$$
\large
\begin{align}
\det{(P_{ij})} &= -1 \\
\det{(P_{i}(c))} &= c \\
\det{(P_{ij}(a))} &= 1
\end{align}
$$

基本行列の定義については下記で詳しく取り扱った。

基本行列の積と行列式

$(1.1)$式と「基本行列の行列式」を元に計算することができる。

行列式の多重線形性・交代性

行列式の還元定理

具体例の確認

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

基本例題$065$

・$[1]$
行列$A$の$i$列目と$j$列目を入れ替えた行列は$AP_{ij}$のように表される。$AP_{ij}$の行列式$\det{(AP_{ij})}$は下記のように計算できる。
$$
\large
\begin{align}
\det{(AP_{ij})} &= \det{(A)} \det{(P_{ij})} \\
&= -\det{(A)}
\end{align}
$$

・$[2]$
行列$A$の$i$列目を$c$倍した行列は$AP_{i}(c)$のように表される。$AP_{i}(c)$の行列式$\det{(AP_{i}(c))}$は下記のように計算できる。
$$
\large
\begin{align}
\det{(AP_{i}(c))} &= \det{(A)} \det{(P_{i}(c))} \\
&= c \det{(A)}
\end{align}
$$

・$[3]$
行列$A$の$i$列目と$j$列目を入れ替えた行列は$AP_{ij}(a)$のように表される。$AP_{ij}(a)$の行列式$\det{(AP_{ij})}$は下記のように計算できる。
$$
\large
\begin{align}
\det{(AP_{ij}(a))} &= \det{(A)} \det{(P_{ij}(a))} \\
&= \det{(A)}
\end{align}
$$

基本例題$066$

(1)

・$[1]$
$$
\large
\begin{align}
A = \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) \quad (2.1)
\end{align}
$$

上記のように定義される$n$次正方行列について行列式の行交代性から下記が成立する。
$$
\large
\begin{align}
\det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } = -\det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } = -\det{(A)}
\end{align}
$$

・$[2]$
$(2.1)$式で定義される$n$次正方行列$A$について行列式の行多重線形性から下記が成立する。
$$
\large
\begin{align}
\det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ c \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } = c\det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } = c \det{(A)}
\end{align}
$$

・$[3]$
$(2.1)$式で定義される$n$次正方行列$A$について行列式の行多重線形性と行交代性から下記が成立する。
$$
\large
\begin{align}
\det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} + c \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } &= \det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } + c \det{ \left(\begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } \\
&= \det{(A)}
\end{align}
$$

(2)

列多重線形性と列交代性を用いて基本例題$066 \, (1)$と同様に示すことができる。

基本例題$067$

・$[1]$
$$
\large
\begin{align}
\left(\begin{array}{ccc} 1 & 2 & 4 \\ 3 & 9 & 27 \\ -1 & 1 & -1 \end{array} \right)
\end{align}
$$

下記のように行列式を計算することができる。
$$
\large
\begin{align}
\left| \begin{array}{ccc} 1 & 2 & 4 \\ 3 & 9 & 27 \\ -1 & 1 & -1 \end{array} \right| &= -\left| \begin{array}{ccc} -1 & 1 & -1 \\ 3 & 9 & 27 \\ 1 & 2 & 4 \end{array} \right| \\
&= \left| \begin{array}{ccc} 1 & -1 & 1 \\ 3 & 9 & 27 \\ 1 & 2 & 4 \end{array} \right| \\
&= \left| \begin{array}{ccc} 1 & -1 & 1 \\ 0 & 12 & 24 \\ 0 & 3 & 3 \end{array} \right| \\
&= 12 \left| \begin{array}{ccc} 1 & -1 & 1 \\ 0 & 1 & 2 \\ 0 & 3 & 3 \end{array} \right| \\
&= 12 \left| \begin{array}{cc} 1 & 2 \\ 3 & 3 \end{array} \right| \\
&= 12 (3-6) = -36
\end{align}
$$

・$[2]$
$$
\large
\begin{align}
\left( \begin{array}{cccc} 1 & 0 & 1 & 0 \\ 2 & 1 & 2 & 1 \\ 3 & 2 & 4 & 2 \\ 4 & 3 & 6 & 4 \end{array} \right)
\end{align}
$$

下記のように行列式を計算することができる。
$$
\large
\begin{align}
\left| \begin{array}{cccc} 1 & 0 & 1 & 0 \\ 2 & 1 & 2 & 1 \\ 3 & 2 & 4 & 2 \\ 4 & 3 & 6 & 4 \end{array} \right| &= \left| \begin{array}{cccc} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 2 & 1 & 2 \\ 0 & 3 & 2 & 4 \end{array} \right| \\
&= \left| \begin{array}{ccc} 1 & 0 & 1 \\ 2 & 1 & 2 \\ 3 & 2 & 4 \end{array} \right| \\
&= \left| \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 0 & 2 & 1 \end{array} \right| \\
&= \left| \begin{array}{cc} 1 & 0 \\ 2 & 1 \end{array} \right| \\
&= 1
\end{align}
$$

・$[3]$
$$
\large
\begin{align}
\left( \begin{array}{cccc} 2 & 3 & 0 & 1 \\ 0 & 1 & 5 & 0 \\ 1 & 0 & 0 & 4 \\ 0 & 7 & 2 & 0 \end{array} \right)
\end{align}
$$

下記のように行列式を計算することができる。
$$
\large
\begin{align}
\left| \begin{array}{cccc} 2 & 3 & 0 & 1 \\ 0 & 1 & 5 & 0 \\ 1 & 0 & 0 & 4 \\ 0 & 7 & 2 & 0 \end{array} \right| &= -\left| \begin{array}{cccc} 1 & 0 & 0 & 4 \\ 0 & 1 & 5 & 0 \\ 2 & 3 & 0 & 1 \\ 0 & 7 & 2 & 0 \end{array} \right| \\
&= -\left| \begin{array}{cccc} 1 & 0 & 0 & 4 \\ 0 & 1 & 5 & 0 \\ 0 & 3 & 0 & -7 \\ 0 & 7 & 2 & 0 \end{array} \right| \\
&= -\left| \begin{array}{ccc} 1 & 5 & 0 \\ 3 & 0 & -7 \\ 7 & 2 & 0 \end{array} \right| \\
&= -\left| \begin{array}{ccc} 1 & 5 & 0 \\ 0 & -15 & -7 \\ 0 & -33 & 0 \end{array} \right| \\
&= \left| \begin{array}{ccc} 1 & 5 & 0 \\ 0 & 15 & 7 \\ 0 & -33 & 0 \end{array} \right| \\
&= \left| \begin{array}{cc} 15 & 7 \\ -33 & 0 \end{array} \right| \\
&= -231
\end{align}
$$

行列の対角化② 交代行列の対角化(diagonalization)

行列の各固有値(eigen value)に対応する固有ベクトル(eigen vector)を用いることで行列の対角化(diagonalization)を行うことが可能です。当記事では交代行列の対角化(diagonalization)について取りまとめました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$8$章「固有値問題と行列の対角化」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

交代行列の対角化

交代行列

行列$A=(a_{ij})$の随伴行列(Adjoint matrix)を$A^{*}=(\overline{a_{ji}})$のように定義する。

このとき、正規行列・歪エルミート行列・交代行列は下記のような式で定義される。

名称 成立する式
正規行列$AA^{*}=A^{*}A$
歪エルミート行列 $A^{*}=-A$
交代行列 $A^{\mathrm{T}}=-A$

$A^{*}=-A$のとき$AA^{*}=-A^{2}=A^{*}A$、$A^{\mathrm{T}}=-A$のとき$AA^{*}=AA^{\mathrm{T}}=-A^{2}=A^{\mathrm{T}}A=A^{*}A$が成立するので歪エルミート行列と交代行列はどちらも正規行列である。また、交代行列は要素が実数の歪エルミート行列である1

テプリッツの定理

$A$が$n$次正方行列であるとき、「行列$A$がユニタリ行列で対角化できる」$\iff$「$A$が正規行列である」が成立する。

標準内積と複素計量ベクトル空間

$\mathbb{C}$上の$n$次元ベクトル空間$\mathbb{C}^{n}$の$2$つのベクトル$\mathbf{u}, \mathbf{v} \in \mathbb{C}^{n}$を下記のように定義する。
$$
\large
\begin{align}
\mathbf{u} = \left( \begin{array}{c} u_{1} \\ \vdots \\ u_{n} \end{array} \right), \quad \mathbf{v} = \left( \begin{array}{c} v_{1} \\ \vdots \\ v_{n} \end{array} \right)
\end{align}
$$

このとき、$\mathbb{C}^{n}$上の標準内積$(\mathbf{u}, \mathbf{v})$が下記のように定義される。
$$
\large
\begin{align}
(\mathbf{u}, \mathbf{v}) = \overline{\mathbf{v}}^{\mathrm{T}} \mathbf{u} = \sum_{i=1}^{n} u_{i} \overline{v_{i}}
\end{align}
$$

上記のように$\mathbb{C}^{n}$上の標準内積$(\mathbf{u}, \mathbf{v})$を元に定義される計量ベクトル空間を複素計量ベクトル空間という。標準内積$(\mathbf{u}, \mathbf{v})$はベクトルの正規化を行う際などに用いられる。

計算例

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

基本例題$158$

$$
\large
\begin{align}
A = \left( \begin{array}{ccc} 0 & -1 & -2 \\ 1 & 0 & -2 \\ 2 & 2 & 0 \end{array} \right)
\end{align}
$$

行列$A$の固有多項式を$F_{A}(t)=\det{(tI_{3} \, – \, A)}$とおくと、$F_{A}(t)$は下記のように計算できる。
$$
\large
\begin{align}
F_{A}(t) &= \det{(tI_{3} \, – \, A)} = \left| \begin{array}{ccc} t & 1 & 2 \\ -1 & t & 2 \\ -2 & -2 & t \end{array} \right| \\
&= \left| \begin{array}{ccc} 0 & t^{2}+1 & 2(t+1) \\ -1 & t & 2 \\ 0 & -2(t+1) & t-4 \end{array} \right| \\
&= (-1) \cdot (-1)^{2+1} \left| \begin{array}{cc} t^{2}+1 & 2(t+1) \\ -2(t+1) & t-4 \end{array} \right| \\
&= (t^{2}+1)(t \, – \, 4) + 4(t+1)^{2} \\
&= t^{3} \, – \, \cancel{4t^{2}} + t \, – \, \cancel{4} + \cancel{4t^{2}} + 8t + \cancel{4} \\
&= t^{3} + 9t \\
&= t(t^{2} + 9)
\end{align}
$$

上記より、行列$A$の固有値は$0, \, 3i \, -3i$である。以下、それぞれの固有値について長さ$1$の固有ベクトルを求める。

・固有値が$0$のとき
行列$A$は下記のように行基本変形を行うことができる。
$$
\large
\begin{align}
\left( \begin{array}{ccc} 0 & -1 & -2 \\ 1 & 0 & -2 \\ 2 & 2 & 0 \end{array} \right) & \longrightarrow \left( \begin{array}{ccc} 1 & 0 & -2 \\ 0 & -1 & -2 \\ 2 & 2 & 0 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & -2 \\ 0 & 1 & 2 \\ 0 & 2 & 4 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & -2 \\ 0 & 1 & 2 \\ 0 & 0 & 0 \end{array} \right)
\end{align}
$$

上記より固有値$0$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{1}$とおくと、$\mathbf{u}_{1}$は下記のように表せる。
$$
\large
\begin{align}
\mathbf{u}_{1} = \frac{1}{3} \left( \begin{array}{c} 2 \\ -2 \\ 1 \end{array} \right)
\end{align}
$$

・固有値が$3i$のとき
行列$A \, – \, 3i I_{3}$は下記のように行基本変形を行うことができる。
$$
\large
\begin{align}
\left( \begin{array}{ccc} -3i & -1 & -2 \\ 1 & -3i & -2 \\ 2 & 2 & -3i \end{array} \right) & \longrightarrow \left( \begin{array}{ccc} 1 & -3i & -2 \\ 2 & 2 & -3i \\ -3i & -1 & -2 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & -3i & -2 \\ 0 & 2+6i & 4 \, – \, 3i \\ 0 & 8 & -2 \, – \, 6i \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & -3i & -2 \\ 0 & 8 & -2 \, – \, 6i \\ 0 & 2+6i & 4 \, – \, 3i \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & -3i & -2 \\ 0 & 1 & \displaystyle -\frac{1+3i}{4} \\ 0 & 2+6i & 4 \, – \, 3i \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & -3i & -2 \\ 0 & 1 & \displaystyle -\frac{1+3i}{4} \\ 0 & 0 & 0 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & \displaystyle \frac{1 \, – \, 3i}{4} \\ 0 & 1 & \displaystyle -\frac{1+3i}{4} \\ 0 & 0 & 0 \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 0 & \displaystyle -\frac{3i \, – \, 1}{4} \\ 0 & 1 & \displaystyle -\frac{3i+1}{4} \\ 0 & 0 & 0 \end{array} \right)
\end{align}
$$

上記より固有値$3i$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{2}$とおくと、$\mathbf{u}_{2}$は下記のように表せる。
$$
\large
\begin{align}
\mathbf{u}_{2} = \frac{1}{6} \left( \begin{array}{c} 3i \, – \, 1 \\ 3i+1 \\ 4 \end{array} \right)
\end{align}
$$

・固有値が$-3i$のとき
行列$A + 3i I_{3}$は下記のように行基本変形を行うことができる。
$$
\large
\begin{align}
\left( \begin{array}{ccc} 3i & -1 & -2 \\ 1 & 3i & -2 \\ 2 & 2 & 3i \end{array} \right) & \longrightarrow \left( \begin{array}{ccc} 1 & 3i & -2 \\ 3i & -1 & -2 \\ 2 & 2 & 3i \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 3i & -2 \\ 0 & 8 & -2+6i \\ 0 & 2 \, – \, 6i & 4+3i \end{array} \right) \\
& \longrightarrow \left( \begin{array}{ccc} 1 & 3i & -2 \\ 0 & 1 & \displaystyle \frac{3i \, – \, 1}{4} \\ 0 & 0 & 0 \end{array} \right) \longrightarrow \left( \begin{array}{ccc} 1 & 0 & \displaystyle \frac{3i+1}{4} \\ 0 & 1 & \displaystyle \frac{3i \, – \, 1}{4} \\ 0 & 0 & 0 \end{array} \right)
\end{align}
$$

上記より固有値$-3i$に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{3}$とおくと、$\mathbf{u}_{3}$は下記のように表せる。
$$
\large
\begin{align}
\mathbf{u}_{3} = \frac{1}{6} \left( \begin{array}{c} -3i \, – \, 1 \\ -3i+1 \\ 4 \end{array} \right)
\end{align}
$$

ここで$\displaystyle U=\left( \begin{array}{ccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} \end{array} \right)$とおくと、$U$は下記のように表すことができる。
$$
\large
\begin{align}
U = \left( \begin{array}{ccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} \end{array} \right) = \frac{1}{6} \left( \begin{array}{ccc} 4 & 3i \, – \, 1 & -3i \, – \, 1 \\ -4 & 3i+1 & -3i+1 \\ 2 & 4 & 4 \end{array} \right)
\end{align}
$$

上記の$U$はユニタリ行列であり、下記のように$A$の対角化を行うことができる。
$$
\large
\begin{align}
U^{*} A U &= \frac{1}{6^{2}} \left( \begin{array}{ccc} 4 & -3i \, – \, 1 & 3i \, – \, 1 \\ -4 & -3i+1 & 3i+1 \\ 2 & 4 & 4 \end{array} \right)^{\mathrm{T}} \left( \begin{array}{ccc} 0 & -1 & -2 \\ 1 & 0 & -2 \\ 2 & 2 & 0 \end{array} \right) \left( \begin{array}{ccc} 4 & 3i \, – \, 1 & -3i \, – \, 1 \\ -4 & 3i+1 & -3i+1 \\ 2 & 4 & 4 \end{array} \right) \\
&= \frac{1}{6^{2}} \left( \begin{array}{ccc} 4 & -4 & 2 \\ -3i \, – \, 1 & -3i+1 & 4 \\ 3i \, – \, 1 & 3i+1 & 4 \end{array} \right) \left( \begin{array}{ccc} 0 & -1 & -2 \\ 1 & 0 & -2 \\ 2 & 2 & 0 \end{array} \right) \left( \begin{array}{ccc} 4 & 3i \, – \, 1 & -3i \, – \, 1 \\ -4 & 3i+1 & -3i+1 \\ 2 & 4 & 4 \end{array} \right) \\
&= \frac{1}{6^{2}} \left( \begin{array}{ccc} 0 & 0 & 0 \\ -3i+9 & 3i+9 & 12i \\ 3i+9 & -3i+9 & -12i \end{array} \right) \left( \begin{array}{ccc} 4 & 3i \, – \, 1 & -3i \, – \, 1 \\ -4 & 3i+1 & -3i+1 \\ 2 & 4 & 4 \end{array} \right) \\
&= \frac{1}{2^{2} \times 3} \left( \begin{array}{ccc} 0 & 0 & 0 \\ -i+3 & i+3 & 4i \\ i+3 & -i+3 & -4i \end{array} \right) \left( \begin{array}{ccc} 4 & 3i \, – \, 1 & -3i \, – \, 1 \\ -4 & 3i+1 & -3i+1 \\ 2 & 4 & 4 \end{array} \right) \\
&= \frac{1}{2^{2} \times 3} \left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 20i+16i & 0 \\ 0 & 0 & -20i \, – \, 16i \end{array} \right) \\
&= \left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 3i & 0 \\ 0 & 0 & -3i \end{array} \right)
\end{align}
$$

  1. 対称行列も$A=A^{\mathrm{T}}$より、$AA^{*}=AA^{\mathrm{T}}=A^{2}=A^{\mathrm{T}}A=A^{*}A$であるので正規行列であることも合わせて抑えておくと良い。 ↩︎

行列式①:行多重線形性・行交代性の導出とそれぞれの計算例の確認

線形代数を学ぶにあたって行列式(determinants)は固有値・固有ベクトルの導出にあたって出てくる固有多項式・固有方程式に出てくるなど、重要な概念です。当記事では行列式の計算によく用いられる行多重線形性と行交代性の式・導出とそれぞれの使用例について取りまとめを行いました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$4$章「行列式」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

行多重線形性と交代性の導出

置換

行多重線形性の式

$n$次正方行列$A$をそれぞれの行に対応するベクトル$\mathbf{a}_{1}, \mathbf{a}_{2}, \cdots , \mathbf{a}_{n}$を用いて下記のように定義する。
$$
\large
\begin{align}
A = \left( \begin{array}{c} \mathbf{a}_{1} \\ \mathbf{a}_{2} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right)
\end{align}
$$

ここで$A$の$i$行目に対応する$\mathbf{a}_{i}$がベクトル$\mathbf{b}, \mathbf{c}$と実数$k, l \in \mathbb{R}$を用いて下記のように表せると仮定する。
$$
\large
\begin{align}
\mathbf{a}_{i} = k \mathbf{b} + l \mathbf{c}
\end{align}
$$

このとき、行列$A$の行列式$\det{(A)}$について下記が成立する。
$$
\large
\begin{align}
\det{(A)} &= \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right)} \\
&= \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ k \mathbf{b} + l \mathbf{c} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right)} \\
&= k \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{b} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right)} + l \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{c} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right)} \quad (1.1)
\end{align}
$$

上記の$(1.1)$式を行列式の行多重線形性という。

行多重線形性の導出

行列$A=(a_{ij})$の行列式$\det{(A)}$は下記のように表される。
$$
\large
\begin{align}
\det{(A)} = \sum_{\sigma} \mathrm{sgn}(\sigma) a_{1 \sigma(1)} \cdots a_{i \sigma(i)} \cdots a_{n \sigma(n)} \quad (1.2)
\end{align}
$$

ここで行列$A$の$i$行の$\mathbf{a}_{i}$が$\mathbf{a}_{i}=k \mathbf{b} + l \mathbf{c}$のように表されるとき、$(1.2)$式は下記のように変形できる。
$$
\large
\begin{align}
\det{(A)} &= \sum_{\sigma} \mathrm{sgn}(\sigma) a_{1 \sigma(1)} \cdots a_{i \sigma(i)} \cdots a_{n \sigma(n)} \quad (1.2) \\
&= \sum_{\sigma} \mathrm{sgn}(\sigma) a_{1 \sigma(1)} \cdots (k b_{\sigma(i)} + l c_{\sigma(i)}) \cdots a_{n \sigma(n)} \\
&= k \sum_{\sigma} \mathrm{sgn}(\sigma) a_{1 \sigma(1)} \cdots b_{\sigma(i)} \cdots a_{n \sigma(n)} + l \sum_{\sigma} \mathrm{sgn}(\sigma) a_{1 \sigma(1)} \cdots c_{\sigma(i)} \cdots a_{n \sigma(n)}
\end{align}
$$

上記より、行多重線形性が成立することが示される。

行多重線形性は$A$の$i$行目を$\mathbf{b}$で置き換えた行列を$B$、$\mathbf{c}$で置き換えた行列を$C$とおいて$\det{(A)} = k\det{(B)} + l\det{(C)}$のように表すこともできる。

交代性の式

$n$次正方行列$A$をそれぞれの行に対応するベクトル$\mathbf{a}_{1}, \mathbf{a}_{2}, \cdots , \mathbf{a}_{n}$を用いて下記のように定義する。
$$
\large
\begin{align}
A = \left( \begin{array}{c} \mathbf{a}_{1} \\ \mathbf{a}_{2} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right)
\end{align}
$$

ここで$i \neq j$である行列$A$の$i$行目と$j$行目を入れ替えた行列の行列式について下記が成立する。
$$
\large
\begin{align}
\det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } = -\det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \vdots \\ \mathbf{a}_{i} \\ \vdots \\ \mathbf{a}_{j} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } \quad (1.3)
\end{align}
$$

上記の$(1.3)$式を行列式の行交代性という。

交代性の導出

置換の使用例

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

基本例題$060$

$$
\large
\begin{align}
A = \left( \begin{array}{c} \mathbf{a}_{1} \\ \mathbf{a}_{2} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right), \quad cA = \left( \begin{array}{c} c\mathbf{a}_{1} \\ c\mathbf{a}_{2} \\ \vdots \\ c\mathbf{a}_{n} \end{array} \right)
\end{align}
$$

$cA$は上記のように表されるので、行多重線形性を用いて$\det{(cA)}$は下記のように変形できる。
$$
\large
\begin{align}
\det{(cA)} &= \det{ \left( \begin{array}{c} c\mathbf{a}_{1} \\ c\mathbf{a}_{2} \\ \vdots \\ c\mathbf{a}_{n} \end{array} \right) } \\
&= c \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ c\mathbf{a}_{2} \\ \vdots \\ c\mathbf{a}_{n} \end{array} \right) } \\
&= c^{2} \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \mathbf{a}_{2} \\ \vdots \\ c\mathbf{a}_{n} \end{array} \right) } \\
&= \cdots \cdots \\
&= c^{n} \det{ \left( \begin{array}{c} \mathbf{a}_{1} \\ \mathbf{a}_{2} \\ \vdots \\ \mathbf{a}_{n} \end{array} \right) } \\
&= c^{n} \det{(A)}
\end{align}
$$

上記より$\det{(cA)} = c^{n} \det{(A)}$が成立する。

行列の対角化① 実対称行列の直交行列による対角化(diagonalization)

行列の各固有値(eigen value)に対応する固有ベクトル(eigen vector)を用いることで行列の対角化(diagonalization)を行うことが可能です。当記事では直交行列を用いた実対称行列の対角化について取りまとめました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$8$章「固有値問題と行列の対角化」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

直交行列を用いた実対象行列の対角化

正方行列の対角化可能条件

$n$次正方行列$A$の固有多項式を$F_{A}(t) = \det{(tI_{n} \, – \, A)}$とおくとき、$F_{A}(t)$が下記のように表せることを仮定する。
$$
\large
\begin{align}
F_{A}(t) &= \det{(tI_{n} \, – \, A)} = \prod_{i=1}^{r} (t \, – \, \lambda_{i})^{m_{i}} \\
i \neq j & \iff \lambda_{i} \neq \lambda_{j}
\end{align}
$$

上記のように$F_{A}(t)$が表せる場合、行列$A$の相異なる固有値は$\lambda_{1}, \cdots , \lambda_{r}$であり、$m_{i}$は対応する固有値$\lambda_{i}$の重複度を表す。このとき下記の$[1]$〜$[3]$が同値となる。
$[1] \,$ 行列$A$が対角化可能であり、$P^{-1} A P$が対角行列となる正則行列$P$が存在する。

$[2] \,$ $i=1, \cdots , r$について固有値$\lambda_{i}$に対応する固有空間を$W_{\lambda_{i}}$とするとき、下記が成立する。
$$
\large
\begin{align}
\mathrm{dim}{W_{\lambda_{1}}} + \cdots + \mathrm{dim}{W_{\lambda_{r}}} = n
\end{align}
$$

$[3] \,$ 各$i=1, \cdots , r$について$\mathrm{dim}{W_{\lambda_{i}}}=m_{i}$が成立する。

対称行列における固有空間の直交性

$A$が$n$次正方行列であることを仮定し、行列$A$の相異なる固有値$\lambda_{i}, \, \lambda_{j}$に対する固有空間を$W_{\lambda_{i}}, \, W_{\lambda_{j}}$のように定義する。

このとき固有空間の$W_{\lambda_{i}}$と$W_{\lambda_{j}}$は$\mathrm{R}^{n}$における標準内積について直交する。すなわち、任意の$\mathbf{v}_{i} \in W_{\lambda_{i}}, \, \mathbf{v}_{j} \in W_{\lambda_{j}}$について下記が成立する。
$$
\large
\begin{align}
(\mathbf{v}_{i}, \mathbf{v}_{j}) = \mathbf{v}_{i} \cdot \mathbf{v}_{j} = 0
\end{align}
$$

対角化の計算例

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

基本例題$157$

$$
\large
\begin{align}
A &= \left( \begin{array}{ccc} 4 & -1 & 1 \\ -1 & 4 & -1 \\ 1 & -1 & 4 \end{array} \right)
\end{align}
$$

行列$A$の固有多項式$\det{(\lambda I_{3} \, – \, A)}$は下記のように計算できる。
$$
\large
\begin{align}
\det{(\lambda I_{3} \, – \, A)} &= \left| \begin{array}{ccc} \lambda \, – \, 4 & 1 & -1 \\ 1 & \lambda \, – \, 4 & 1 \\ -1 & 1 & \lambda \, – \, 4 \end{array} \right| \\
&= \left| \begin{array}{ccc} 0 & 1-(\lambda \, – \, 4)^{2} & -1 \, – \, (\lambda \, – \, 4) \\ 1 & \lambda \, – \, 4 & 1 \\ 0 & 1+(\lambda \, – \, 4) & \lambda \, – \, 4 + 1 \end{array} \right| \\
&= \left| \begin{array}{ccc} 0 & -(\lambda \, – \, 5)(\lambda \, – \, 3) & -(\lambda \, – \, 3) \\ 1 & \lambda \, – \, 4 & 1 \\ 0 & \lambda \, – \, 3 & \lambda \, – \, 3 \end{array} \right| \\
&= (-1)^{2+1} \left| \begin{array}{cc} -(\lambda \, – \, 5)(\lambda \, – \, 3) & -(\lambda \, – \, 3) \\ \lambda \, – \, 3 & \lambda \, – \, 3 \end{array} \right| \\
&= \left| \begin{array}{cc} (\lambda \, – \, 5)(\lambda \, – \, 3) & \lambda \, – \, 3 \\ \lambda \, – \, 3 & \lambda \, – \, 3 \end{array} \right| \\
&= (\lambda \, – \, 3)^{2} (\lambda \, – \, 5) \, – \, (\lambda \, – \, 3)^{2} \\
&= (\lambda \, – \, 3)^{2} (\lambda \, – \, 6)
\end{align}
$$

よって行列$A$の固有値は$\lambda=3, 6$である。以下、それぞれの固有値に対してそれぞれ直交する長さ$1$の固有ベクトルを求める。

・$\lambda = 3$
長さ$1$の固有ベクトルを$\displaystyle \mathbf{u} = \left( \begin{array}{c} x \\ y \\ z \end{array} \right)$とおくとき、$A \mathbf{u} = 3 \mathbf{u}$より下記が得られる。
$$
\large
\begin{align}
A \mathbf{u} &= 3 \mathbf{u} \\
\left( \begin{array}{ccc} 4 & -1 & 1 \\ -1 & 4 & -1 \\ 1 & -1 & 4 \end{array} \right) \left( \begin{array}{c} x \\ y \\ z \end{array} \right) &= 3 \left( \begin{array}{c} x \\ y \\ z \end{array} \right) \\
\left( \begin{array}{c} 4x \, – \, y + z \\ -x + 4y \, – \, z \\ x \, – \, y + 4y \end{array} \right) &= 3 \left( \begin{array}{c} x \\ y \\ z \end{array} \right) \\
\left( \begin{array}{c} x \, – \, y + z \\ -x + y \, – \, z \\ x \, – \, y + y \end{array} \right) &= \left( \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right) \\
x \, – \, y + z &= 0
\end{align}
$$

$x \, – \, y + z = 0$かつ$|\mathbf{u}|=1$であるので、$\mathbf{u}$の$1$つ$\mathbf{u}_{1}$は下記のように表すことができる。
$$
\large
\begin{align}
\mathbf{u}_{1} = \frac{1}{\sqrt{2}} \left( \begin{array}{c} 1 \\ 1 \\ 0 \end{array} \right) \quad (1)
\end{align}
$$

ここで$x \, – \, y + z = 0$より$y = x+z$であるので、$(1)$と直交する$\lambda = 3$に対応する固有ベクトルを$\displaystyle \mathbf{u}_{2} = \left( \begin{array}{c} x \\ x+z \\ z \end{array} \right)$とおくと$\mathbf{u}_{2}$について下記が成立する。
$$
\large
\begin{align}
\mathbf{u}_{1} \cdot \mathbf{u}_{2} &= 0 \\
\frac{1}{\sqrt{2}} \left( \begin{array}{c} 1 \\ 1 \\ 0 \end{array} \right) \cdot \left( \begin{array}{c} x \\ x+z \\ z \end{array} \right) &= 0 \\
x + (x+z) &= 0 \\
2x + z &= 0 \\
z &= -2x
\end{align}
$$

上記と$x \, – \, y + z = 0, \, |\mathbf{u}_{2}|=1$より、$\mathbf{u}_{2}$は下記のように得られる。
$$
\large
\begin{align}
\mathbf{u}_{2} = \frac{1}{\sqrt{6}} \left( \begin{array}{c} -1 \\ 1 \\ 2 \end{array} \right) \quad (2)
\end{align}
$$

・$\lambda = 6$
長さ$1$の固有ベクトルを$\displaystyle \mathbf{u} = \left( \begin{array}{c} x \\ y \\ z \end{array} \right)$とおくとき、$A \mathbf{u} = 6 \mathbf{u}$より下記が得られる。
$$
\large
\begin{align}
A \mathbf{u} &= 6 \mathbf{u} \\
\left( \begin{array}{ccc} 4 & -1 & 1 \\ -1 & 4 & -1 \\ 1 & -1 & 4 \end{array} \right) \left( \begin{array}{c} x \\ y \\ z \end{array} \right) &= 6 \left( \begin{array}{c} x \\ y \\ z \end{array} \right) \\
\left( \begin{array}{c} -2x \, – \, y + z \\ -(x+2y+z) \\ x \, – \, y \, – \, 2z \end{array} \right) &= \left( \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right)
\end{align}
$$

$|\mathbf{u}|=1$であるので、$\mathbf{u}_{3}$は下記のように表すことができる。
$$
\large
\begin{align}
\mathbf{u}_{3} = \frac{1}{\sqrt{3}} \left( \begin{array}{c} 1 \\ -1 \\ 1 \end{array} \right) \quad (3)
\end{align}
$$

$(1)$〜$(3)$に基づいて直交行列$U$を下記のように定義する。
$$
\large
\begin{align}
U = \left( \begin{array}{ccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} \end{array} \right) = \frac{1}{\sqrt{6}} \left( \begin{array}{ccc} \sqrt{3} & -1 & \sqrt{2} \\ \sqrt{3} & 1 & -\sqrt{2} \\ 0 & 2 & \sqrt{2} \end{array} \right)
\end{align}
$$

このとき、$U^{\mathrm{T}} A U$は下記のように計算できる。
$$
\large
\begin{align}
U^{\mathrm{T}} A U &= \frac{1}{6} \left( \begin{array}{ccc} \sqrt{3} & \sqrt{3} & 0 \\ -1 & 1 & 2 \\ \sqrt{2} & -\sqrt{2} & \sqrt{2} \end{array} \right) \left( \begin{array}{ccc} 4 & -1 & 1 \\ -1 & 4 & -1 \\ 1 & -1 & 4 \end{array} \right) \left( \begin{array}{ccc} \sqrt{3} & -1 & \sqrt{2} \\ \sqrt{3} & 1 & -\sqrt{2} \\ 0 & 2 & \sqrt{2} \end{array} \right) \\
&= \frac{1}{6} \left( \begin{array}{ccc} 3 \sqrt{3} & 3 \sqrt{3} & 0 \\ -3 & 3 & 6 \\ 6 \sqrt{2} & -6 \sqrt{2} & 6 \sqrt{2} \end{array} \right) \left( \begin{array}{ccc} \sqrt{3} & -1 & \sqrt{2} \\ \sqrt{3} & 1 & -\sqrt{2} \\ 0 & 2 & \sqrt{2} \end{array} \right) \\
&= \frac{1}{6} \left( \begin{array}{ccc} 18 & 0 & 0 \\ 0 & 18 & 0 \\ 0 & 0 & 36 \end{array} \right) = \left( \begin{array}{ccc} 3 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 6 \end{array} \right)
\end{align}
$$

重要例題$082$

行列式と置換⑦:置換行列(permutation matrix)

線形代数の枠組みで$n$次正方行列の行列式(determinant)を取り扱うにあたっては置換(permutation)という概念を抑えておく必要があります。当記事では置換行列(permutation matrix)について取りまとめを行いました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$4$章「行列式」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

置換行列

置換行列の定義

$K^{n}$の標準基底を$\{ \mathbf{e}_{1}, \mathbf{e}_{2}, \cdots , \mathbf{e}_{n} \}$のように定義する、このとき$1, \cdots , n$の$n$個の文字の置換$\sigma$について$n$次の置換行列$E(\sigma)$は下記のように定義される。
$$
\large
\begin{align}
E(\sigma) = \left[ \mathbf{e}_{\sigma(1)} \,\, \mathbf{e}_{\sigma(2)} \,\, \cdots \,\, \mathbf{e}_{\sigma(n)} \right]
\end{align}
$$

置換行列の性質

$n$次の置換行列$E(\sigma)$について下記が成立する。
$$
\large
\begin{align}
\det{(E(\sigma))} &= \mathrm{sgn}(\sigma) \quad (1) \\
E(\sigma \tau) &= E(\sigma) E(\tau) \quad (2) \\
\left[ E(\sigma) \right]^{-1} &= E \left( \sigma^{-1} \right) \quad (3) \\
AE(\sigma) &= \left[ \mathbf{v}_{\sigma(1)} \,\, \mathbf{v}_{\sigma(2)} \cdots \,\, \mathbf{v}_{\sigma(n)} \right] \quad (4)
\end{align}
$$

ただし上記の$\sigma, \tau$は$1, \cdots , n$の$n$個の置換、行列$A$は下記の$(5)$式の行列を表す。
$$
\large
\begin{align}
A &= \left[ \mathbf{v}_{1} \,\, \mathbf{v}_{2} \,\, \cdots \,\, \mathbf{v}_{n} \right] \quad (5)
\end{align}
$$

例題の確認

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

重要例題$043$

ハウスホルダー行列(Householder Matrix)の定義と成立する定理の導出

ハウスホルダー行列(Householder Matrix)は固有値解析にあたって行列の三重対角化を行う際などに用いられます。当記事ではハウスホルダー行列(Householder Matrix)の定義と成立する定理の導出などについて取りまとめました。
「A Suggestion on Mathematical Materials for Freshman Education Vol.44 ~ Householder Matrix ~」の内容を参考に作成を行いました。

・数学まとめ
https://www.hello-statisticians.com/explain-terms

ハウスホルダー行列の定義と定理

ハウスホルダー行列の定義

長さの等しい$2$つの$n$次元ベクトル$x \in \mathbb{R}^{n}$と$y \in \mathbb{R}^{n}$が与えられるとき、ハウスホルダー行列(Householder Matrix)を$H$とおくと$H$は下記のような式で定義されます。
$$
\large
\begin{align}
H &= I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{v^{\mathrm{T}} v} = I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{||v||^{2}} \\
v &= x \, – \, y \\
||x|| &= ||y||
\end{align}
$$

ハウスホルダー行列について成立する定理

前項「ハウスホルダー行列の定義」で定義したハウスホルダー行列については下記の定理がそれぞれ成立します。
・$(1) \,$ $Hx = y$
・$(2) \,$ $Hy = x$
・$(3) \,$ $H$は対称行列であり、$H^{\mathrm{T}} = H$が成立
・$(4) \,$ $H$は直交行列であり、$H^{\mathrm{T}} H = I_{n}$が成立
・$(5) \,$ $H$は固有値$1$を持つ
・$(6) \,$ $H$は固有値$-1$を持つ

定理の導出

$(1)$の導出

$v = x \, – \, y$、$||x|| = ||y||$より、$||v||^{2}$は下記のように表すことができます。
$$
\large
\begin{align}
||v||^{2} &= v^{\mathrm{T}} v \\
&= (x \, – \, y)^{\mathrm{T}} (x \, – \, y) \\
&= ||x||^{2} + ||y||^{2} – 2x \cdot y \\
&= ||x||^{2} + ||x||^{2} – 2x \cdot y \\
&= 2(||x||^{2} – x \cdot y)
\end{align}
$$

このとき$Hx$は下記のように式変形を行うことが可能です。
$$
\large
\begin{align}
Hx &= \left( I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{||v||^{2}} \right) x \\
&= x \, – \, \frac{2 (v^{\mathrm{T}} x) v}{||v||^{2}} \\
&= x \, – \, \frac{2(x \, – \, y)^{\mathrm{T}} x (x \, – \, y)}{||v||^{2}} \\
&= x \, – \, \frac{2(||x||^{2} \, – \, x \cdot y) (x \, – \, y)}{||v||^{2}} \\
&= x \, – \, \frac{\cancel{||v||^{2}} (x \, – \, y)}{\cancel{||v||^{2}}} \\
&= \cancel{x} \, – \, (\cancel{x} \, – \, y) \\
&= y
\end{align}
$$

上記より、$Hx = y$が成立することが確認できます。

$(2)$の導出

前項「$(1)$の導出」の式変形における$x$と$y$を入れ替えることで$Hy = x$が成立することが確認できます。

$(3)$の導出

$$
\large
\begin{align}
H = I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{||v||^{2}}
\end{align}
$$

上記を元に、$H^{\mathrm{T}}$は下記のように式変形を行うことができます。
$$
\large
\begin{align}
H^{\mathrm{T}} &= \left( I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{||v||^{2}} \right)^{\mathrm{T}} \\
&= I_{n}^{\mathrm{T}} \, – \, \frac{2(v v^{\mathrm{T}})^{\mathrm{T}}}{||v||^{2}} \\
&= I_{n} \, – \, \frac{2 (v^{\mathrm{T}})^{\mathrm{T}} v^{\mathrm{T}}}{||v||^{2}} \\
&= I_{n} \, – \, \frac{2 v v^{\mathrm{T}}}{||v||^{2}} \\
&= H
\end{align}
$$

上記より$H^{\mathrm{T}} = H$が成立するので、ハウスホルダー行列$H$が対称行列であることが確認できます。

$(4)$の導出

$H^{\mathrm{T}} = H$より、$H^{\mathrm{T}} H = H^{2}$が成立することに基づいて$H^{\mathrm{T}} H$は下記のように式変形を行うことができます。
$$
\large
\begin{align}
H^{\mathrm{T}} H &= H^{2} \\
&= \left( I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{||v||^{2}} \right)^{2} \\
&= I_{n} \, – \, \frac{4v v^{\mathrm{T}}}{||v||^{2}} + \frac{4(v v^{\mathrm{T}}) (v v^{\mathrm{T}})}{||v||^{4}} \\
&= I_{n} \, – \, \frac{4v v^{\mathrm{T}}}{||v||^{2}} + \frac{4v (v^{\mathrm{T}} v) v^{\mathrm{T}}}{||v||^{4}} \\
&= I_{n} \, – \, \frac{4v v^{\mathrm{T}}}{||v||^{2}} + \frac{4v ||v||^{2} v^{\mathrm{T}}}{||v||^{4}} \\
&= I_{n} \, – \, \frac{4v v^{\mathrm{T}}}{||v||^{2}} + \frac{4||v||^{2} v v^{\mathrm{T}}}{||v||^{4}} \\
&= I_{n} \, – \, \cancel{\frac{4v v^{\mathrm{T}}}{||v||^{2}}} + \cancel{\frac{4v v^{\mathrm{T}}}{||v||^{2}}} \\
&= I_{n}
\end{align}
$$

上記より$H^{\mathrm{T}} H = I_{n}$が成立するのでハウスホルダー行列$H$が直交行列であることが確認できます。

$(5)$の導出

$H(x+y)$は定理の$(1)$と$(2)$を用いて下記のように計算できます。
$$
\large
\begin{align}
H(x+y) &= Hx + Hy \\
&= y + x \\
&= (x+y)
\end{align}
$$

上記より$H$は固有値$1$を持つことが確認できます。

$(6)$の導出

$H(x \,- \, y)$は下記のように計算できます。
$$
\large
\begin{align}
H(x \, – \, y) &= \left( I_{n} \, – \, \frac{2v v^{\mathrm{T}}}{||v||^{2}} \right)v \\
&= v \, – \, \frac{2v (v^{\mathrm{T}}v)}{||v||^{2}} \\
&= v \, – \, \frac{2v \cancel{||v||^{2}}}{\cancel{||v||^{2}}} \\
&= v \, – \, 2v \\
&= – \, v \\
&= -(x \, – \, y)
\end{align}
$$

上記より$H$は固有値$-1$を持つことが確認できます。

二次形式(quadratic form)と対称行列(symmetric matrix)の対応

$n$個の変数についての$2$次の単項式$x_i, x_j$の実数係数の$1$次結合の式を$2$次形式といいます。当記事では二次形式(quadratic form)と対称行列(symmetric matrix)の対応について取りまとめを行いました。
作成にあたっては「チャート式シリーズ 大学教養 線形代数」の第$8$章「固有値問題と行列の対角化」を主に参考にしました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

$2$次形式と対称行列の対応

対称行列

行列$A$と$A$の転置行列$A^{\mathrm{T}}$について$A^{\mathrm{T}} = A$が成立するとき、$A$を対称行列という。下記に対称行列の具体例を表した。
$$
\large
\begin{align}
\left( \begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right), \, \left( \begin{array}{cc} 2 & 1 & 2 \\ 1 & 2 & 1 \\ 2 & 1 & 2 \end{array} \right)
\end{align}
$$

$2$次形式の定義

$n$個の変数$x_1, \cdots , x_{n}$についての$2$次の単項式$x_i, x_j$の実数係数の$1$次結合の式を$2$次形式という。$2$次形式を$q(x_1, \cdots , x_n)$とおくと、$q(x_1, \cdots , x_n)$は係数$a_{ij}$を用いて下記のように定義できる。
$$
\large
\begin{align}
q(x_1, \cdots , x_n) = \sum_{i=1}^{n} \sum_{j=1}^{n} a_{ij} x_{i} x_{j}
\end{align}
$$

ここで上記に対し、係数$a_{ij}$を$(i,j)$成分に持つ行列$A=(a_{ij})$、変数を成分に持つ列ベクトルを$\displaystyle \mathbf{x} = \left( \begin{array}{c} x_1 \\ \vdots \\ x_n \end{array} \right)$とおくと、$q(x_1, \cdots , x_n)$は下記のように表せる。
$$
\large
\begin{align}
q(x_1, \cdots , x_n) &= \sum_{i=1}^{n} \sum_{j=1}^{n} a_{ij} x_{i} x_{j} \\
&= \left( \begin{array}{ccc} x_1 & \cdots & x_n \end{array} \right) \left( \begin{array}{ccc} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \cdots & a_{nn} \end{array} \right) \left( \begin{array}{c} x_1 \\ \vdots \\ x_n \end{array} \right) \\
&= \mathbf{x}^{\mathrm{T}} A \mathbf{x} \quad (1)
\end{align}
$$

$2$次形式と対称行列の対応

行列$A$で決まる$2$次形式$q_{A}(x_1, \cdots , x_n)$の$x_{i}^{2}$の係数は$A$の対角成分$a_{ii}$に対応する。また、$x_{i}x_{j} = x_{j}x_{i}$に対応する係数は$a_{ij}+a_{ji}$であり、$a_{ij}, a_{ji}$は$\displaystyle \frac{a_{ij}+a_{ji}}{2}$で取り替えることもできる。よって、$2$次形式を前項の$(1)$のような表記に直すとき、$A$が対称行列であるという前提をおいても良い。

$2$次形式(quadratic form)と対称行列(symmetric matrix)の対応の例

以下、「チャート式シリーズ 大学教養 線形代数」の例題の確認を行う。

基本例題$161$

・$[1]$
$$
\large
\begin{align}
A = \left( \begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right)
\end{align}
$$

変数$x_1, x_2$についての$2$次形式を$q(x_1, x_2)$とおくと、$2$次形式の定義より、$q(x_1, x_2)$は下記のように求めることができる。
$$
\large
\begin{align}
q(x_1, x_2) &= \left( \begin{array}{cc} x_1 & x_2 \end{array} \right) \left( \begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \end{array} \right) \\
&= \left( \begin{array}{cc} x_1 & x_2 \end{array} \right) \left( \begin{array}{c} 2x_1 + x_2 \\ x_1 + 2x_2 \end{array} \right) \\
&= 2x_1^{2} + x_1 x_2 + x_1 x_2 + 2x_2^{2} = 2(x_1^{2} + x_1 x_2 + x_2^{2})
\end{align}
$$

・$[2]$
$$
\large
\begin{align}
A = \left( \begin{array}{ccc} 0 & 1 & 2 \\ 1 & 0 & 1 \\ 2 & 1 & 0 \end{array} \right)
\end{align}
$$

変数$x_1, x_2, x_3$についての$2$次形式を$q(x_1, x_2, x_3)$とおくと、$2$次形式の定義より$q(x_1, x_2, x_3)$は下記のように求めることができる。
$$
\large
\begin{align}
q(x_1, x_2, x_3) &= \left( \begin{array}{ccc} x_1 & x_2 & x_3 \end{array} \right) \left( \begin{array}{ccc} 0 & 1 & 2 \\ 1 & 0 & 1 \\ 2 & 1 & 0 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) \\
&= x_1 x_2 + 2 x_1 x_3 + x_2 x_1 + x_2 x_3 + 2 x_3 x_1 + x_3 x_2 \\
&= 2x_1x_2 + 4x_1 x_3 + 2 x_2 x_3
\end{align}
$$

$[1]$ではベクトル・行列の積に基づいて計算を行ったが、$[2]$では前節の$(1)$式に基づいて行列の係数のインデックスに基づいて各単項式の書き下しを行なった。

基本例題$162$