計測工学の統計学


計測工学のPythonコードをまとめました。

https://igses.kyushu-u.ac.jp/tanimoto/wp/wp-content/uploads/2022/03/KeisokuGijutsuGairon.pdf


```Python
# measurement_error_simulation.py
# プログラム名: measurement_error_simulation.py
# Simulate measurement errors including systematic, random, and mistake errors / 系統誤差・偶然誤差・過失誤差のシミュレーション

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameter settings ---
true_value = 100.0              # 真の値 / True value
systematic_bias = 2.0           # 系統誤差（バイアス） / Systematic error (bias)
random_sigma = 1.5              # 偶然誤差の標準偏差 / Standard deviation of random error
outlier_prob = 0.02             # 過失誤差（外れ値）発生確率 / Probability of mistake error (outlier)
outlier_magnitude = 10.0        # 過失誤差の大きさ / Magnitude of mistake error
num_samples = 1000              # 測定回数 / Number of measurements

# --- 測定値生成 / Generate measurements ---
measurements = []
for _ in range(num_samples):
    # 偶然誤差を生成 / Generate random error
    random_error = np.random.normal(loc=0.0, scale=random_sigma)
    # 過失誤差（ランダムに発生） / Mistake error (random outlier)
    if np.random.rand() < outlier_prob:
        mistake_error = np.random.choice([-1, 1]) * outlier_magnitude
    else:
        mistake_error = 0.0
    # 測定値 = 真の値 + 系統誤差 + 偶然誤差 + 過失誤差
    measured = true_value + systematic_bias + random_error + mistake_error
    measurements.append(measured)

measurements = np.array(measurements)
errors = measurements - true_value  # 誤差 = 測定値 - 真の値 / Error = Measured - True

# --- 誤差統計量を表示 / Print error statistics ---
print(f"Mean Error: {np.mean(errors):.2f}")       # 誤差の平均 / Mean of errors
print(f"Std of Error: {np.std(errors):.2f}")      # 誤差の標準偏差 / Std of errors
print(f"Max Error: {np.max(errors):.2f}")         # 最大誤差 / Max error
print(f"Min Error: {np.min(errors):.2f}")         # 最小誤差 / Min error

# --- プロット: 測定値の分布 / Plot distribution of measurements ---
plt.figure()
plt.hist(measurements, bins=30, edgecolor='black')
plt.axvline(true_value, color='red', linestyle='dashed', linewidth=2, label='True Value')
plt.xlabel('Measurement Value')  # 英語ラベル / English label
plt.ylabel('Frequency')          # 英語ラベル / English label
plt.title('Measured Values Distribution')  # 英語タイトル / English title
plt.legend()
plt.show()

# --- プロット: 誤差の分布 / Plot distribution of errors ---
plt.figure()
plt.hist(errors, bins=30, edgecolor='black')
plt.xlabel('Error (Measured - True)')  # 英語ラベル / English label
plt.ylabel('Frequency')                # 英語ラベル / English label
plt.title('Measurement Error Distribution')  # 英語タイトル / English title
plt.show()

```

```Python
# measurement_statistics.py
# プログラム名: measurement_statistics.py
# Simulate and analyze measurement statistics: error, deviation, residual
# 測定の統計解析：誤差、偏差、残差をシミュレーションし計算・可視化

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameter settings ---
true_value = 100.0              # 真の値 (True value)
systematic_bias = 2.0           # 系統誤差 (Systematic bias)
random_sigma = 2.0              # 偶然誤差の標準偏差 (Random error σ)
num_samples = 500               # 測定回数 (Number of measurements)

# --- 測定値生成 / Generate measurements ---
# 真の母平均 m = true_value + systematic_bias
population_mean = true_value + systematic_bias

measurements = true_value + systematic_bias + np.random.normal(0, random_sigma, num_samples)
# 測定値配列 / Measurement array

# --- 統計量計算 / Compute statistics ---
sample_mean = np.mean(measurements)                 # 試料平均 M̄ (Sample mean)
errors = measurements - true_value                   # 誤差 Mi - X (Error)
deviations = measurements - population_mean          # 偏差 Mi - m (Deviation)
residuals = measurements - sample_mean               # 残差 Mi - M̄ (Residual)

# 結果を表示 / Print results
print(f"Population Mean (m): {population_mean:.2f}")
print(f"Sample Mean (M̄): {sample_mean:.2f}")
print(f"Mean Error: {np.mean(errors):.2f}")
print(f"Mean Deviation: {np.mean(deviations):.2f}")
print(f"Mean Residual: {np.mean(residuals):.2f}")

# --- プロット: 測定値の分布 / Plot distribution of measurements ---
plt.figure()
plt.hist(measurements, bins=30, edgecolor='black')
plt.axvline(true_value, color='red', linestyle='dashed', linewidth=2, label='True Value')
plt.axvline(population_mean, color='blue', linestyle='dashdot', linewidth=2, label='Population Mean')
plt.axvline(sample_mean, color='green', linestyle='solid', linewidth=2, label='Sample Mean')
plt.xlabel('Measurement Value')   # 英語ラベル / English label
plt.ylabel('Frequency')           # 英語ラベル / English label
plt.title('Measurement Value Distribution')  # 英語タイトル / English title
plt.legend()
plt.show()

# --- プロット: 誤差・偏差・残差の分布 / Plot distributions of Error, Deviation, Residual ---
fig, axes = plt.subplots(3, 1, figsize=(6, 10))
axes[0].hist(errors, bins=30, edgecolor='black')
axes[0].set_xlabel('Error (Mi - X)')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Error Distribution')

axes[1].hist(deviations, bins=30, edgecolor='black')
axes[1].set_xlabel('Deviation (Mi - m)')
axes[1].set_ylabel('Frequency')
axes[1].set_title('Deviation Distribution')

axes[2].hist(residuals, bins=30, edgecolor='black')
axes[2].set_xlabel('Residual (Mi - M̄)')
axes[2].set_ylabel('Frequency')
axes[2].set_title('Residual Distribution')

plt.tight_layout()
plt.show()

```

