ホーム » 統計検定1級 2014年 統計数理

「統計検定1級 2014年 統計数理」カテゴリーアーカイブ

投稿一覧

2014 Q5(4)

投稿日:2024年7月24日 　最終更新日:2025年1月18日

テレビ番組の視聴満足回数の調査で満足と答えた回数のそれぞれ人数の適合度検定を実施し適合を判断する手順を述べました。

　

コード

多項確率q_jが、q_j = q_j(p) = 4C_j p^j (1 - p)^{4-j}, \quad j = 0, \ldots, 4に忠実に基づく場合と、各番組または各回答者の満足度が異なるモデルにおいて、カイ二乗統計量、尤度比統計量、およびp-値を計算し、モデルの適合性を比較します。

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

# シミュレーションの設定
n = 100  # 回答者数
p_true = 0.6  # 真の成功確率（忠実なモデル）
categories = np.arange(5)  # カテゴリ (0～4)
n_simulations = 1000  # シミュレーション回数

# 忠実なモデルの確率 q_j
qj_true = [np.math.comb(4, j) * (p_true ** j) * ((1 - p_true) ** (4 - j)) for j in categories]

# 統計量を記録するリスト
chi_squared_true = []
likelihood_ratio_true = []
chi_squared_noisy_program = []
likelihood_ratio_noisy_program = []
chi_squared_noisy_responder = []
likelihood_ratio_noisy_responder = []

# 二項分布に従わない (qj_noisy_program): q_j をランダムに変更してズラす
qj_noisy_program = np.array(qj_true) + np.random.uniform(-0.1, 0.1, size=len(qj_true))
qj_noisy_program = np.clip(qj_noisy_program, 0, None)  # 確率が0未満にならないようにする
qj_noisy_program = qj_noisy_program / np.sum(qj_noisy_program)  # 確率なので正規化

# 各回答者の満足度が異なる場合（qj_noisy_responder）をランダム設定
p_responders = np.random.uniform(0.35, 0.85, size=n)
qj_noisy_responder = np.zeros(len(categories))
for p in p_responders:
    qj_responder = [np.math.comb(4, j) * (p ** j) * ((1 - p) ** (4 - j)) for j in categories]
    qj_noisy_responder += qj_responder
qj_noisy_responder /= np.sum(qj_noisy_responder)

# シミュレーション開始
for _ in range(n_simulations):
    # 忠実なモデル（均一）
    obs_true = np.random.multinomial(n, qj_true)
    exp_true = n * np.array(qj_true)
    chi2_stat_true = np.sum((obs_true - exp_true) ** 2 / exp_true)  # カイ二乗統計量
    likelihood_stat_true = 2 * np.sum(obs_true * np.log(obs_true / exp_true, where=(obs_true > 0)))  # 尤度比統計量
    chi_squared_true.append(chi2_stat_true)
    likelihood_ratio_true.append(likelihood_stat_true)

    # 各番組の満足度が異なる場合
    obs_noisy_program = np.random.multinomial(n, qj_noisy_program)
    exp_noisy_program = n * np.array(qj_true)
    chi2_stat_noisy_program = np.sum((obs_noisy_program - exp_noisy_program) ** 2 / exp_noisy_program)  # カイ二乗統計量
    likelihood_stat_noisy_program = 2 * np.sum(obs_noisy_program * np.log(obs_noisy_program / exp_noisy_program, where=(obs_noisy_program > 0)))  # 尤度比統計量
    chi_squared_noisy_program.append(chi2_stat_noisy_program)
    likelihood_ratio_noisy_program.append(likelihood_stat_noisy_program)

    # 各回答者の満足度が異なる場合
    obs_noisy_responder = np.random.multinomial(n, qj_noisy_responder)
    exp_noisy_responder = n * np.array(qj_true)
    chi2_stat_noisy_responder = np.sum((obs_noisy_responder - exp_noisy_responder) ** 2 / exp_noisy_responder)  # カイ二乗統計量
    likelihood_stat_noisy_responder = 2 * np.sum(obs_noisy_responder * np.log(obs_noisy_responder / exp_noisy_responder, where=(obs_noisy_responder > 0)))  # 尤度比統計量
    chi_squared_noisy_responder.append(chi2_stat_noisy_responder)
    likelihood_ratio_noisy_responder.append(likelihood_stat_noisy_responder)

# 平均値を計算
mean_chi2_true = np.mean(chi_squared_true)
mean_likelihood_true = np.mean(likelihood_ratio_true)
mean_chi2_noisy_program = np.mean(chi_squared_noisy_program)
mean_likelihood_noisy_program = np.mean(likelihood_ratio_noisy_program)
mean_chi2_noisy_responder = np.mean(chi_squared_noisy_responder)
mean_likelihood_noisy_responder = np.mean(likelihood_ratio_noisy_responder)

