0
1

実験の統計学

Last updated at Posted at 2024-08-14
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()



0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1