```Python
# -*- coding: utf-8 -*-
# プログラム名: statistical_measurement_distribution.py
# Program Name: statistical_measurement_distribution.py

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

# -----------------------------
# パラメータ定義 / Define parameters
# -----------------------------
mu = 50          # 母平均 Population mean
sigma = 10       # 母標準偏差 Population std. deviation
sample_size = 30 # 標本サイズ Sample size
num_samples = 1000  # サンプル数 Number of sampling trials

# -----------------------------
# 無作為標本抽出 / Random sampling
# -----------------------------
np.random.seed(0)
samples = np.random.normal(mu, sigma, (num_samples, sample_size))
sample_means = samples.mean(axis=1)  # 試料平均 Sample means

# -----------------------------
# 標本平均の分布可視化 / Visualization of sample mean distribution
# -----------------------------
x = np.linspace(mu - 4*sigma, mu + 4*sigma, 1000)
pdf = norm.pdf(x, loc=mu, scale=sigma)  # 元の正規分布 PDF
pdf_sample_mean = norm.pdf(x, loc=mu, scale=sigma/np.sqrt(sample_size))  # 標本平均分布のPDF

# -----------------------------
# プロット / Plot
# -----------------------------
plt.figure(figsize=(10, 6))

# 元の母集団分布 Original normal distribution
plt.plot(x, pdf, label='Population Distribution $N(\\mu, \\sigma^2)$', linestyle='--')

# 標本平均の分布 Sample mean distribution
plt.plot(x, pdf_sample_mean, label=f'Sample Mean Distribution $N(\\mu, \\sigma^2/n)$ (n={sample_size})')

# 実際の標本平均ヒストグラム Histogram of sample means
plt.hist(sample_means, bins=30, density=True, alpha=0.4, label='Sample Mean Histogram')

# ラベルとタイトル / Labels and title
plt.title("Distribution of Measurement Values and Sample Means")
plt.xlabel("Measured Value")
plt.ylabel("Probability Density")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# -----------------------------
# 測定精度と精密度の確認 / Accuracy and Precision
# -----------------------------
true_value = mu
measured_means = np.mean(sample_means)
bias = measured_means - true_value  # 系統誤差 Bias (Accuracy error)
random_error = np.std(sample_means) # ランダム誤差 Random error (Precision)

# 結果出力 / Print results
print("===== Statistical Summary =====")
print(f"True Population Mean (μ): {mu}")
print(f"Measured Mean of Sample Means: {measured_means:.2f}")
print(f"Bias (Accuracy error): {bias:.2f}")
print(f"Standard Deviation of Sample Means (Precision): {random_error:.2f}")

```