# 各統計量に対応するp値を計算
p_value_chi2_true = 1 - chi2.cdf(mean_chi2_true, df=3)
p_value_likelihood_true = 1 - chi2.cdf(mean_likelihood_true, df=3)
p_value_chi2_noisy_program = 1 - chi2.cdf(mean_chi2_noisy_program, df=3)
p_value_likelihood_noisy_program = 1 - chi2.cdf(mean_likelihood_noisy_program, df=3)
p_value_chi2_noisy_responder = 1 - chi2.cdf(mean_chi2_noisy_responder, df=3)
p_value_likelihood_noisy_responder = 1 - chi2.cdf(mean_likelihood_noisy_responder, df=3)

# 結果の出力
print(f"忠実なモデル（均一）の平均カイ二乗統計量: {mean_chi2_true:.4f}, 平均p値: {p_value_chi2_true:.4f}")
print(f"忠実なモデル（均一）の平均尤度比統計量: {mean_likelihood_true:.4f}, 平均p値: {p_value_likelihood_true:.4f}")
print(f"各番組の満足度が異なる場合の平均カイ二乗統計量: {mean_chi2_noisy_program:.4f}, 平均p値: {p_value_chi2_noisy_program:.4f}")
print(f"各番組の満足度が異なる場合の平均尤度比統計量: {mean_likelihood_noisy_program:.4f}, 平均p値: {p_value_likelihood_noisy_program:.4f}")
print(f"各回答者の満足度が異なる場合の平均カイ二乗統計量: {mean_chi2_noisy_responder:.4f}, 平均p値: {p_value_chi2_noisy_responder:.4f}")
print(f"各回答者の満足度が異なる場合の平均尤度比統計量: {mean_likelihood_noisy_responder:.4f}, 平均p値: {p_value_likelihood_noisy_responder:.4f}")

# カイ二乗分布のPDFを計算
x = np.linspace(0, 20, 500)  # x軸の範囲
pdf = chi2.pdf(x, df=3)  # 自由度3のPDF

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label="自由度3のカイ二乗分布", linewidth=2)
plt.axvline(mean_chi2_true, color="blue", linestyle="--", label=f"忠実モデルのカイ二乗: {mean_chi2_true:.2f}")
plt.axvline(mean_likelihood_true, color="cyan", linestyle="--", label=f"忠実モデルの尤度比: {mean_likelihood_true:.2f}")
plt.axvline(mean_chi2_noisy_program, color="green", linestyle="--", label=f"番組ごとに異なるカイ二乗: {mean_chi2_noisy_program:.2f}")
plt.axvline(mean_likelihood_noisy_program, color="lime", linestyle="--", label=f"番組ごとに異なる尤度比: {mean_likelihood_noisy_program:.2f}")
plt.axvline(mean_chi2_noisy_responder, color="red", linestyle="--", label=f"回答者ごとに異なるカイ二乗: {mean_chi2_noisy_responder:.2f}")
plt.axvline(mean_likelihood_noisy_responder, color="magenta", linestyle="--", label=f"回答者ごとに異なる尤度比: {mean_likelihood_noisy_responder:.2f}")
plt.xlabel("$X^2$ 値", fontsize=12)
plt.ylabel
plt.ylabel("確率密度", fontsize=12)
plt.title("自由度3のカイ二乗分布と統計量", fontsize=14)
plt.legend(fontsize=12)  # ここで凡例を表示
plt.grid(alpha=0.5)
plt.show()
忠実なモデル（均一）の平均カイ二乗統計量: 4.1751, 平均p値: 0.2432
忠実なモデル（均一）の平均尤度比統計量: 4.3247, 平均p値: 0.2285
各番組の満足度が異なる場合の平均カイ二乗統計量: 13.2002, 平均p値: 0.0042
各番組の満足度が異なる場合の平均尤度比統計量: 16.4412, 平均p値: 0.0009
各回答者の満足度が異なる場合の平均カイ二乗統計量: 10.3066, 平均p値: 0.0161
各回答者の満足度が異なる場合の平均尤度比統計量: 8.5302, 平均p値: 0.0362

忠実なモデルの平均p-値は大きく、モデルによく適合していることが分かります。一方、各番組または各回答者の満足度が異なるモデルではp-値が0.05以下となり、モデルがデータに適合していないことが分かります。

2014 Q5(3)

投稿日:2024年7月23日 　最終更新日:2025年1月17日

テレビ番組の視聴満足回数の調査で満足と答えた回数のそれぞれ人数の適合度検定のための尤度比検定量と自由度を求めました。

　

コード

多項確率q_jが、q_j = q_j(p) = 4C_j p^j (1 - p)^{4-j}, \quad j = 0, \ldots, 4に忠実に基づく場合と、ランダムにノイズを加えた（ズレた）場合についてシミュレーションを行い、尤度比統計量とp-値を計算して比較します。

# 2014 Q5(3)  2025.1.17

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

# シミュレーションの設定
n = 100  # 人数
p_true = 0.6  # 真の成功確率
categories = np.arange(5)  # カテゴリ (0～4)
n_simulations = 1000  # シミュレーション回数

# 二項分布モデルに基づく q_j の計算 (忠実な設定)
qj_true = [np.math.comb(4, j) * (p_true ** j) * ((1 - p_true) ** (4 - j)) for j in categories]

