見出し画像

6-3 推定量の分散 ～ ばらつきを小さくする重さの計測方法

ネイピア DS

今回の統計トピック

「推定量」の概念と「推定量の分散」について考えます。
重さの計測・推定方法を工夫することで「重さの推定時の誤差（分散）を抑制できる」ことも学びます！

公式問題集の準備

「公式問題集」の問題を利用します。お手元に公式問題集をご用意ください。
公式問題集が無い場合もご安心ください！
「知る」「実践する」の章で、のんびり統計をお楽しみください！

問題を解く

📘公式問題集のカテゴリ

標本分布の分野
問3 推定量の分散（天秤はかりの分散）

試験実施年月
調査中

問題

公式問題集をご参照ください。

解き方

題意
$${(X+Y)/2, \ (X-Y)/2}$$の意味をクリアにして、シンプルに確率変数の分散を求めます。

コイン$${A}$$の重さは$${a}$$、コイン$${B}$$の重さは$${b}$$です。
天秤はかりでコインの重さを次の２通りの方法で計測します。
①片方の受け皿にコイン$${A,B}$$を置き、重さ$${X=a+b}$$を測る方法
②両方の受け皿に別々にコイン$${A,B}$$を置き、重さ$${Y=a-b}$$を測る方法

計測は誤差を伴っています。
誤差は平均$${0}$$、分散$${\sigma^2}$$の独立な確率変数$${\varepsilon_1,\ \varepsilon_2}$$です。
$${X=a+b+\varepsilon_1}$$、$${Y=a-b+\varepsilon_2}$$になります。

$${(X+Y)/2}$$と$${(X-Y)/2}$$でコイン$${A,\ B}$$の重さを推定できます。
「コイン$${B}$$の重さの推定量の分散」を解きましょう。

公式問題集の記述を改変

まずは、$${X=a+b+\varepsilon_1}$$、$${Y=a-b+\varepsilon_2}$$の図です。
計測誤差の$${\varepsilon_1,\ \varepsilon_2}$$は、ひとまず外して図示します。

この図で$${X}$$・$${Y}$$と釣り合っているコイン$${A,\ B}$$を足したり引いたりして、$${(X+Y)/2}$$と$${(X-Y)/2}$$の秘密を探りましょう！

$${\boldsymbol{(X+Y)/2}}$$と$${\boldsymbol{(X-Y)/2}}$$って？

$${(X+Y)/2}$$でコイン$${A}$$の重さを推定できます。

$${(X-Y)/2}$$でコイン$${B}$$の重さを推定できます。

$${X=a+b+\varepsilon_1}$$と$${Y=a-b+\varepsilon_2}$$を代入してみましょう。

$${\boldsymbol{(X+Y)/2}}$$とコイン$${\boldsymbol{A}}$$の重さ$${\boldsymbol{a}}$$

$$
\begin{align*}
(X+Y)/2&=(a+b+\varepsilon_1 + a-b+\varepsilon_2)/2 \\
&=(2a+\varepsilon_1+\varepsilon_2)/2 \\
&=a+(\varepsilon_1+\varepsilon_2)/2 \\
\end{align*}
$$

コイン$${A}$$の重さ$${a}$$に誤差$${(\varepsilon_1+\varepsilon_2)/2}$$を加えた値を測定できます。

$${\boldsymbol{(X-Y)/2}}$$とコイン$${\boldsymbol{B}}$$の重さ$${\boldsymbol{b}}$$

$$
\begin{align*}
(X-Y)/2&=(a+b+\varepsilon_1 - a+b-\varepsilon_2)/2 \\
&=(2b+\varepsilon_1-\varepsilon_2)/2 \\
&=b+(\varepsilon_1-\varepsilon_2)/2 \\
\end{align*}
$$

コイン$${B}$$の重さ$${b}$$に誤差$${(\varepsilon_1-\varepsilon_2)/2}$$を加えた値を測定できます。

天秤座のイラスト（星座）：「いらすとや」さんより

コイン$${\boldsymbol{B}}$$の重さの推定量は確率変数