```Python
# -*- coding: utf-8 -*-
# プログラム名: normal_distribution_error_precision_accuracy.py
# Program Name: normal_distribution_error_precision_accuracy.py

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

# -----------------------------
# 母平均・母標準偏差 / True parameters
# -----------------------------
mu_true = 0       # 真の値 (True value)
sigma_true = 1    # 母標準偏差 (True std)

# -----------------------------
# 測定データシミュレーション / Simulate measurements
# -----------------------------
np.random.seed(42)
num_measurements = 1000

# シナリオ (a) 正確で精密 / Accurate & Precise
data_a = np.random.normal(mu_true, 0.2, num_measurements)

# シナリオ (b) 精密だが正確でない / Precise but not accurate
data_b = np.random.normal(mu_true + 1.0, 0.2, num_measurements)

# シナリオ (c) 正確だが精密でない / Accurate but not precise
data_c = np.random.normal(mu_true, 1.0, num_measurements)

# シナリオ (d) 精度も精密度も悪い / Neither accurate nor precise
data_d = np.random.normal(mu_true + 1.0, 1.0, num_measurements)

# -----------------------------
# 描画 / Plot
# -----------------------------
x = np.linspace(-3, 4, 1000)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
scenarios = [data_a, data_b, data_c, data_d]
titles = ['(a) Accurate and Precise',
          '(b) Precise but not Accurate',
          '(c) Accurate but not Precise',
          '(d) Neither Accurate nor Precise']

for ax, data, title in zip(axes.ravel(), scenarios, titles):
    ax.hist(data, bins=50, density=True, alpha=0.6, label='Measured')
    mean = np.mean(data)
    std = np.std(data)
    ax.plot(x, norm.pdf(x, mean, std), label=f'N({mean:.2f}, {std:.2f}²)', linestyle='--')
    ax.axvline(mu_true, color='red', linestyle=':', label='True Value (μ)')
    ax.set_title(title)
    ax.legend()
    ax.grid(True)

plt.suptitle("Measurement Accuracy and Precision Scenarios")
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()

# -----------------------------
# 標準正規分布と確率計算 / Standard normal distribution
# -----------------------------
z_lower = -2
z_upper = 2
prob_within_2sigma = norm.cdf(z_upper) - norm.cdf(z_lower)

print("===== Standard Normal Distribution =====")
print(f"Probability within ±2σ (Z between -2 and 2): {prob_within_2sigma * 100:.2f}%")

```

```Python
# -*- coding: utf-8 -*-
# プログラム名: significant_figures_error_propagation_statistics.py
# Program Name: significant_figures_error_propagation_statistics.py

import numpy as np

# -----------------------------
# 有効数字の表示関数 / Format with significant figures
# -----------------------------
def format_sigfig(value, sigfigs):
    from decimal import Decimal
    return f"{Decimal(value):.{sigfigs}g}"

# -----------------------------
# 加減乗除の誤差伝播 / Error propagation for arithmetic
# -----------------------------
def add_sub_error(val1, err1, val2, err2):
    result = val1 + val2
    total_error = err1 + err2
    return result, total_error

def mul_div_error(val1, err1, val2, err2, operation='*'):
    # 相対誤差を足し合わせる（乗除共通）
    rel_err1 = err1 / val1
    rel_err2 = err2 / val2
    if operation == '*':
        result = val1 * val2
    elif operation == '/':
        result = val1 / val2
    total_rel_error = rel_err1 + rel_err2
    total_error = result * total_rel_error
    return result, total_error

# -----------------------------
# 算術平均と標準偏差 / Arithmetic mean and sample std
# -----------------------------
def compute_statistics(data):
    mean = np.mean(data)
    std = np.std(data, ddof=1)  # 標本標準偏差 (sample standard deviation)
    return mean, std

# -----------------------------
# 【演習1】有効数字を含む球の体積計算
# -----------------------------
def calc_sphere_volume(diameter_mm):
    radius = diameter_mm / 2 / 1000  # convert mm to m
    volume = (4/3) * np.pi * radius**3
    return volume  # [m^3]

d1 = 20.05  # mm, 有効数字4桁
d2 = 20.1   # mm, 有効数字3桁

v1 = calc_sphere_volume(d1)
v2 = calc_sphere_volume(d2)

print("===== Exercise 1: Sphere Volume =====")
print(f"Diameter = {d1} mm → Volume = {format_sigfig(v1, 4)} m³ (4 sigfigs)")
print(f"Diameter = {d2} mm → Volume = {format_sigfig(v2, 3)} m³ (3 sigfigs)")

# -----------------------------
# 【演習2】20回測定値から平均と標準偏差
# -----------------------------
measurements = [
    11.45, 11.15, 11.17, 12.10, 12.05, 12.01, 12.02, 12.19, 11.84, 12.54,
    12.21, 11.99, 12.11, 12.01, 11.85, 12.00, 11.90, 12.08, 12.02, 12.03
]

mean_val, std_val = compute_statistics(measurements)
print("\n===== Exercise 2: Sample Mean and Std Dev =====")
print(f"Sample Mean (M̄) = {mean_val:.3f}")
print(f"Sample Standard Deviation (s) = {std_val:.3f}")

# -----------------------------
# 【拡張】誤差伝播の例
# -----------------------------
m1, e1 = 1.23, 0.01  # 質量 [kg], 誤差 [kg]
m2, e2 = 0.456, 0.001
m3, e3 = 7.8, 0.1

sum_mass, err_mass = add_sub_error(m1, e1, m2, e2)
sum_mass, err_mass = add_sub_error(sum_mass, err_mass, m3, e3)

print("\n===== Error Propagation (Addition) =====")
print(f"Sum of masses = {sum_mass:.3f} ± {err_mass:.3f} kg")

```