# 二項分布に従わない (qj_noisy): q_j をランダムに変更してズラす
qj_noisy = np.array(qj_true) + np.random.uniform(-0.1, 0.1, size=len(qj_true))
qj_noisy = np.clip(qj_noisy, 0, None)  # 確率が0未満にならないようにする
qj_noisy = qj_noisy / np.sum(qj_noisy)  # 確率なので正規化

# 尤度比統計量を記録するリスト
likelihood_ratio_true = []
p_values_true = []
likelihood_ratio_noisy = []
p_values_noisy = []

for _ in range(n_simulations):
    # 忠実なモデルの観測データ
    obs_true = np.random.multinomial(n, qj_true)
    exp_true = n * np.array(qj_true)
    likelihood_stat_true = 2 * np.sum(obs_true * np.log(obs_true / exp_true, where=(obs_true > 0)))
    p_val_true = 1 - chi2.cdf(likelihood_stat_true, df=3)  # 自由度 = カテゴリ数 - 1 - 推定パラメータ数
    likelihood_ratio_true.append(likelihood_stat_true)
    p_values_true.append(p_val_true)

    # ズレたモデルの観測データ
    obs_noisy = np.random.multinomial(n, qj_noisy)
    exp_noisy = n * np.array(qj_true)  # 忠実な期待値を使用
    likelihood_stat_noisy = 2 * np.sum(obs_noisy * np.log(obs_noisy / exp_noisy, where=(obs_noisy > 0)))
    p_val_noisy = 1 - chi2.cdf(likelihood_stat_noisy, df=3)
    likelihood_ratio_noisy.append(likelihood_stat_noisy)
    p_values_noisy.append(p_val_noisy)

# 平均尤度比統計量を計算
mean_likelihood_true = np.mean(likelihood_ratio_true)
mean_likelihood_noisy = np.mean(likelihood_ratio_noisy)

# カイ二乗分布のPDFを計算
x = np.linspace(0, 20, 500)  # x軸の範囲
pdf = chi2.pdf(x, df=3)  # 自由度3のPDF

# 結果の出力
print(f"忠実なモデルの平均尤度比統計量: {mean_likelihood_true:.4f}")
print(f"忠実なモデルの平均p値: {np.mean(p_values_true):.4f}")
print(f"ズレたモデルの平均尤度比統計量: {mean_likelihood_noisy:.4f}")
print(f"ズレたモデルの平均p値: {np.mean(p_values_noisy):.4f}")

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label="自由度3のカイ二乗分布", linewidth=2)
plt.axvline(mean_likelihood_true, color="blue", linestyle="--", label=f"忠実モデルの平均: {mean_likelihood_true:.2f}")
plt.axvline(mean_likelihood_noisy, color="red", linestyle="--", label=f"ズレモデルの平均: {mean_likelihood_noisy:.2f}")
plt.xlabel("$X^2$ 値", fontsize=12)
plt.ylabel("確率密度", fontsize=12)
plt.title("自由度3のカイ二乗分布と平均尤度比統計量", fontsize=14)
plt.legend(fontsize=12)
plt.grid(alpha=0.5)
plt.show()
忠実なモデルの平均尤度比統計量: 3.9418
忠実なモデルの平均p値: 0.3909
ズレたモデルの平均尤度比統計量: 16.7360
ズレたモデルの平均p値: 0.0081

忠実なモデルの平均p-値は大きく、モデルによく適合していることが分かります。一方、ランダムにノイズを加えた（ズレた）モデルではp-値が0.05以下となり、モデルがデータに適合していないことが分かります。

2014 Q5(2)

投稿日:2024年7月22日 　最終更新日:2025年1月17日

テレビ番組の視聴満足回数の調査で満足と答えた回数のそれぞれ人数の適合度カイ二乗統計量と自由度を求めました。

　

コード

多項確率q_jが、q_j = q_j(p) = 4C_j p^j (1 - p)^{4-j}, \quad j = 0, \ldots, 4に忠実に基づく場合と、ランダムにノイズを加えた（ズレた）場合についてシミュレーションを行い、適合度カイ二乗統計量とp-値を計算して比較します。

# 2014 Q5(2)  2025.1.16

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

# シミュレーションの設定
n = 100  # 人数
p_true = 0.6  # 真の成功確率
categories = np.arange(5)  # カテゴリ (0～4)
n_simulations = 1000  # シミュレーション回数

# 二項分布モデルに基づく q_j の計算 (忠実な設定)
qj_true = [np.math.comb(4, j) * (p_true ** j) * ((1 - p_true) ** (4 - j)) for j in categories]

# 二項分布に従わない (qj_noisy): q_j をランダムに変更してズラす
qj_noisy = np.array(qj_true) + np.random.uniform(-0.1, 0.1, size=len(qj_true))
qj_noisy = np.clip(qj_noisy, 0, None)  # 確率が0未満にならないようにする
qj_noisy = qj_noisy / np.sum(qj_noisy)  # 確率なので正規化

# カイ二乗統計量と p 値を記録するリスト
chi_squared_true = []
p_values_true = []
chi_squared_noisy = []
p_values_noisy = []