コイン$${B}$$の重さの推定量は$${(X-Y)/2=b+(\varepsilon_1-\varepsilon_2)/2}$$です。
確率変数$${\varepsilon_1,\ \varepsilon_2}$$を$${2}$$で割って定数$${b}$$を足している「$${b+(\varepsilon_1-\varepsilon_2)/2}$$は確率変数」です。
確率変数$${b+(\varepsilon_1-\varepsilon_2)/2}$$の分散$${V[b+(\varepsilon_1-\varepsilon_2)/2]}$$を確認しましょう。

分散$${\boldsymbol{V[b+(\varepsilon_1-\varepsilon_2)/2]}}$$を解く  

分散$${V[b+(\varepsilon_1-\varepsilon_2)/2]}$$の$${[ \ ]}$$の中をほぐします。

$$
\begin{align*}
V[b+(\varepsilon_1-\varepsilon_2)/2]&=V[ \ b+(\varepsilon_1-\varepsilon_2)/2 \ ] \\
&=V \left[ \cfrac{1}{2} \ (\varepsilon_1-\varepsilon_2)+b \right] \\
&=V \left[ \cfrac{1}{2} \ \varepsilon_1-\cfrac{1}{2} \ \varepsilon_2+b \right] \\
\end{align*}
$$

ここで分散の演算の公式の登場です。

確率変数$${X, Y}$$が独立のとき、
$${V[aX+bY+c]=a^2V[X]+b^2V[X]}$$

$${\varepsilon_1,\ \varepsilon_2}$$は独立な確率変数です。
$${1/2}$$は確率変数$${\varepsilon_1}$$に掛かる係数です。
$${-1/2}$$は確率変数$${\varepsilon_2}$$に掛かる係数です。
$${b}$$は定数です。

$${V[aX+bY+c]=a^2V[X]+b^2V[Y]}$$の公式に「$${X=\varepsilon_1,\ Y=\varepsilon_2,\ a=1/2,\ b=-1/2,\ c=b}$$」を当てはめます。

$$
\begin{align*}
V[b+(\varepsilon_1-\varepsilon_2)/2]&=V \left[ \cfrac{1}{2} \ \varepsilon_1-\cfrac{1}{2} \ \varepsilon_2+b \right] \\
 \\
&=\ V \left[ \cfrac{1}{2} \varepsilon_1 \right]+ \ V \left[-\cfrac{1}{2}\varepsilon_2 \right] +V[b] \\
  \\
&= \left(\cfrac{1}{2} \right)^2\ V[\varepsilon_1]+ \left(-\cfrac{1}{2} \right)^2 \ V[\varepsilon_2] +V[b] \\
  \\
&=\cfrac{1}{4} \ V[\varepsilon_1]+\cfrac{1}{4} \ V[\varepsilon_2] \\
  \\
&=\cfrac{1}{4}  \sigma^2+ \cfrac{1}{4}  \sigma^2 \\
  \\
&=\cfrac{\sigma^2}{2}
\end{align*}
$$

出題の「コイン$${B}$$の重さの推定量の分散」は$${\cfrac{\sigma^2}{2}}$$です。

勉強で使ういろいろなマーク：「いらすとや」さんより

コイン$${\boldsymbol{A}}$$の重さの推定量の分散 
せっかくですので、コイン$${A}$$も計算しましょう。
コイン$${A}$$の重さの推定量は$${(X+Y)/2=a+(\varepsilon_1+\varepsilon_2)/2}$$です。
分散$${V[a+(\varepsilon_1+\varepsilon_2)/2]}$$をほぐします。

$$
\begin{align*}
V[a+(\varepsilon_1+\varepsilon_2)/2]&=V \left [\cfrac{1}{2}\ \varepsilon_1 + \cfrac{1}{2}\ \varepsilon_2 + a\right]\\
 \\
&=\ V \left[ \cfrac{1}{2} \varepsilon_1 \right]+ \ V \left[\cfrac{1}{2}\varepsilon_2 \right] +V[a] \\
  \\