```Python
# -*- coding: utf-8 -*-
# プログラム名: error_propagation_examples.py
# Program Name: error_propagation_examples.py

import numpy as np

# -----------------------------
# ガウスの誤差伝播公式 / Gaussian error propagation
# -----------------------------
def propagate_error(q_func, variables, errors, partials):
    """
    q_func: 計算されたqの値 (float)
    variables: 各変数の値のリスト
    errors: 各変数の誤差のリスト
    partials: 各変数に関する偏微分係数のリスト (∂q/∂xi)
    """
    squared_terms = [(df_dx * err)**2 for df_dx, err in zip(partials, errors)]
    total_error = np.sqrt(np.sum(squared_terms))
    return q_func, total_error

# -----------------------------
# 演習4: 円柱体積 V = π * (d/2)^2 * h
# -----------------------------
print("===== Exercise 4: Volume of Cylinder =====")
d = 14.38  # mm
delta_d = 0.43
h = 6.42   # mm
delta_h = 0.57

# 体積の計算 / Volume in mm^3
V = np.pi * (d / 2)**2 * h

# 偏微分 ∂V/∂d = π * h * d / 2, ∂V/∂h = π * (d/2)^2
dV_dd = np.pi * h * d / 2
dV_dh = np.pi * (d / 2)**2

_, delta_V = propagate_error(V, [d, h], [delta_d, delta_h], [dV_dd, dV_dh])

print(f"Volume V = {V:.1f} ± {delta_V:.1f} mm³")

# -----------------------------
# 演習5: 時間 t = L / v
# -----------------------------
print("\n===== Exercise 5: Travel Time =====")
v = 54.2        # m/s
delta_v = 0.4
L_km = 1.05     # km
delta_L_km = 0.03
L = L_km * 1000       # convert to meters
delta_L = delta_L_km * 1000

# 時間計算 / Time
t = L / v

# 偏微分 ∂t/∂L = 1/v, ∂t/∂v = -L/v^2
dt_dL = 1 / v
dt_dv = -L / v**2

_, delta_t = propagate_error(t, [L, v], [delta_L, delta_v], [dt_dL, dt_dv])

print(f"Time t = {t:.2f} ± {delta_t:.2f} s")

# -----------------------------
# 演習3: 平均値の標準偏差
# -----------------------------
print("\n===== Exercise 3: Std. Dev. of Mean =====")
sigma = 0.5
n = 20
sigma_mean = sigma / np.sqrt(n)
print(f"If σ = {sigma}, then σ̄ (standard error of mean) = {sigma_mean:.3f}")

```

```Python
# -*- coding: utf-8 -*-
# プログラム名: least_squares_regression_analysis.py
# Program Name: least_squares_regression_analysis.py

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# -----------------------------
# 仮の実験データセット / Sample experimental data (x, y)
# -----------------------------
# 実験式: y = a * x^2 + b という形を想定し、log変換で線形化なども可能
x = np.array([1, 2, 3, 4, 5])
y = np.array([2.3, 4.1, 9.3, 15.5, 25.1])  # 曲線的な関係（例: y = x^2 + ノイズ）

# -----------------------------
# 線形回帰モデル: y = b0 + b1 * x（単回帰）/ Linear regression
# -----------------------------
X = x.reshape(-1, 1)  # 説明変数を列ベクトルに変換
model = LinearRegression()
model.fit(X, y)
b0 = model.intercept_
b1 = model.coef_[0]

# -----------------------------
# 回帰結果の表示 / Print regression result
# -----------------------------
print("===== Least Squares Linear Regression =====")
print(f"Regression Equation: y = {b0:.3f} + {b1:.3f} * x")
print(f"R² score: {model.score(X, y):.4f}")

# -----------------------------
# プロット / Plot
# -----------------------------
x_range = np.linspace(min(x), max(x), 100)
y_pred = model.predict(x_range.reshape(-1, 1))

plt.figure(figsize=(8, 5))
plt.scatter(x, y, color='blue', label='Measured Data')
plt.plot(x_range, y_pred, color='red', label='Least Squares Fit')
plt.title("Linear Regression using Least Squares Method")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# -----------------------------
# 拡張：非線形関係の線形化（実験式の適用例）
# 例：Y = A * X^B → ln(Y) = ln(A) + B ln(X)
# -----------------------------
print("\n===== Power Law Fit (log-log linearization) =====")
log_x = np.log(x)
log_y = np.log(y)
log_model = LinearRegression()
log_model.fit(log_x.reshape(-1, 1), log_y)

lnA = log_model.intercept_
B = log_model.coef_[0]
A = np.exp(lnA)

print(f"Power law model: y = {A:.3f} * x^{B:.3f}")
print(f"R² (log-log): {log_model.score(log_x.reshape(-1, 1), log_y):.4f}")

```