for _ in range(n_simulations):
    # 忠実なモデルの観測データ
    obs_true = np.random.multinomial(n, qj_true)
    exp_true = n * np.array(qj_true)
    chi2_stat_true = np.sum((obs_true - exp_true) ** 2 / exp_true)
    p_val_true = 1 - chi2.cdf(chi2_stat_true, df=3)  # 自由度 = カテゴリ数 - 1 - 推定パラメータ数
    chi_squared_true.append(chi2_stat_true)
    p_values_true.append(p_val_true)

    # ズレたモデルの観測データ
    obs_noisy = np.random.multinomial(n, qj_noisy)
    exp_noisy = n * np.array(qj_true)  # 忠実な期待値を使用
    chi2_stat_noisy = np.sum((obs_noisy - exp_noisy) ** 2 / exp_noisy)
    p_val_noisy = 1 - chi2.cdf(chi2_stat_noisy, df=3)
    chi_squared_noisy.append(chi2_stat_noisy)
    p_values_noisy.append(p_val_noisy)

# 平均カイ二乗統計量を計算
mean_chi2_true = np.mean(chi_squared_true)
mean_chi2_noisy = np.mean(chi_squared_noisy)

# カイ二乗分布のPDFを計算
x = np.linspace(0, 20, 500)  # x軸の範囲
pdf = chi2.pdf(x, df=3)  # 自由度3のPDF

# 結果の出力
print(f"忠実なモデルの平均カイ二乗統計量: {mean_chi2_true:.4f}")
print(f"忠実なモデルの平均p値: {np.mean(p_values_true):.4f}")
print(f"ズレたモデルの平均カイ二乗統計量: {mean_chi2_noisy:.4f}")
print(f"ズレたモデルの平均p値: {np.mean(p_values_noisy):.4f}")

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label="自由度3のカイ二乗分布", linewidth=2)
plt.axvline(mean_chi2_true, color="blue", linestyle="--", label=f"忠実モデルの平均: {mean_chi2_true:.2f}")
plt.axvline(mean_chi2_noisy, color="red", linestyle="--", label=f"ズレモデルの平均: {mean_chi2_noisy:.2f}")
plt.xlabel("$X^2$ 値", fontsize=12)
plt.ylabel("確率密度", fontsize=12)
plt.title("自由度3のカイ二乗分布と平均カイ二乗統計量", fontsize=14)
plt.legend(fontsize=12)
plt.grid(alpha=0.5)
plt.show()
忠実なモデルの平均カイ二乗統計量: 3.8112
忠実なモデルの平均p値: 0.4026
ズレたモデルの平均カイ二乗統計量: 13.4317
ズレたモデルの平均p値: 0.0236

忠実なモデルの平均p-値は大きく、モデルによく適合していることが分かります。一方、ランダムにノイズを加えた（ズレた）モデルではp-値が0.05以下となり、モデルがデータに適合していないことが分かります。

2014 Q5(1)

投稿日:2024年7月20日 　最終更新日:2025年1月15日

テレビ番組の視聴満足回数の調査で満足と答えた回数のそれぞれ人数との組みが多項分布に従うのを利用し多項確率を構成する二項分布の確率pの最尤推定量を求めました。

　

コード

真の確率pに基づいてシミュレーションを行い、最尤推定値\hat{p} = \frac{1}{4n} \sum_{j=0}^4 j n_jを計算し、真の確率pと比較します。

# 2014 Q5(1)  2025.1.15

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションのパラメータ
n_simulations = 10000  # シミュレーション回数
n = 100  # 調査対象の人数
p_true = 0.6  # 真の満足確率

# nj のシミュレーション (j = 0, 1, 2, 3, 4)
data = np.random.binomial(4, p_true, size=n_simulations * n).reshape(n_simulations, n)

# 各シミュレーションでの最尤推定値を計算
p_estimates = np.sum(data, axis=1) / (4 * n)

# 平均値（推定されたp^の期待値）を計算
mean_p_estimate = np.mean(p_estimates)

# 推定されたp^の平均値と真の値を出力
print(f"推定された最尤推定値の平均: {mean_p_estimate:.4f}")
print(f"真の値: {p_true:.4f}")

# グラフの描画
plt.hist(p_estimates, bins=30, density=True, alpha=0.75, label="推定された $\\hat{p}$")
plt.axvline(p_true, color='red', linestyle='dashed', label=f"真の値 $p = {p_true}$")
plt.xlabel("推定された $\\hat{p}$")
plt.ylabel("密度")
plt.title("最尤推定値 $\\hat{p}$ の分布")
plt.legend()
plt.grid(alpha=0.5)
plt.show()
推定された最尤推定値の平均: 0.5997
真の値: 0.6000

シミュレーションの結果、最尤推定値\hat{p} = \frac{1}{4n} \sum_{j=0}^4 j n_jは、真の確率pに近い値を示しました。

2014 Q4(5)

投稿日:2024年7月19日 　最終更新日:2025年1月14日

5個の物体の重さを、1つずつ2回量る場合と、2つずつ10通り量る場合で、物体の重さが全て同じか、そうでないかF検定するとき、どちらの場合が効率的かを示しました。

　

