計測工学の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()
# 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()
# -*- 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}")
# -*- 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}%")
# -*- 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")
# -*- 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}")
# -*- 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}")
# -*- 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³")
# -*- 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}")
# -*- 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()
# -*- 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()
# -*- 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()
# 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()
# -*- 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()
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}")