```Python
# -*- coding: utf-8 -*-
# プログラム名: measurement_uncertainty_analysis.py
# Program Name: measurement_uncertainty_analysis.py

import numpy as np
from scipy.stats import t

# -----------------------------
# 1. 単純な標本平均の95%信頼区間 / CI for average with known std
# -----------------------------
def ci_known_sigma(xbar, sigma, n, conf=0.95):
    t_val = t.ppf((1 + conf) / 2, df=n-1)
    U = t_val * sigma / np.sqrt(n)
    return xbar, U

# -----------------------------
# 2. 不偏分散による信頼区間 / CI from sample
# -----------------------------
def ci_sample_data(data, conf=0.95):
    n = len(data)
    xbar = np.mean(data)
    s = np.std(data, ddof=1)
    t_val = t.ppf((1 + conf) / 2, df=n-1)
    U = t_val * s / np.sqrt(n)
    return xbar, s, U

# -----------------------------
# 3. 間接測定の不確かさ伝播（ガウス） / Gaussian propagation
# -----------------------------
def indirect_uncertainty(value_func, values, errors, partials):
    squared_terms = [(df_dx * err)**2 for df_dx, err in zip(partials, errors)]
    uc = np.sqrt(np.sum(squared_terms))  # 合成標準不確かさ
    return value_func(*values), uc

# -----------------------------
# 4. 感度係数の計算 / Sensitivity coefficients
# -----------------------------
def numerical_sensitivity(f, vals, i, delta=1e-6):
    v_plus = vals.copy()
    v_minus = vals.copy()
    v_plus[i] += delta
    v_minus[i] -= delta
    return (f(*v_plus) - f(*v_minus)) / (2 * delta)

# -----------------------------
# 5. 使用例
# -----------------------------

# [例1] 標準偏差既知の場合の平均値の95%信頼区間
print("===== Example 1: Known σ confidence interval =====")
mean_val, U1 = ci_known_sigma(xbar=24.985, sigma=1.2, n=5)
print(f"Mean = {mean_val:.3f} ± {U1:.3f} mm (95% CI)")

# [例2] 実測データによる平均と標準偏差からの95%信頼区間
print("\n===== Example 2: Sample-based CI =====")
lengths = np.array([26.04, 24.27, 26.30, 26.09, 26.31, 30.28, 30.91, 29.31, 30.28, 30.69])
mean2, std2, U2 = ci_sample_data(lengths)
print(f"Mean = {mean2:.3f} ± {U2:.3f} mm (95% CI, sample-based)")

# [例3] 円柱体積 V = π(d/2)^2 * h の不確かさ
print("\n===== Example 3: Indirect uncertainty (cylinder volume) =====")
d = 14.38  # mm
h = 6.42   # mm
sd_d = 0.43
sd_h = 0.57

def volume(d, h):
    return np.pi * (d / 2)**2 * h

# 数値微分による感度係数
partial_d = numerical_sensitivity(volume, [d, h], i=0)
partial_h = numerical_sensitivity(volume, [d, h], i=1)

V, uc = indirect_uncertainty(volume, [d, h], [sd_d, sd_h], [partial_d, partial_h])
print(f"Volume = {V:.1f} ± {uc:.1f} mm³")

# -----------------------------
# 6. 包括不確かさ（expanded uncertainty）
# -----------------------------
print("\n===== Example 4: Expanded Uncertainty =====")
coverage_factor = 2  # 95%信頼水準
U_exp = coverage_factor * uc
print(f"Expanded uncertainty (U95) = {U_exp:.1f} mm³")

```