コード

(1)の計画行列X1と、(2)の計画行列X2を、それぞれ0と1を反転させた計画行列X3,X4を加え、それぞれの分散と非心度λを計算して比較してみます。

# 2014 Q4(5)  2025.1.14

import numpy as np

# 設計行列 (1)～(4)
X1 = np.array([
    [1, 0, 0, 0, 0],
    [1, 0, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 1]
])

X2 = np.array([
    [1, 1, 0, 0, 0],
    [1, 0, 1, 0, 0],
    [1, 0, 0, 1, 0],
    [1, 0, 0, 0, 1],
    [0, 1, 1, 0, 0],
    [0, 1, 0, 1, 0],
    [0, 1, 0, 0, 1],
    [0, 0, 1, 1, 0],
    [0, 0, 1, 0, 1],
    [0, 0, 0, 1, 1]
])

X3 = 1 - X2  # (2) の反転
X4 = 1 - X1  # (1) の反転

# 設計行列リスト
design_matrices = [X1, X2, X3, X4]

# 真の値 θ
true_theta = np.array([10, 20, 30, 40, 50])

# 結果を格納するリスト
results = []

for i, X in enumerate(design_matrices, start=1):
    # 分散行列を計算
    variance_matrix = np.linalg.inv(X.T @ X)
    # 平均分散を計算
    mean_variance = np.trace(variance_matrix) / len(variance_matrix)
    
    # 非心度 λ を計算
    theta_mean = np.mean(true_theta)
    ones_vector = np.ones((X.shape[0], 1))  # X の行数に応じた列ベクトル
    theta_centered = true_theta - theta_mean  # θ から平均を引いたベクトル
    numerator = theta_centered.T @ X.T @ X @ theta_centered
    denominator = ones_vector.T @ X @ X.T @ ones_vector
    non_centrality_lambda = numerator.item()  # スカラーに変換
    
    # 結果を保存
    results.append((i, mean_variance, non_centrality_lambda))

# 結果の出力
for i, mean_variance, non_centrality_lambda in results:
    print(f"(X{i}) の平均分散: {mean_variance:.6f}, 非心度 λ: {non_centrality_lambda:.6f}")
(X1) の平均分散: 0.500000, 非心度 λ: 2000.000000
(X2) の平均分散: 0.291667, 非心度 λ: 3000.000000
(X3) の平均分散: 0.277778, 非心度 λ: 3000.000000
(X4) の平均分散: 0.406250, 非心度 λ: 2000.000000

(2)の計画行列X2は、(1)の計画行列X1に比べ、平均分散が小さく、非心度 λが大きいため、より優れた計画行列と言えます。また、X2の0と1を反転したX3は、X2よりも分散がさらに小さく、非心度 λが変化しないことから、この中で最も優れた計画行列と考えられます。

X3 = np.array([
    [0, 0, 1, 1, 1],
    [0, 1, 0, 1, 1],
    [0, 1, 1, 0, 1],
    [0, 1, 1, 1, 0],
    [1, 0, 0, 1, 1],
    [1, 0, 1, 0, 1],
    [1, 0, 1, 1, 0],
    [1, 1, 0, 0, 1],
    [1, 1, 0, 1, 0],
    [1, 1, 1, 0, 0]
])

2014 Q4(4)

投稿日:2024年7月18日 　最終更新日:2025年1月14日

5個の物体の重さを2つずつ10通り量る場合で、誤差の分散を未知とし、全ての物体の重さが等しいとはいえないとする対立仮説のF検定量の非心度を求めました。

　

コード

真の値\theta_1, \dots, \theta_5=[10,10,10,10,10]から\theta_1, \dots, \theta_5=[10,20,30,40,50]へ徐々に変化させて、非心度λを計算し、その変化をグラフにプロットしてみます。測定法(1)と測定法(2)を同時に表示し、分散の違いを確認してみます。

# 2014 Q4(4)  2025.1.13

import numpy as np
import matplotlib.pyplot as plt

# 初期設定
theta_initial = np.array([10, 10, 10, 10, 10])  # 初期の真の値
theta_final = np.array([10, 20, 30, 40, 50])  # 最終的な真の値
steps = 50  # 変化のステップ数

# θを徐々に変化させる
theta_values = [theta_initial + (theta_final - theta_initial) * t / steps for t in range(steps + 1)]

# 非心度 λ を計算
# 測定法(1): 係数は2
lambdas_1 = [2 * np.sum((theta - np.mean(theta)) ** 2) for theta in theta_values]
# 測定法(2): 係数は3
lambdas_2 = [3 * np.sum((theta - np.mean(theta)) ** 2) for theta in theta_values]

# サンプルのインデックス（途中でいくつかの点を表示）
sample_indices = [0, 10, 20, 30, 40, 50]

# 真の値と対応する λ を取得 (測定法(1) と測定法(2))
sample_theta_values = [theta_values[i] for i in sample_indices]
sample_lambdas_1 = [lambdas_1[i] for i in sample_indices]
sample_lambdas_2 = [lambdas_2[i] for i in sample_indices]

# 視覚化
plt.figure(figsize=(12, 8))