&= \left(\cfrac{1}{2} \right)^2\ V[\varepsilon_1]+ \left(\cfrac{1}{2} \right)^2 \ V[\varepsilon_2] +V[a] \\
 \\
&=\cfrac{1}{4} \ V[ \varepsilon_1] +\cfrac{1}{4} \ V[ \varepsilon_2]  \\
  \\
&=\cfrac{1}{4} \ \sigma^2 +\cfrac{1}{4} \ \sigma^2  \\
  \\
&=\cfrac{\sigma^2}{2}
\end{align*}
$$

コイン$${A}$$の重さの推定量の分散は$${\cfrac{\sigma^2}{2}}$$です。
コイン$${B}$$の重さの推定量の分散と同じになりました。


勉強で使ういろいろなマーク：「いらすとや」さんより

問題の文章をパッと見たときに「何だか難しそうだ」という印象を受けました。
しかし、丁寧に問題文章を追いかけて、次の４つの流れを見つけると、サクッと解答できました。
① コインＢの重さの推定式は$${\cfrac{X-Y}{2}}$$です。
② 分散$${\sigma^2}$$が出現する誤差$${\varepsilon_1, \varepsilon_2}$$を含んだ$${X,Y}$$の重さの計算式を見出します。
③ 式①に式②を代入します。
④ 式③を分散の線形結合の公式に当てはめて解きます。

解答

③ $${\cfrac{\sigma^2}{2}}$$ です。

難易度　ふつう

・知識：分散
・計算力：数式組み立て（中）
・時間目安：1分

知る

おしながき

公式問題集の問題に接近してみましょう！
今回は、推定量の概要を確認して、公式問題集の計測方法 $${(X+Y)/2,\ (X-Y)/2}$$ の秘密に迫ります。

推定量

📕公式テキスト：3.3.1 点推定（106ページ～）
参考書籍「統計学入門」による推定量の説明を要約して記載いたします。

母集団の確率分布を決めている定数を母数（parameter）と呼びます。
母数の例は、母平均$${\mu}$$、母分散$${\sigma^2}$$です。

母数を推定するために標本から求めた統計量を推定量（estimator）と呼びます。
例えば、標本平均$${\bar{X}}$$は母平均$${\mu}$$の推定量、標本分散$${s^2}$$は母分散$${\sigma^2}$$の推定量です。

推定量には$${\hat{}}$$を付けます。
例えば、母平均の推定量は$${\hat{\mu}}$$、母分散の推定量は$${\hat{\sigma^2}}$$です。

具体的な観測値にもとづいて計算される値を推定値（estimate）と呼びます。

参考書籍「統計学入門」第11章の記述を改変

重さの推定量

推定量は確率変数です。
したがって、推定量にも統計量としての期待値や分散があります。

公式問題集の問題では、$${(X+Y)/2, (X-Y)/2}$$をコイン$${A,\ B}$$の重さ$${a,\ b}$$の推定量とし、これらの分散を求めていました。
振り返ります。

まず次の２つの重さを得ます。

・コイン$${A, B}$$を天秤の同じ皿に乗せて計測して得た重さ$${X}$$
・コイン$${A, B}$$を天秤の別々の皿に乗せて計測して得た重さ$${Y}$$

そして、次の計算でコインの重さの推定量を得ます。

・コイン$${A}$$の重さ$${a}$$の推定量は$${(X+Y)/2}$$
・コイン$${B}$$の重さ$${a}$$の推定量は$${(X-Y)/2}$$

このような推定によって、計測誤差の分散を$${\sigma^2/2}$$に抑えています。

コインをバラで計測する
では、コイン$${A}$$だけを天秤に乗せて計測し、コイン$${B}$$だけを天秤に乗せて計測する場合の分散はどうなるでしょう？
計測した$${A}$$の推定量を$${a+\varepsilon_3}$$、$${B}$$の推定量を$${b+\varepsilon_4}$$としましょう。
$${\varepsilon_3,\ \varepsilon_4}$$は誤差項であり、平均$${\mu}$$、分散$${\sigma^2}$$の独立な確率変数です。