```Python
# -*- coding: utf-8 -*-
# プログラム名: final_uncertainty_analysis_cylinder.py
# Program Name: final_uncertainty_analysis_cylinder.py

import numpy as np
from scipy.stats import t

# -----------------------------
# 1. 円柱の体積 V = (π/4)*d^2*h
# -----------------------------
def cylinder_volume(d, h):
    return (np.pi / 4) * d**2 * h

# 偏微分項（感度係数） Sensitivity coefficients
def partials_cylinder(d, h):
    dV_dd = (np.pi / 2) * d * h
    dV_dh = (np.pi / 4) * d**2
    return dV_dd, dV_dh

# -----------------------------
# 2. 不確かさ合成関数（ガウス）
# -----------------------------
def gaussian_uncertainty(V, dV_dd, dV_dh, sd_d, sd_h):
    u2 = (dV_dd * sd_d)**2 + (dV_dh * sd_h)**2
    return np.sqrt(u2)

# -----------------------------
# 3. 精度・精密度・不確かさ評価（スライド式）
# -----------------------------
def total_uncertainty(B, S, coverage='95'):
    t_val = 2 if coverage == '95' else 3
    return np.sqrt(B**2 + (t_val * S)**2)

# -----------------------------
# 4. 修正トンプソンτ法（異常値検出）
# -----------------------------
def modified_thompson_tau(data, alpha=0.05):
    N = len(data)
    mean = np.mean(data)
    std = np.std(data, ddof=1)
    # τ表（95%）から抜粋：最大25件まで対応
    tau_table = {
        3: 1.15, 4: 1.48, 5: 1.71, 6: 1.89, 7: 2.02, 8: 2.13,
        9: 2.21, 10: 2.29, 11: 2.34, 12: 2.41, 13: 2.46, 14: 2.51,
        15: 2.55, 16: 2.59, 17: 2.62, 18: 2.65, 19: 2.68, 20: 2.71,
        25: 2.82
    }
    tau = tau_table.get(N, 2.5)
    threshold = tau * std / np.sqrt(N)
    outliers = [i for i, x in enumerate(data) if abs(x - mean) > threshold]
    return outliers, threshold

# -----------------------------
# 実行ブロック
# -----------------------------

# 入力値（スライドの例）：h=10.05 cm, d=2.03 cm, σh=0.07 cm, σd=0.05 cm
h = 10.05
d = 2.03
sh = 0.07
sd = 0.05

# 絶対精度の B 値（スライド例より）
Bh = 0.05
Bd = 0.005

# 計算
V = cylinder_volume(d, h)
dV_dd, dV_dh = partials_cylinder(d, h)
S = gaussian_uncertainty(V, dV_dd, dV_dh, sd, sh)

# 絶対精度と合成不確かさ
B_total = np.sqrt((dV_dd * Bd)**2 + (dV_dh * Bh)**2)
U_95 = total_uncertainty(B_total, S, coverage='95')

# 結果表示
print("===== Final Uncertainty Evaluation (Cylinder Volume) =====")
print(f"Volume V = {V:.3f} cm³")
print(f"Combined Standard Uncertainty (S) = {S:.3f} cm³")
print(f"Total Accuracy (B) = {B_total:.3f} cm³")
print(f"Expanded Uncertainty (95%): U95 = {U_95:.3f} cm³")

# -----------------------------
# 修正トンプソンτ法の例
# -----------------------------
print("\n===== Outlier Detection (Modified Thompson τ method) =====")
lengths = [10.02, 9.99, 10.01, 10.00, 10.05, 10.09, 10.01, 9.97, 10.03, 10.12]
outliers, tau_threshold = modified_thompson_tau(lengths)
print(f"Threshold = ±{tau_threshold:.4f}")
print(f"Outlier indices: {outliers}")

```


```Python
# -*- coding: utf-8 -*-
# Program Name: simulate_drift_noise_transmission.py

import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# 時間軸 / Time axis
# -----------------------------
t = np.linspace(0, 10, 1000)  # 0〜10秒を1000分割

# -----------------------------
# 信号成分 / Signal components
# -----------------------------
true_signal = np.ones_like(t) * 1.0  # 入力信号（定数値）
drift = 0.05 * t                     # 緩やかなドリフト
noise = np.random.normal(0, 0.03, size=t.shape)  # 雑音（ノイズ）
output_signal = true_signal + drift + noise

# -----------------------------
# 描画 / Plot
# -----------------------------
plt.figure(figsize=(10, 4))
plt.plot(t, output_signal, label='Output signal (with drift and noise)', color='gray')
plt.plot(t, true_signal, label='True input signal', color='black', linestyle='--')
plt.plot(t, true_signal + drift, label='Drift component', color='red', linestyle=':')
plt.xlabel("Time")
plt.ylabel("Signal")
plt.title("Drift and Noise Simulation")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

```



