実験の統計学



```
import math

# Constants
k = 2  # Coverage factor for expanded uncertainty (k=2 for ~95% confidence level)

# Direct Measurement

# X(Y): Measured value (change this to your value)
X_Y = 5.0  # Example value

# Experimental standard deviation (change this to your value)
std_dev = 0.1  # Example value

# Instrument resolution (change this to your value)
resolution = 0.01  # Example value

# Standard uncertainty due to resolution
u_resolution = resolution / math.sqrt(12)

# Other factors' uncertainty (change this to your value)
u_other = 0.02  # Example value

# Combined standard uncertainty
u_combined = math.sqrt(std_dev**2 + u_resolution**2 + u_other**2)

# Expanded uncertainty
u_expanded = k * u_combined

print("Direct Measurement:")
print(f"Measured Value (X(Y)): {X_Y}")
print(f"Standard Deviation: {std_dev}")
print(f"Standard Uncertainty (Resolution): {u_resolution}")
print(f"Other Factors' Uncertainty: {u_other}")
print(f"Combined Standard Uncertainty: {u_combined}")
print(f"Expanded Uncertainty: {u_expanded}")

# Indirect Measurement

# Function relating measured quantities (example: Z = X*Y, change this to your relation)
def indirect_measurement(x, y):
    return x * y

# Measured values for X and Y (change these to your values)
X = 2.0  # Example value
Y = 3.0  # Example value

# Standard uncertainties for X and Y (change these to your values)
u_X = 0.05  # Example value
u_Y = 0.07  # Example value

# Partial derivatives (change to match your function)
partial_X = Y
partial_Y = X

# Combined standard uncertainty for indirect measurement (propagation of uncertainty)
u_combined_indirect = math.sqrt((partial_X * u_X)**2 + (partial_Y * u_Y)**2)

# Expanded uncertainty for indirect measurement
u_expanded_indirect = k * u_combined_indirect

print("\nIndirect Measurement:")
print(f"Measured Value (Z): {indirect_measurement(X, Y)}")
print(f"Combined Standard Uncertainty (Indirect): {u_combined_indirect}")
print(f"Expanded Uncertainty (Indirect): {u_expanded_indirect}")

```

```

import numpy as np
import matplotlib.pyplot as plt

# Data: Measured values and associated uncertainties
# Example data for 5 measurements
measurements = np.array([5.2, 5.5, 5.3, 5.4, 5.6])  # Measured values
uncertainties = np.array([0.1, 0.15, 0.12, 0.1, 0.14])  # Uncertainties

# Weighted Average Method
weights = 1 / uncertainties**2
weighted_avg = np.sum(weights * measurements) / np.sum(weights)
combined_uncertainty_weighted = np.sqrt(1 / np.sum(weights))

# Equal Interval Averaging Method
equal_interval_avg = np.mean(measurements)
combined_uncertainty_equal = np.sqrt(np.sum((uncertainties**2)) / len(measurements))

# Least Squares Method (Linear Fit: y = ax + b)
x_values = np.arange(len(measurements))  # Example x-values (e.g., time intervals)
coefficients, covariance_matrix = np.polyfit(x_values, measurements, 1, cov=True)
slope, intercept = coefficients
slope_uncertainty, intercept_uncertainty = np.sqrt(np.diag(covariance_matrix))

# Sum of Squared Residuals
residuals = measurements - (slope * x_values + intercept)
sum_squared_residuals = np.sum(residuals**2)

# Best Estimate from Linear Fit
best_estimate = slope * x_values + intercept

# Plotting the results
plt.figure(figsize=(10, 6))

# Plot measured values with uncertainties
plt.errorbar(x_values, measurements, yerr=uncertainties, fmt='o', label='Measured Values')

# Plot weighted average
plt.axhline(weighted_avg, color='red', linestyle='--', label=f'Weighted Average: {weighted_avg:.2f}')

# Plot equal interval average
plt.axhline(equal_interval_avg, color='green', linestyle=':', label=f'Equal Interval Avg: {equal_interval_avg:.2f}')

# Plot least squares fit line
plt.plot(x_values, best_estimate, label=f'Linear Fit: y = {slope:.2f}x + {intercept:.2f}', color='blue')

# Labels and legend
plt.xlabel('Measurement Index')
plt.ylabel('Measured Value')
plt.title('Analysis of Measurements and Uncertainties')
plt.legend()
plt.grid(True)
plt.show()

# Output results
print("Results:")
print(f"Weighted Average: {weighted_avg} ± {combined_uncertainty_weighted}")
print(f"Equal Interval Average: {equal_interval_avg} ± {combined_uncertainty_equal}")
print(f"Linear Fit Slope: {slope} ± {slope_uncertainty}")
print(f"Linear Fit Intercept: {intercept} ± {intercept_uncertainty}")
print(f"Sum of Squared Residuals: {sum_squared_residuals}")

```