# 測定法(1) のプロット
plt.plot(range(steps + 1), lambdas_1, label="測定法(1): λ", color="blue")
plt.scatter(sample_indices, sample_lambdas_1, color="blue", label="測定法(1) サンプル点")
for i, (theta, lam) in enumerate(zip(sample_theta_values, sample_lambdas_1)):
    plt.text(sample_indices[i], lam, f"{np.round(theta, 1)}\nλ={lam:.1f}", 
             fontsize=10, ha="center", color="blue")

# 測定法(2) のプロット
plt.plot(range(steps + 1), lambdas_2, label="測定法(2): λ", color="red")
plt.scatter(sample_indices, sample_lambdas_2, color="red", label="測定法(2) サンプル点")
for i, (theta, lam) in enumerate(zip(sample_theta_values, sample_lambdas_2)):
    plt.text(sample_indices[i], lam, f"{np.round(theta, 1)}\nλ={lam:.1f}", 
             fontsize=10, ha="center", color="red")

# グラフの設定
plt.xlabel("変化のステップ数")
plt.ylabel("非心度 λ")
plt.title("真の値の変化による非心度 λ の変化 (測定法(1) と 測定法(2))")
plt.legend()
plt.grid()
plt.show()

真の値\theta_1, \dots, \theta_5のばらつきが大きくなるにつれて、非心度λも非線形的に増加することが確認できました。また測定法(1)のほうが測定法(2)より分散が大きいことが確認できました。

2014 Q4(4)

投稿日:2024年7月17日 　最終更新日:2025年1月14日

5個の物体の重さを1つずつ2回量る場合で、誤差の分散を未知とし、全ての物体の重さが等しいとはいえないとする対立仮説のF検定量の非心度を求めました。

　

コード

真の値\theta_1, \dots, \theta_5=[10,10,10,10,10]から\theta_1, \dots, \theta_5=[10,20,30,40,50]へ徐々に変化させて、非心度λを計算し、その変化をグラフにプロットしてみます。

# 2014 Q4(4)  2025.1.12

import numpy as np
import matplotlib.pyplot as plt

# 初期設定
theta_initial = np.array([10, 10, 10, 10, 10])  # 初期の真の値
theta_final = np.array([10, 20, 30, 40, 50])  # 最終的な真の値
steps = 50  # 変化のステップ数

# θを徐々に変化させる
theta_values = [theta_initial + (theta_final - theta_initial) * t / steps for t in range(steps + 1)]

# 非心度 λ を計算 (測定法(1): 係数は2)
lambdas_1 = [2 * np.sum((theta - np.mean(theta)) ** 2) for theta in theta_values]

# サンプルのインデックス（途中でいくつかの点を表示）
sample_indices = [0, 10, 20, 30, 40, 50]

# 真の値と対応する λ を取得 (測定法(1))
sample_theta_values = [theta_values[i] for i in sample_indices]
sample_lambdas_1 = [lambdas_1[i] for i in sample_indices]

# 視覚化
plt.figure(figsize=(10, 6))
plt.plot(range(steps + 1), lambdas_1, label="測定法(1): λ", color="blue")
plt.scatter(sample_indices, sample_lambdas_1, color="red", label="サンプル点")
for i, (theta, lam) in enumerate(zip(sample_theta_values, sample_lambdas_1)):
    plt.text(sample_indices[i], lam, f"{np.round(theta, 1)}\nλ={lam:.1f}", 
             fontsize=10, ha="center", color="blue")

# グラフの設定
plt.xlabel("変化のステップ数")
plt.ylabel("非心度 λ")
plt.title("真の値の変化による非心度 λ の変化 (測定法(1))")
plt.legend()
plt.grid()
plt.show()

真の値\theta_1, \dots, \theta_5のばらつきが大きくなるにつれて、非心度λも非線形的に増加することが確認できました。

2014 Q4(3)

投稿日:2024年7月15日 　最終更新日:2025年1月11日

5個の物体の重さを、1つずつ2回量る場合と、2つずつ10通り量る場合で、誤差の分散を未知とし、全ての物体の重さが等しいとする帰無仮説を検定するF検定量を求めました。

　

コード

真の値\theta_1, \dots, \theta_5=10とし、二つの測定方法(1)と(2)においてシミュレーションを行い、F統計量を計算して、帰無仮説H0が棄却されるかを検定を行います。

# 2014 Q4(3)  2025.1.11

import numpy as np
from scipy.stats import f

# パラメータ設定
n = 10  # 観測データ数
p_1 = 5  # 自由度（対立仮説モデルのパラメータ数）
p_0 = 1  # 自由度（帰無仮説モデルのパラメータ数）
nu_1 = p_1 - p_0  # 分子の自由度
nu_2 = n - p_1  # 分母の自由度

# 設計行列
X1 = np.array([  # 測定法 (1)
    [1, 0, 0, 0, 0],
    [1, 0, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 1]
])