コイン$${\boldsymbol{A}}$$の重さの推定量の分散 
$${V[a+\varepsilon_3]}$$をほぐしましょう。

$${V[a+\varepsilon_3]=V[\varepsilon_3]=\sigma^2}$$
コイン$${A}$$の重さの推定量の分散は$${\sigma^2}$$です。

コイン$${\boldsymbol{B}}$$の重さの推定量の分散 
$${V[b+\varepsilon_4]}$$をほぐしましょう。

$${V[b+\varepsilon_4]=V[\varepsilon_4]=\sigma^2}$$
コイン$${B}$$の重さの推定量の分散は$${\sigma^2}$$です。

コインをバラで計測するときの分散は$${\sigma^2}$$になりました。

なんと！
公式問題集の計測方法による推定量$${(X+Y)/2,\ (X-Y)/2}$$の分散$${\sigma^2/2}$$の方が値が小さいです。
つまり、公式問題集の計測方法のほうがバラツキが少ないということであり、推定の精度が高いというのです。
計測方法によって計測値（推定値）の分散が変わるのって、面白いですね！

ということで、公式問題集の推定方法と、コインをバラで計測する方法について、バラツキ方をPythonで実験してみましょう！

保護具を付けて実験をする人のイラスト（男性）：「いらすとや」さんより

重さの推定量の分散の実験

■重さの計測方法
次の４つの方法でそれぞれ５０回、重さを計測します。

  1. 重さ$${X=a+b+\varepsilon_1}$$：コイン$${A, B}$$を天秤の同じ皿に乗せて計測（公式問題集の計測方法）

  2. 重さ$${Y=a-b+\varepsilon_2}$$：コイン$${A, B}$$を天秤の別々の皿に乗せて計測（公式問題集の計測方法）

  3. コイン$${A}$$をバラで計測 $${a+\varepsilon_3}$$

  4. コイン$${B}$$をバラで計測 $${b+\varepsilon_4}$$

天秤のイラスト：「いらすとや」さんより

■重さの推定方法
次の４つの方法でコインの重さの推定します。

  1. コイン$${A}$$の重さ$${a}$$の推定値1（公式問題集の推定方法）
    $${(X+Y)/2=a+(\varepsilon_1+\varepsilon_2)/2}$$

  2. コイン$${A}$$の重さ$${a}$$の推定値2（コインをバラで計測）
    $${a+\varepsilon_3}$$

  3. コイン$${B}$$の重さ$${b}$$の推定値1（公式問題集の推定方法）
    $${(X-Y)/2=b+(\varepsilon_1-\varepsilon_3)/2}$$

  4. コイン$${B}$$の重さ$${b}$$の推定値2（コインをバラで計測）
    $${b+\varepsilon_4}$$

予想屋のイラスト：「いらすとや」さんより

■設定値
コインの真の重さを$${a=10,\ b=8}$$にします。
誤差の$${ \varepsilon_1, \varepsilon_2, \varepsilon_3, \varepsilon_4}$$は、連続一様分布$${U(-3,\ 3)}$$の乱数で取得します。
連続一様分布$${U(-3,\ 3)}$$の期待値は$${0}$$、分散は$${3}$$です。

実験の結果
🔵青丸が公式問題集の方法、🔴赤丸がバラで計測する方法です。

んー、バラツキの差がぼんやりしています。
🔵青丸と🔴赤丸の見た目のバラツキは、バラで計測する🔴赤丸の方がバラついている感じがしないでもない、こんな感じです。
実は、乱数シードの値を変えると、バラで計測する🔴赤丸の方が重さの真値に近くなることもあります。

試行回数を 10回 にしてみました。
こちらのほうが、バラで計測する🔴赤丸のバラツキが大きい感じがします。

実験は以上です。
少しモヤモヤが残りました。
Pythonのコードは「Pythonで作成してみよう！」の章で紹介いたします。