```


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

# データのセットアップ (例としてのデータを使用)
data = np.array([10.2, 9.8, 10.0, 9.7, 10.1, 10.3, 9.9])

### 1. 信頼区間 (Confidence Interval)
def confidence_interval(data, confidence=0.95):
    n = len(data)
    mean = np.mean(data)
    std_err = stats.sem(data)
    h = std_err * stats.t.ppf((1 + confidence) / 2., n-1)
    return mean, mean - h, mean + h

mean, ci_lower, ci_upper = confidence_interval(data)
print(f"信頼区間 (95%): 平均 = {mean:.2f}, CI = [{ci_lower:.2f}, {ci_upper:.2f}]")

### 2. 仮説検定 (Hypothesis Testing)
# 例として、平均が10であるかを検定
hypothesized_mean = 10.0
t_stat, p_value = stats.ttest_1samp(data, hypothesized_mean)
alpha = 0.05

print("\n仮説検定:")
print(f"t統計量 = {t_stat:.2f}, p値 = {p_value:.2f}")
if p_value < alpha:
    print("帰無仮説を棄却します (効果あり)")
else:
    print("帰無仮説を棄却しません (効果なし)")

### 3. 相関と回帰分析 (Correlation and Regression)
# 相関係数の計算
data_x = np.array([1, 2, 3, 4, 5, 6, 7])
corr_coeff = np.corrcoef(data_x, data)[0, 1]
print(f"\n相関係数 (Correlation Coefficient) = {corr_coeff:.2f}")

# 回帰分析 (線形回帰)
model = LinearRegression()
data_x_reshaped = data_x.reshape(-1, 1)
model.fit(data_x_reshaped, data)
slope = model.coef_[0]
intercept = model.intercept_

print(f"\n回帰分析 (Regression Analysis): y = {slope:.2f}x + {intercept:.2f}")

# 回帰直線をプロット
plt.scatter(data_x, data, color='blue', label='データポイント')
plt.plot(data_x, model.predict(data_x_reshaped), color='red', label=f'回帰直線: y = {slope:.2f}x + {intercept:.2f}')
plt.xlabel('独立変数 X')
plt.ylabel('従属変数 Y')
plt.legend()
plt.grid(True)
plt.show()

### 4. 不確かさの伝播 (Propagation of Uncertainty)
def propagation_of_uncertainty(u_X, u_Y, X, Y):
    Z = X * Y  # 例としてZ = X * Yの関数を仮定
    partial_X = Y
    partial_Y = X
    u_combined = math.sqrt((partial_X * u_X)**2 + (partial_Y * u_Y)**2)
    return Z, u_combined

# 入力データの例
X = 5.0
Y = 10.0
u_X = 0.1
u_Y = 0.2

Z, u_combined = propagation_of_uncertainty(u_X, u_Y, X, Y)
print(f"\n不確かさの伝播: Z = {Z}, 合成標準不確かさ = {u_combined:.2f}")

### 5. サンプリング (Sampling)
def sampling(data, sample_size):
    return np.random.choice(data, sample_size, replace=False)

sample_size = 3
sample = sampling(data, sample_size)
print(f"\nサンプリングされた標本 (Sample): {sample}")

# サンプリング誤差の計算
sample_mean = np.mean(sample)
sampling_error = np.abs(np.mean(data) - sample_mean)
print(f"サンプリングエラー (Sampling Error): {sampling_error:.2f}")
```


```



import numpy as np
import matplotlib.pyplot as plt

# 1. Parameters and Definitions
true_value = 100.0  # The true value of the measurement
measured_values = np.array([99.8, 100.1, 99.9, 100.2, 100.0])  # Sampled measurement values
calibration_uncertainty = 0.05  # Calibration uncertainty from a JCSS calibration

# 2. Calculate Measurement Errors
errors = measured_values - true_value

# 3. Evaluate Standard Uncertainty (Type A)
mean_value = np.mean(measured_values)
standard_deviation = np.std(measured_values, ddof=1)
standard_uncertainty_type_A = standard_deviation / np.sqrt(len(measured_values))

# 4. Combine with Calibration Uncertainty (Type B)
combined_uncertainty = np.sqrt(standard_uncertainty_type_A**2 + calibration_uncertainty**2)

# 5. Coverage Factor and Expanded Uncertainty
coverage_factor = 2  # For approximately 95% confidence level
expanded_uncertainty = coverage_factor * combined_uncertainty

# 6. Display Results
print(f"True Value: {true_value}")
print(f"Measured Values: {measured_values}")
print(f"Mean Value: {mean_value}")
print(f"Standard Uncertainty (Type A): {standard_uncertainty_type_A}")
print(f"Combined Uncertainty: {combined_uncertainty}")
print(f"Expanded Uncertainty: {expanded_uncertainty}")

# 7. Plotting the Results
plt.errorbar(np.arange(len(measured_values)), measured_values, yerr=expanded_uncertainty, fmt='o', label='Measurements with Uncertainty')
plt.axhline(y=true_value, color='r', linestyle='--', label='True Value')
plt.xlabel('Measurement Index')
plt.ylabel('Measured Value')
plt.title('Measurement Values with Expanded Uncertainty')
plt.legend()
plt.show()

```



