線形モデルのうち、最も簡単な例のひとつとして、ある物体(一つしかない)の長さ$\mu$の計測を考えます。
y = \mu + \epsilon
ある物体(一つしかない)の長さ$\mu$は一定ですが、操作あるいは観測などによって誤差$\epsilon$が生じます。実験においては、物体に応じてノギスやマイクロメータ、メジャーなどにより複数回測定を行い、平均と不確かさを算出します。
例
物体(一つしかない)の長さを3回計測した結果が、2.25 mm、2.27 mm、2.24 mmとすると、平均は2.25 mm(有効数字3桁)、標本標準偏差は0.015 mmとなります。
よって、標本標準偏差を不確かさとして採用する場合、物体の長さの推定値は2.25 ± 0.02 mmとなります。
s=\sqrt{{((2.25-2.25333)^2+(2.27-2.25333)^2+(2.24-2.25333)^2})/(3-1)}=0.015275
なお、同じ 3 点のデータから
u=\frac{s}{\sqrt{n}}=0.009
として 平均値の標準不確かさを求めることもできます。これは「母平均の推定誤差」を表す量であり、複数個の物体の平均的な長さを推定する場合(工場のロット検査など)に適しています。しかし本例では、対象は 単一の物体であり、推定したいのはその物体固有の長さ$\mu$であるため、測定値のばらつきを表す標本標準偏差(0.015 mm)を不確かさとして採用するのが適切です。
個人的には、±(数値)のような表現をした時にはそれが何を計算したものなのか書いてくれると嬉しいです。
学生時代に物理学実験において、一つの物体の長さ(や重さ)がある物理量Xを文字で表現したときのパラメータの一つとなっており、ある物理量Xの不確かさを誤差伝播(不確かさの合成)で求める際に、物体の長さ(や重さ)の不確かさが必要となったので、この例はそういったとき参考になると思います。
本題に戻ります。理想的な状況ではなく、現実を反映させたシミュレーションをする際、誤差、ノイズ等を考慮するのはよくあることです。計測を独立に$n$回行うとして、$i$回目の計測結果を
y_i = \mu + \epsilon_i, i=1,\dots,n
と表します。今回はこれをコンピュータでどんなふうに表現するか、考えてみましょう。
シミュレーションデータの作成
問題は自然科学の統計学75ページ 2.1 i)繰り返し測定のモデルを参考にしています。
[0,1)の一様分布(離散ではなく連続)から乱数を生成し、$u_1、u_2、u_3、\dots$として、$\epsilon_i(i=1,2,3,\dots)$を
\epsilon_1=\sum_{i=1}^{12} u_i-6,\epsilon_2=\sum_{i=13}^{24} u_i-6,\epsilon_3=\sum_{i=25}^{36} u_i-6,\dots
のように定義して、次の確率変数の列$y_i$を計算しましょう。
y_i = 22.5 + \epsilon_i\
(i=1,2,3,...,10)
# 75ページ 2.1 i)繰り返し測定のモデル を参考にした!!
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
block = 12 # 1区間の長さ
N = 10 # y1〜y10 を求める
u = np.random.rand(block * N)
baseline = 22.5
ys = []
for k in range(N):
ep = np.sum(u[k*block:(k+1)*block]) - 6
y = baseline + ep
ys.append(y)
# 結果の表示
#for i, y in enumerate(ys, start=1):
#print(f"y{i}: {y:.5f}")
# --- ヒストグラム作成 ---
plt.figure(figsize=(6,4))
plt.hist(ys, bins=5, edgecolor='black', alpha=0.7)
plt.xlabel("y value")
plt.ylabel("Frequency")
plt.title("Histogram of generated y values")
plt.grid(alpha=0.3)
plt.show()
簡単にコードを解説します。
ライブラリの話
numpyは行列や数値計算で大変便利なライブラリです。FFTや乱数生成もできて便利です。matplotlibはグラフの表示に大変便利なライブラリです。
乱数のシードの話
乱数のシードを指定しなくてもプログラムは動きますが、シミュレーションにおいて再現性を確保するためにシードを指定しています。例えば、画像処理のプログラムを軽量に美しく書き換えて、書き換え前後で結果(画像の行列)を比較したときに、乱数(画像処理では必要に応じてノイズなどを入れます)が実行のたびに違っていたら、自分のコードの書き換えが間違っているのか、乱数の違いのせいなのか判別ができません。
yのデータ数Nの話
問題ではN=10ですが、N=1000とかにすると、ヒストグラムが正規分布に近い形になって面白いです。