背景
統計検定２級の問題で、分散を計算し、その大小を比較させる問題がありました。
（統計検定２級受験予定で、この問題を解いてない方は、答えを見ないでくださいね。）
$ X_{1},X_{2},,,,X_{n} $ （n>=3 の整数）が独立に平均μ、分散 $ σ^2 $ の分布に従うとき、
\begin{align}
μ1 &= \frac{1}{n}\sum_{i=1}^{n}X_{i} \\
μ2 &= \frac{2}{n(n+1)}\sum_{i=1}^{n}iX_{i}
\end{align}
の分散はどちらが小さいか、という問題です。
（オリジナル問題とは違いますが、ポイントのみ抽出して、問題を単純化しています。）
μ1 は、単純な平均。
μ2 は、i番目には、iを掛けているので、iが大きくなるほど値が大きくなるのでばらつきが大きくなりそうだが、n(n+1)と$ n^2 $ で割っているので、、、どうなるのやら、
難しい問題ですね。
ちなみに、期待値は、どちらも、μになります。
これを、プログラムで検証してみました。
準備
まずは、準備
- numpy のインポート
- n個のXを出力する関数（正規分布をランダム出力）
- μ1、μ2 を計算する関数
- 平均、標準偏差、サンプルサイズn（３以上の整数）を定義
# 準備
import numpy as np
rng = np.random.default_rng()
# 平均μ、標準偏差σ の正規分布に従う、サイズnのサンプルデータを取得
def get_X( mu, sigma, n ):
return rng.normal(mu,sigma,n)
# μ1の計算
def calc_mu1(X):
return X.mean()
# μ2の計算
def calc_mu2(X):
mu2 = 0
for i,x in enumerate(X):
mu2 += (i+1) * x
n = X.shape[0]
mu2 *= 2/(n*(n+1))
return mu2
mu = 10 # 平均10
sigma = 3 # 標準偏差3
n_list = range(3,101) # nは3以上の整数
平均と分散を計算
次に、平均と分散を計算するために、
- n個のXを出力し、μ1とμ2を計算する。
- この処理を1000回繰り返し、リスト化しておく。
- そのリストから、μ1とμ2の平均と分散を求める。
- この処理を、3以上の整数nについて繰り返す。
# 平均のリスト
mu1_mean_list = list()
mu2_mean_list = list()
# 分散のリスト
mu1_var_list = list()
mu2_var_list = list()
# 3以上の整数nのループ
for n in n_list:
mu1_list = list()
mu2_list = list()
# 1000回繰り返す
for i in range(1000):
# Xの取得
X = get_X(mu,sigma,n)
# μ1の計算
mu1_list.append(calc_mu1(X))
# μ2の計算
mu2_list.append(calc_mu2(X))
# μ1の平均/分散
mu1_mean_list.append(np.array(mu1_list).mean())
mu1_var_list.append(np.array(mu1_list).var())
# μ2の平均
mu2_mean_list.append(np.array(mu2_list).mean())
mu2_var_list.append(np.array(mu2_list).var())
平均と分散の可視化
平均/分散のリストを単純に出力しても、数字の羅列ではわからないので可視化する。
- μ1 と μ2 の平均、もしくは分散をプロットする関数を定義
- 理論値が簡単に（直点的に）プロットできる場合は、理論値も
import matplotlib.pyplot as plt
import japanize_matplotlib
# μ1 と μ2 の平均、もしくは分散をプロットする関数
def plot_mu1_mu2(title,list1,list2,label1='μ1',label2='μ2', y_line=None, line_label=None):
plt.plot(n_list,list1, label=label1)
plt.plot(n_list,list2, label=label2)
# 理論値を直線表示する場合
if y_line!=None:
plt.plot([n_list[0],n_list[-1]],[y_line,y_line], label=line_label,color='red')
plt.title(title)
plt.xlabel('n')
plt.legend()
plt.show()
# 平均の可視化
plot_mu1_mu2('平均',mu1_mean_list,mu2_mean_list,y_line=mu,line_label='理論値')
平均は、理論値の10に収束しているように見える。
続いて、分散の可視化
# 分散の可視化
plot_mu1_mu2('分散',mu1_var_list,mu2_var_list)
問題が、「μ1とμ2の分散はどちらが小さいか」なので、μ1<μ2 のように見える。
理論値を求める
ここで、ズバリ回答（理論値）を求めると、
\begin{align}
V[μ1] &= \frac{1}{n}σ \\
V[μ2] &= \frac{2(2n+1)}{3n(n+1)}σ
\end{align}
となります。
自然数の２乗和
自然数の２乗和が計算できないと解答できない、超難しい問題です。（統計検定２級の範囲なんですかね、ここまでの数学知識？）
\begin{align}
\sum_{i=1}^{n}X_{i}^2 = \frac{n(n+1)(2n+1)}{6}
\end{align}
分散の比
単純に分散を可視化しただけではわかりにくいので、違いが分かるように、分散の比 $ \frac{V[μ2]}{V[μ1]} $ を可視化します。
分散の比の理論値は、
\begin{align}
\frac{V[μ2]}{V[μ1]} = \frac{\frac{2(2n+1)}{3n(n+1)}σ}{\frac{1}{n}σ} = \frac{2(2n+1)}{3(n+1)}
\end{align}
となります。この理論値の極限も、
\begin{align}
\lim_{n \to \infty} \frac{V[μ2]}{V[μ1]} = \frac{4+\frac{2}{n}}{3+\frac{3}{n}} = \frac{4}{3}
\end{align}
となるので、一緒に可視化します。
# 分散の比 V[μ2]/V[μ1]を求める
var_ratio = np.array(mu2_var_list) / np.array(mu1_var_list)
# 理論値
theoretical_value = [ 2*(2*n+1) / (3*(n+1)) for n in n_list]
# 極限値
limit_value = 4/3
plot_mu1_mu2('分散の比 V[μ2]/V[μ1]',var_ratio,theoretical_value,'比','理論値', limit_value, '極限値')
多少のばらつきはあるものの、理論値に近づいており、nが十分大きくなると、極限値になるようです。
まとめ
理論値通り、V[μ1] < V[μ2] となり、その比は、$ \frac{4}{3} $(1.3333)に収束していく、ということがプログラムでも検証できました。