X2 = np.array([  # 測定法 (2)
    [1, 1, 0, 0, 0],
    [1, 0, 1, 0, 0],
    [1, 0, 0, 1, 0],
    [1, 0, 0, 0, 1],
    [0, 1, 1, 0, 0],
    [0, 1, 0, 1, 0],
    [0, 1, 0, 0, 1],
    [0, 0, 1, 1, 0],
    [0, 0, 1, 0, 1],
    [0, 0, 0, 1, 1]
])

X0 = np.ones((n, 1))  # 帰無仮説モデルの設計行列

# シミュレーション設定
true_theta = np.array([10, 10, 10, 10, 10])  # 真の値 (帰無仮説下)
sigma_squared = 1  # 誤差分散
n_simulations = 10000  # シミュレーション回数

# 統計量を保存するリスト
f_statistics_1 = []
f_statistics_2 = []

# シミュレーション
for _ in range(n_simulations):
    # 観測データ生成
    epsilon = np.random.normal(0, np.sqrt(sigma_squared), n)
    x1 = X1 @ true_theta + epsilon
    x2 = X2 @ true_theta + epsilon

    # 最小二乗推定量
    theta_hat_1 = np.linalg.inv(X1.T @ X1) @ X1.T @ x1
    theta_hat_2 = np.linalg.inv(X2.T @ X2) @ X2.T @ x2
    theta_hat_0_1 = np.linalg.inv(X0.T @ X0) @ X0.T @ x1
    theta_hat_0_2 = np.linalg.inv(X0.T @ X0) @ X0.T @ x2

    # 分子と分母の計算
    numerator_1 = np.sum((X1 @ theta_hat_1 - X0 @ theta_hat_0_1) ** 2) / nu_1
    denominator_1 = np.sum((x1 - X1 @ theta_hat_1) ** 2) / nu_2
    f_stat_1 = numerator_1 / denominator_1

    numerator_2 = np.sum((X2 @ theta_hat_2 - X0 @ theta_hat_0_2) ** 2) / nu_1
    denominator_2 = np.sum((x2 - X2 @ theta_hat_2) ** 2) / nu_2
    f_stat_2 = numerator_2 / denominator_2

    f_statistics_1.append(f_stat_1)
    f_statistics_2.append(f_stat_2)

# 平均 F 統計量を計算
f_stat_mean_1 = np.mean(f_statistics_1)
f_stat_mean_2 = np.mean(f_statistics_2)

# 理論 F 分布の平均計算
f_theoretical_mean = nu_2 / (nu_2 - 2) if nu_2 > 2 else None

# P値の計算
p_value_1 = 1 - f.cdf(f_stat_mean_1, nu_1, nu_2)
p_value_2 = 1 - f.cdf(f_stat_mean_2, nu_1, nu_2)

# 有意水準
alpha = 0.05
reject_null_1 = p_value_1 < alpha
reject_null_2 = p_value_2 < alpha

# 結果を出力
print("理論 F 分布の平均 (自由度 ν1 = {0}, ν2 = {1}): {2:.4f}".format(nu_1, nu_2, f_theoretical_mean))

print("\n測定法 (1):")
print("  平均 F 統計量: {:.4f}".format(f_stat_mean_1))
print("  P値: {:.4e}".format(p_value_1))
print("  帰無仮説を棄却するか: {}".format("はい" if reject_null_1 else "いいえ"))

print("\n測定法 (2):")
print("  平均 F 統計量: {:.4f}".format(f_stat_mean_2))
print("  P値: {:.4e}".format(p_value_2))
print("  帰無仮説を棄却するか: {}".format("はい" if reject_null_2 else "いいえ"))
理論 F 分布の平均 (自由度 ν1 = 4, ν2 = 5): 1.6667

測定法 (1):
  平均 F 統計量: 1.6255
  P値: 3.0067e-01
  帰無仮説を棄却するか: いいえ

測定法 (2):
  平均 F 統計量: 1.6490
  P値: 2.9569e-01
  帰無仮説を棄却するか: いいえ

二つの測定方法(1)と(2)においてシミュレーションを行い、F検定の結果、どちらの場合も帰無仮説H0が棄却されませんでした。

2014 Q4(2)

投稿日:2024年7月14日 　最終更新日:2025年1月10日

5個の物体の重さを、1つずつ2回量る場合と、2つずつ10通り量る場合で得た最小二乗法による各物体の重さの推定値の分散を求めました。

　

コード

シミュレーションにより、推定量\hat{\theta}_1, \dots, \hat{\theta}_5の分散を求め、理論分散と比較してみます。

# 2014 Q4(2)  2025.1.10

import numpy as np

# シミュレーション設定
np.random.seed(42)  # 再現性のための乱数シード
n_simulations = 10000  # シミュレーション回数
sigma_squared = 1  # 誤差分散
true_theta = np.array([10, 20, 30, 40, 50])  # 真の値 (任意設定)

# (1) 各物体を2回測定する場合
X1 = np.array([
    [1, 0, 0, 0, 0],
    [1, 0, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 1, 0],
    [0, 0, 0, 0, 1],
    [0, 0, 0, 0, 1]
])