```Python
# -*- coding: utf-8 -*-
# Program Name: sampling_theorem_simulation.py

import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# 元信号（アナログ信号）/ Original analog signal
# -----------------------------
t_cont = np.linspace(0, 1, 1000)
f = 10  # [Hz]
x_cont = np.sin(2 * np.pi * f * t_cont)

# -----------------------------
# サンプリング条件
# -----------------------------
fs_over = 100  # オーバーサンプリング
fs_nyquist = 2 * f  # 理想的サンプリング
fs_under = 8   # アンダーサンプリング（< 2f）

def sample_and_plot(fs, title, subplot_index):
    t_sample = np.arange(0, 1, 1/fs)
    x_sample = np.sin(2 * np.pi * f * t_sample)
    plt.subplot(3, 1, subplot_index)
    plt.plot(t_cont, x_cont, label='Original signal', color='lightgray')
    plt.stem(t_sample, x_sample, linefmt='r-', markerfmt='ro', basefmt=' ', label=f'Sampled (fs={fs} Hz)')
    plt.title(title)
    plt.xlabel("Time [s]")
    plt.ylabel("Amplitude")
    plt.grid(True)
    plt.legend()

# -----------------------------
# 描画 / Plot
# -----------------------------
plt.figure(figsize=(10, 8))
sample_and_plot(fs_over, "Oversampling (fs = 100 Hz)", 1)
sample_and_plot(fs_nyquist, "Nyquist Sampling (fs = 20 Hz)", 2)
sample_and_plot(fs_under, "Undersampling (fs = 8 Hz)", 3)
plt.tight_layout()
plt.show()


```

```Python
# -*- coding: utf-8 -*-
# Program Name: autocorrelation_spectrum_analysis.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate
from scipy.fftpack import fft, fftfreq

# -----------------------------
# 信号の生成（周期 + 雑音）
# -----------------------------
fs = 100  # サンプリング周波数 [Hz]
t = np.linspace(0, 2, 2 * fs, endpoint=False)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.random.randn(len(t))

# -----------------------------
# 自己相関関数（biased方式）/ Autocorrelation
# -----------------------------
autocorr = correlate(signal, signal, mode='full') / len(signal)
lags = np.arange(-len(signal)+1, len(signal))

# -----------------------------
# パワースペクトル密度（PSD） via FFT of autocorr
# Wiener-Khinchinの定理
# -----------------------------
S = np.abs(fft(autocorr))  # 自己相関関数のFFT
freqs = fftfreq(len(S), 1/fs)

# -----------------------------
# プロット / Plot
# -----------------------------
plt.figure(figsize=(12, 6))

# 自己相関関数
plt.subplot(2, 1, 1)
plt.plot(lags/fs, autocorr)
plt.title("Autocorrelation Function")
plt.xlabel("Lag [s]")
plt.ylabel("Correlation")
plt.grid(True)

# パワースペクトル
plt.subplot(2, 1, 2)
plt.plot(freqs[:len(freqs)//2], S[:len(S)//2])  # 正の周波数のみ
plt.title("Power Spectral Density (from Autocorrelation)")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude")
plt.grid(True)

plt.tight_layout()
plt.show()


```



```Python
# Program Name: exponential_decay_signal_analysis.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate
from scipy.fftpack import fft, fftfreq

# 信号生成：指数減衰雑音信号
fs = 1000
t = np.linspace(0, 1, fs)
alpha = 5
signal = np.exp(-alpha * t) * np.random.randn(len(t))

# 自己相関
auto_corr = correlate(signal, signal, mode='full') / len(signal)
lags = np.arange(-len(signal)+1, len(signal))

# スペクトル
spec = np.abs(fft(auto_corr))
freq = fftfreq(len(spec), d=1/fs)

# プロット
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(lags/fs, auto_corr)
plt.title("Autocorrelation (Exponential Decay)")
plt.xlabel("Lag [s]")
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(freq[:len(freq)//2], spec[:len(freq)//2])
plt.title("Power Spectrum")
plt.xlabel("Frequency [Hz]")
plt.grid(True)
plt.tight_layout()
plt.show()

```