アイディアの共有のイラスト：「いらすとや」さんより


実践する

重さの推定量の分散を実験してみよう

「知る」の重さの推定量の分散を実験しましょう。
Pythonを活用します。


電卓・手作業で作成してみよう！

「知る」の章の推定量の分散の計算を実施しましょう！


EXCELで作成してみよう！

今回はお休みです。

EXCELサンプルファイルのダウンロード
今回はファイル提供はありません。


Pythonで作成してみよう！

プログラムコードを読んで、データを流したりデータを変えてみたりして、データを追いかけることで、作表ロジックを把握する方法も効果的でしょう。
サンプルコードを揃えておけば、類似する作表作業を自動化して素早く結果を得ることができます。

今回は、「知る」の章で用いた「重さの推定量の分散」実験のPythonコードです。

①インポート
誤差の$${ \varepsilon_1, \varepsilon_2, \varepsilon_3, \varepsilon_4}$$の乱数取得に、scipy.statsのuniform（連続一様分布）を利用します。

import numpy as np
from scipy.stats import uniform
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'MS Gothic'
%matplotlib inline

②重さの推定量の分散の実験
コードを細かく分解して書きます。

■設定
次の値を変更することで、さまざまなケースの実験結果を得ることができます。

  • a, b はコインの重さです。

  • low, high は連続一様分布$${U(low, high)}$$のパラメータである区間です。
    例では、-3 から 3 の乱数を生成します。

  • size は標本サイズです。ｎです。

  • np.random.seed( ) は乱数シードです。

### 設定 
a, b = 10, 8         # コインA,Bの真の重さa,b
low, high = -3, 3    # 一様分布のパラメータ U(low, high)
size= 10             # 標本サイズ
np.random.seed(35)   # 乱数シード # size=10→35, size=50→10

■重さの計測の試行と重さの推定値の計算

  • X, Y は公式問題集の方法で計測した重さです。
    uniform.rvsで生成した -3 から 3 の範囲の連続一様分布乱数を取得して、誤差epsilon1, epsilon2に設定しています。

  • a_est1, b_est1 は公式問題集の推定方法$${(X+Y)/2,\ (X-Y)/2}$$による重さの推定値です。

  • a_est2, b_est2 はバラで計測するときの推定値$${a+\varepsilon_3, b+\varepsilon4}$$です。
    uniform.rvsで生成した -3 から 3 の範囲の連続一様分布乱数を取得して、誤差epsilon3, epsilon4に設定しています。

### 重さ計測の試行

## 公式問題集の方法による計測
#  X,Yの重さの計測 X=a+b+ε1, Y=a-b+ε2
X = np.array([a + b + epsilon1 for epsilon1 in uniform.rvs(low, high-low, size)])
Y = np.array([a - b + epsilon2 for epsilon2 in uniform.rvs(low, high-low, size)])

#  コインの重さの推定量 a_est1=(X+Y)/2, b_est1=(X-Y)/2
a_est1 = (X + Y) / 2
b_est1 = (X - Y) / 2

## コインをバラで計測による推定量 a_est2=a+ε3, b_est2=b+ε4
a_est2 = np.array([a + epsilon3 for epsilon3 in uniform.rvs(low, high-low, size)])
b_est2 = np.array([b + epsilon4 for epsilon4 in uniform.rvs(low, high-low, size)])

■グラフのプロット
①プロットの設定
グラフの領域 fig, ax の設定、散布図の x 軸の値、グラフの線幅・点の大きさ・透明度、グラフの点の色を設定します。

### プロット
## プロットの設定
fig, ax = plt.subplots(1, 2, figsize=(12,5))
x = np.arange(1, size+1, 1)
linewidth, dotsize, alpha = 1, 200, 0.3   # 線幅、点の大きさ、透明度
est1_color, est2_color, true_color = 'blue', 'red', 'black'