# (2) 2つずつ異なる組み合わせで測定する場合
X2 = np.array([
    [1, 1, 0, 0, 0],
    [1, 0, 1, 0, 0],
    [1, 0, 0, 1, 0],
    [1, 0, 0, 0, 1],
    [0, 1, 1, 0, 0],
    [0, 1, 0, 1, 0],
    [0, 1, 0, 0, 1],
    [0, 0, 1, 1, 0],
    [0, 0, 1, 0, 1],
    [0, 0, 0, 1, 1]
])

# 理論分散
theoretical_variance_1 = sigma_squared / 2
theoretical_variance_2 = 7 * sigma_squared / 24

# シミュレーション結果を格納するリスト
variances_1 = []
variances_2 = []

for _ in range(n_simulations):
    # 誤差を生成
    epsilon = np.random.normal(0, np.sqrt(sigma_squared), 10)

    # 測定値を生成
    x1 = X1 @ true_theta + epsilon  # (1) の場合
    x2 = X2 @ true_theta + epsilon  # (2) の場合

    # 最小二乗推定量を計算
    theta_hat_1 = np.linalg.inv(X1.T @ X1) @ X1.T @ x1
    theta_hat_2 = np.linalg.inv(X2.T @ X2) @ X2.T @ x2

    # 推定量の分散を計算
    variances_1.append(np.mean((theta_hat_1 - true_theta)**2))
    variances_2.append(np.mean((theta_hat_2 - true_theta)**2))

# シミュレーション結果の平均
simulated_variance_1 = np.mean(variances_1)
simulated_variance_2 = np.mean(variances_2)

# 結果の出力
print("(1) 理論分散:", theoretical_variance_1)
print("(1) シミュレーション分散:", simulated_variance_1)
print("(2) 理論分散:", theoretical_variance_2)
print("(2) シミュレーション分散:", simulated_variance_2)
(1) 理論分散: 0.5
(1) シミュレーション分散: 0.5003499949049185
(2) 理論分散: 0.2916666666666667
(2) シミュレーション分散: 0.2926572303618685

シミュレーションによる推定量\hat{\theta}_1, \dots, \hat{\theta}_5の分散は、理論分散とよく一致しました。

2014 Q4(1)-2

投稿日:2024年7月13日 　最終更新日:2025年1月9日

5個の物体の重さを、2つずつ10通り量る場合で、最小二乗法による各物体の重さの推定量を求めました。

　

コード

最小二乗推定量\hat{\theta} = (X^T X)^{-1} X^T xに基づいて計算してみます。

# 2014 Q4(1)-2  2025.1.9

import numpy as np

# 測定値
x = np.array([10.1, 12.3, 13.5, 14.6, 9.8, 11.7, 15.3, 12.8, 13.1, 14.2])

# 設計行列 X を構築 (2つずつ異なる組み合わせで10通り測定)
X = np.array([
    [1, 1, 0, 0, 0],
    [1, 0, 1, 0, 0],
    [1, 0, 0, 1, 0],
    [1, 0, 0, 0, 1],
    [0, 1, 1, 0, 0],
    [0, 1, 0, 1, 0],
    [0, 1, 0, 0, 1],
    [0, 0, 1, 1, 0],
    [0, 0, 1, 0, 1],
    [0, 0, 0, 1, 1]
])

# 最小二乗推定量の計算
X_transpose = X.T
theta_hat = np.linalg.inv(X_transpose @ X) @ X_transpose @ x

# 結果を出力
print("最小二乗推定量 (θ̂):", theta_hat)
最小二乗推定量 (θ̂): [6.21666667 5.01666667 5.38333333 6.78333333 8.45      ]

最小二乗推定量\hat{\theta}が計算されました。

次に、手計算で導いた式に基づいて計算してみます。

# 2014 Q4(1)-2  2025.1.9

import numpy as np

# 測定値
x = np.array([10.1, 12.3, 13.5, 14.6, 9.8, 11.7, 15.3, 12.8, 13.1, 14.2])

# θ^ を問題で導いた式に基づいて計算
theta_hat = np.array([
    (3 * (x[0] + x[1] + x[2] + x[3]) - (x[4] + x[5] + x[6] + x[7] + x[8] + x[9])) / 12,
    (3 * (x[0] + x[4] + x[5] + x[6]) - (x[1] + x[2] + x[3] + x[7] + x[8] + x[9])) / 12,
    (3 * (x[1] + x[4] + x[7] + x[8]) - (x[0] + x[2] + x[3] + x[5] + x[6] + x[9])) / 12,
    (3 * (x[2] + x[5] + x[7] + x[9]) - (x[0] + x[1] + x[3] + x[4] + x[6] + x[8])) / 12,
    (3 * (x[3] + x[6] + x[8] + x[9]) - (x[0] + x[1] + x[2] + x[4] + x[5] + x[7])) / 12
])

# 結果を出力
print("最小二乗推定量 (θ̂):", theta_hat)

最小二乗推定量 (θ̂): [6.21666667 5.01666667 5.38333333 6.78333333 8.45      ]

最小二乗推定量\hat{\theta}が計算され、\hat{\theta} = (X^T X)^{-1} X^T xで計算された値と一致しました。

投稿ナビゲーション

← 過去の投稿