```Python
# -*- coding: utf-8 -*-
# プログラム名: prior_uncertainty_evaluation.py
# Program Name: prior_uncertainty_evaluation.py
#
# 事前分析による不確かさ評価をシミュレート
# Simulate a priori uncertainty evaluation for products sampled from a manufacturing line

import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# 1️ 母集団モデル / Population model
# -----------------------------
# 2つの製造条件 (machine A, machine B) を混合正規分布で表現
np.random.seed(0)
N_population = 100000
mu_A, sigma_A = 50.0, 1.0   # Machine A
mu_B, sigma_B = 49.5, 1.5   # Machine B
ratio_A = 0.6               # 60% from machine A

pop_A = np.random.normal(mu_A, sigma_A, int(N_population * ratio_A))
pop_B = np.random.normal(mu_B, sigma_B, int(N_population * (1 - ratio_A)))
population = np.concatenate([pop_A, pop_B])

# 母平均 μ, 母分散 σ²（真値） / True population parameters
mu_true  = np.mean(population)
sigma2_true = np.var(population, ddof=0)

# -----------------------------
# 2️ 事前測定 (サンプリング) / Preliminary sampling
# -----------------------------
n = 20                                    # 測定回数 / sample size
sample = np.random.choice(population, n)  # 無作為抽出 / random sampling
x_bar = np.mean(sample)                   # 標本平均 / sample mean
s2 = np.var(sample, ddof=1)               # 標本分散 / sample variance

# 推定された母平均・母分散 / Estimated population parameters
mu_hat  = x_bar
sigma2_hat = s2

# 標準不確かさ (標本平均) / Standard uncertainty of mean
u_xbar = np.sqrt(sigma2_hat / n)

# -----------------------------
# 3️ 結果表示 / Print summary
# -----------------------------
print("===== Prior Uncertainty Evaluation =====")
print(f"True population mean μ            : {mu_true:.3f}")
print(f"True population variance σ²        : {sigma2_true:.3f}")
print(f"Sample mean x̄ (estimate μ̂)         : {x_bar:.3f}")
print(f"Sample variance s² (estimate σ̂²)    : {s2:.3f}")
print(f"Standard uncertainty u(x̄)          : {u_xbar:.4f}")

# -----------------------------
# 4️ 可視化 / Visualization
# -----------------------------
plt.figure(figsize=(9, 4))

# 母集団ヒストグラム / Population histogram
plt.hist(population, bins=60, alpha=0.3, density=True, label='Population')

# 標本値 / Sample points
plt.scatter(sample, np.zeros_like(sample), color='red', marker='x', label='Sample data')

# 標本平均ライン / Sample mean line
plt.axvline(x_bar, color='red', linestyle='--', label=f'Sample mean x̄ = {x_bar:.2f}')

plt.title("Population vs Sample (Prior Uncertainty Evaluation)")
plt.xlabel("Measurement Value")
plt.ylabel("Probability Density")
plt.legend()
plt.tight_layout()
plt.show()

```

```Python
import numpy as np
from scipy.stats import t
# uncertainty_evaluation_combined.py


# --- Type A標準不確かさ (標本標準偏差 / √n) ---
# --- Type A standard uncertainty (sample std / √n) ---
def typeA_uncertainty(data):
    n = len(data)
    std = np.std(data, ddof=1)
    return std / np.sqrt(n), std

# --- Type B標準不確かさ（例: 一様分布）---
# --- Type B standard uncertainty (e.g. uniform distribution: a/√3) ---
def rectangular_distribution(a):
    return a / np.sqrt(3)

# --- 合成標準不確かさ ---
# --- Combined standard uncertainty ---
def combine_uncertainty(u_list):
    return np.sqrt(sum(u**2 for u in u_list))

# --- 拡張不確かさ（信頼度95%）---
# --- Expanded uncertainty (95% confidence level) ---
def expanded_uncertainty(uc, dof, confidence_level=0.95):
    t_factor = t.ppf((1 + confidence_level) / 2.0, dof)
    return t_factor * uc

# --- 実験データ / Sample measurement data ---
data = np.array([
    11.45, 11.15, 11.17, 12.10, 12.05,
    12.01, 12.02, 12.32, 12.19, 11.84
])

# Type Aの不確かさと標本標準偏差 / Type A uncertainty and std
uA, s = typeA_uncertainty(data)

# Type Bの不確かさ（±0.1仮定） / Type B uncertainty (assume ±0.1)
uB = rectangular_distribution(0.1)

# 合成不確かさ / Combined uncertainty
uc = combine_uncertainty([uA, uB])

# 自由度 / Degrees of freedom
dof = len(data) - 1

# 拡張不確かさ / Expanded uncertainty (95% level)
U95 = expanded_uncertainty(uc, dof)

# 結果出力 / Print results
print("===== 不確かさ評価 / Uncertainty Evaluation =====")
print(f"測定値の平均値 / Sample Mean           : {np.mean(data):.3f}")
print(f"標本標準偏差 / Sample Std Dev           : {s:.3f}")
print(f"Type A 不確かさ / Type A Uncertainty     : {uA:.5f}")
print(f"Type B 不確かさ / Type B Uncertainty     : {uB:.5f}")
print(f"合成不確かさ / Combined Uncertainty uc   : {uc:.5f}")
print(f"拡張不確かさ (95%) / Expanded Uncertainty: {U95:.5f}")
print(f"信頼区間 (95%) / Confidence Interval     : {np.mean(data):.3f} ± {U95:.3f}")

```