```

import numpy as np
import matplotlib.pyplot as plt

# Constants and Parameters
standard_resistor_value = 10000  # Reference standard resistor value in ohms
measured_values = np.array([9998, 10002, 10001, 9999, 10000])  # Measured resistance values in ohms
leakage_current_effect = 0.002  # Leakage current effect from Johnson terminal in ohms
traceability_uncertainty = 0.05  # Traceability uncertainty in ohms
calibration_uncertainty = 0.03  # Calibration uncertainty in ohms

# Correct measured values by accounting for leakage current
corrected_values = measured_values - leakage_current_effect

# Calculate Measurement Errors
errors = corrected_values - standard_resistor_value

# Evaluate Standard Uncertainty (Type A)
mean_value = np.mean(corrected_values)
standard_deviation = np.std(corrected_values, ddof=1)
standard_uncertainty_type_A = standard_deviation / np.sqrt(len(corrected_values))

# Combine with other Uncertainty Factors (Type B)
combined_uncertainty = np.sqrt(standard_uncertainty_type_A**2 + traceability_uncertainty**2 + calibration_uncertainty**2)

# Coverage Factor and Expanded Uncertainty
coverage_factor = 2  # For approximately 95% confidence level
expanded_uncertainty = coverage_factor * combined_uncertainty

# Display Results
print(f"Standard Resistor Value: {standard_resistor_value} ohms")
print(f"Corrected Measured Values: {corrected_values} ohms")
print(f"Mean Value: {mean_value} ohms")
print(f"Standard Uncertainty (Type A): {standard_uncertainty_type_A} ohms")
print(f"Combined Uncertainty: {combined_uncertainty} ohms")
print(f"Expanded Uncertainty: {expanded_uncertainty} ohms")

# Plotting the Results
plt.errorbar(np.arange(len(corrected_values)), corrected_values, yerr=expanded_uncertainty, fmt='o', label='Corrected Measurements with Uncertainty')
plt.axhline(y=standard_resistor_value, color='r', linestyle='--', label='Standard Resistor Value')
plt.xlabel('Measurement Index')
plt.ylabel('Corrected Measured Value (ohms)')
plt.title('Corrected Measurement Values with Expanded Uncertainty')
plt.legend()
plt.show()

```



```
import numpy as np
import matplotlib.pyplot as plt

# 1. Prepare the dataset
x = np.array([1, 2, 3, 4])
y = np.array([2, 3, 5, 7])

# 2. Calculate the mean
mean_x = np.mean(x)
mean_y = np.mean(y)
print(f"Mean of x: {mean_x:.2f}")
print(f"Mean of y: {mean_y:.2f}")

# 3. Calculate the correlation coefficient (r)
numerator = np.sum((x - mean_x) * (y - mean_y))
denominator = np.sqrt(np.sum((x - mean_x) ** 2)) * np.sqrt(np.sum((y - mean_y) ** 2))
r = numerator / denominator
print(f"Correlation coefficient (r): {r:.2f}")

# 4. Calculate the slope (beta_1)
std_x = np.std(x, ddof=0)
std_y = np.std(y, ddof=0)
beta_1 = r * (std_y / std_x)
print(f"Slope (beta_1): {beta_1:.2f}")

# Calculate the intercept (beta_0)
beta_0 = mean_y - beta_1 * mean_x
print(f"Intercept (beta_0): {beta_0:.2f}")

# 5. Create the regression line
y_pred = beta_0 + beta_1 * x

# 6. Evaluate the model
ss_total = np.sum((y - mean_y) ** 2)
ss_residual = np.sum((y - y_pred) ** 2)
r_squared = 1 - (ss_residual / ss_total)
print(f"R-squared (R^2): {r_squared:.2f}")

# 7. Plotting
plt.scatter(x, y, color='blue', label='Data points')
plt.plot(x, y_pred, color='red', label='Regression line')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Linear Regression')
plt.show()



```