②コイン$${A}$$の重さのプロット
次の５つの要素をプロットしています。

  1. 重さ$${a}$$を公式問題集の方法で計測して推定するときの推定値$${(X+Y)/2}$$の散布図

  2. 重さ$${a}$$をバラで計測して推定するときの推定値$${a+\varepsilon_3}$$の散布図

  3. 重さ$${a}$$の真値の水平線

  4. 番号 1 の「公式問題集の方法による重さ$${a}$$の推定値の平均値」の水平線

  5. 番号 2 の「バラで計測する方法による重さ$${a}$$の推定値の平均値」の水平線

## コインAの重さのプロット
# a_est1=(X+Y)/2=a+(epsilon1+epsilon2)/2の標本の散布図
ax[0].scatter(x, a_est1, c=est1_color, s=dotsize, alpha=alpha,
              label='$a$の推定値 $a+(\\varepsilon_1+\\varepsilon_2)/2$')
# a_est2=a+epsilon3の散布図の標本の散布図
ax[0].scatter(x, a_est2, c=est2_color, s=dotsize, alpha=alpha,
              label='$a$の推定値 $a+\\varepsilon_3$')
# 重さaの真値の水平線
ax[0].hlines(a, 0, size, lw=linewidth, ec=true_color, label=f'$a$の真値:{a}')
# 標本平均の水平線　(a_est1=(X+Y)/2=a+(epsilon1+epsilon2)/2)/size
ax[0].hlines(np.mean(a_est1), 0, size, lw=linewidth, ec=est1_color,
             label=f'$(X+Y)/2$の平均:{np.mean(a_est1):.2f}, '
                   f'$\Delta$ {np.mean(a_est1)-a:.2f}')
# 標本平均の水平線　(a_est2=a+epsilon3)/size
ax[0].hlines(np.mean(a_est2), 0, size, lw=linewidth, ec=est2_color,
             label=f'$a+\\varepsilon_3$の平均:{np.mean(a_est2):.2f}, '
                   f'$\Delta$ {np.mean(a_est2)-a:.2f}')
ax[0].legend(loc='best')

③コイン$${B}$$の重さのプロットと描画
コイン$${A}$$と同様にコイン$${B}$$の重さをプロットします。

## コインBの重さのプロット
# b_est1=(X-Y)/2=b+(epsilon1-epsilon2)/2の標本の散布図
ax[1].scatter(x, b_est1, c=est1_color, s=dotsize, alpha=alpha,
              label='$b$の推定値 $b+(\\varepsilon_1-\\varepsilon_2)/2$')
# b_est2=b+epsilon4の散布図の標本の散布図
ax[1].scatter(x, b_est2, c=est2_color, s=dotsize, alpha=alpha, 
              label='$b$の推定値 $b+\\varepsilon_4$')
# 重さbの真値の水平線
ax[1].hlines(b, 0, size, lw=linewidth, ec=true_color, label=f'$b$の真値:{b}')
# 標本平均の水平線　(b_est1=(X-Y)/2=b+(epsilon1-epsilon2)/2)/size
ax[1].hlines(np.mean(b_est1), 0, size, lw=linewidth, ec=est1_color,
             label=f'$(X-Y)/2$の平均:{np.mean(b_est1):.2f}, '
                   f'$\Delta$ {np.mean(b_est1)-b:.2f}')
# 標本平均の水平線　(b_est2=b+epsilon4)/size
ax[1].hlines(np.mean(b_est2), 0, size, lw=linewidth, ec=est2_color,
             label=f'$b+\\varepsilon_4$の平均:{np.mean(b_est2):.2f}, '
                   f'$\Delta$ {np.mean(b_est2)-b:.2f}')
ax[1].legend(loc='best')

plt.tight_layout()
plt.show()
出力イメージ

Pythonサンプルファイルのダウンロード
こちらのリンクからJupyter Notebook形式のサンプルファイルをダウンロードできます。


おわりに

今回の実験～シミュレートはすっきりしない結果でした。
気分を取り直して、また記事執筆に励みます。

最後までお読みいただきまして、ありがとうございました。

のんびり統計シリーズの記事

次の記事

前の記事

目次

この記事が気に入ったらサポートをしてみませんか？