自动计算不确定性传播：回归替代方案？

Question

我最近完成了一个光学实验室。我测量了入射角 i1 和折射角 i2。我将这些数据记录在电子表格中并绘制了 sinu2061i2=fn(sinu2061i1)。我的不确定性为 1∘。

在实验室说明中，我们被要求使用Regressi，一个有点过时的软件，界面不直观（但非常强大的软件）

两列用于角度，两列用于不确定性，两列用于计算的正弦值，并且......等等：两列自动生成的列用于正弦值的不确定性！这很棒，据我所知，有没有其他带有 GUI 的免费软件可以提供此功能。只需点击几下，我就可以定义回归表达式（这里是一个简单的 y=a⋅x，但它可以更复杂），就是这样！

在 Python 中，等效的设置如下所示：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# python -m venv venv
# source venv/bin/activate
# pip install uncertainties numpy matplotlib scipy
import numpy as np
import matplotlib.pyplot as plt
import uncertainties.unumpy as unp
from scipy.optimize import curve_fit

# Initial data (angles of incidence and refraction)
i1 = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80])
i2 = np.array([0, 6.5, 13, 19.5, 25.5, 31, 35.5, 39, 41.5])
σ = 1

# Convert to objects with uncertainties
sin_i1 = unp.sin(unp.radians(unp.uarray(i1, σ)))
sin_i2 = unp.sin(unp.radians(unp.uarray(i2, σ)))

# i₁ in degrees
# i₂ in degrees
# Uncertainty in degrees
# Calculate sin(i₁) with uncertainties
# Calculate sin(i₂) with uncertainties

# Extract nominal values and uncertainties
sin_i1_vals = unp.nominal_values(sin_i1)
# Central values of sin(i₁)
sin_i2_vals = unp.nominal_values(sin_i2)
# Central values of sin(i₂)
sin_i1_err = unp.std_devs(sin_i1)
# Uncertainties of sin(i₁)
sin_i2_err = unp.std_devs(sin_i2)
# Uncertainties of sin(i₂)

# Plot data with error bars
plt.errorbar(sin_i1_vals, sin_i2_vals, xerr=sin_i1_err, yerr=sin_i2_err, fmt=' g', label='Measurements')

# Linear regression model: f(x) = a*x (intercept b = 0)
def linear_model(x, a):
    return a * x

# Regression with curve_fit
popt, pcov = curve_fit(linear_model, sin_i1_vals, sin_i2_vals) # pcov: parameter covariance, not used but noted...
a_opt = popt[0] # Regression coefficient

# Plot the regression
plt.plot(sin_i1_vals, linear_model(sin_i1_vals, a_opt), 'tab:blue', label=f"Regression: sin i₂ = {a_opt:.4f} * sin i₁")

# Formatting and display
plt.title('sin i₂ = f(sin i₁)')
plt.xlabel('sin i₁')
plt.ylabel('sin i₂')
plt.xlim([0, 1])
plt.ylim([0, 0.7])
plt.xticks(np.arange(0, 1.1, 0.1))
plt.yticks(np.arange(0, 0.8, 0.1))
plt.grid(True)
plt.legend()
plt.show()

（顺便说一句，如果您对我如何改进这个 Python 脚本有任何建议，我愿意接受想法！）

在 R（我对 R 的了解甚至比 Python 还要少）或 Octave 中也可能是可能的。

但是有件事已经困扰我好几天了。回归似乎主要在法国使用，那么其他国家的学生使用什么软件来自动计算实验室工作的不确定性传播？研究人员使用什么？ （我假设他们使用 Origin Pro，因为他们不介意它是否免费。）

所以，我一直在寻找替代方案。

首先，我尝试了SciDaVis/Alphaplot，我非常喜欢它......但我找不到如何自动传播不确定性。当然，我可以手动为 σsinu2061θ=∣cosu2061(θ)∣×σθ 创建一列，但是要找到这个公式需要工作，现在我知道软件会自动完成它，我不想不再自己做！

然后我尝试了labplot2，Gnuplot，LibreOffice Calc，GeoGebra，Dataplot，JASP，Jamovi。在这些软件选项中我都找不到这个功能🤔???

我不敢相信 Regressi 是世界上唯一提供自动不确定性传播的免费软件。如果是这样的话，为什么？？？非法国学生使用什么？为什么回归没有变得更广为人知并被认为是必不可少的？

我一定错过了什么，你能帮我吗？

Answer 1

设法找到了脚本，但是我遇到了一个小问题，因为我写的脚本考虑了

a*x + b

并且错误主要来自于传播梯度的行

J

。

我还通过添加随机噪声（以 0 为中心）稍微夸大了误差。

函数和导入在这里：

%matplotlib qt
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import numdifftools as nd
import uncertainties.unumpy as unp
def LinearRegressionModel(x, a, b):
    ''' Simple linear regression as a*x + b '''
    return a * x  + b # Define the linear regression model
def proxy(r, x):
    ''' Proxy function for LinearRegressionModel '''
    return LinearRegressionModel(x, *r)  # Apply the model to x with parameters r
def propagation(J, Cp):
    ''' Define the propagation function to compute the covariance '''
    return J @ Cp @ J.T  # Covariance propagation
def error(x, p, Cp):
    ''' Define the error function '''
    gradient = nd.Gradient(lambda r: proxy(r, x))  # gradient function of proxy w.r.t. parameters
    J = gradient(p)  # get gradient at model parameters p
    Cy = np.apply_along_axis(lambda row: propagation(row, Cp), 1, J) # propagate each row of the gradient
    sy = np.sqrt(Cy) # Calculate standard errors
    return sy

这就是例子：

i1 = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80])
i2 = np.array([0, 6.5, 13, 19.5, 25.5, 31, 35.5, 39, 41.5])
σ = 1
sin_i1 = unp.sin(unp.radians(unp.uarray(i1, σ)))
sin_i2 = unp.sin(unp.radians(unp.uarray(i2, σ)))
sin_i2 += np.random.normal(0, 0.1, len(sin_i2)) # inject some noise to exagerate, when removed it's not possible to see the results...
sin_i1_vals = unp.nominal_values(sin_i1)
sin_i2_vals = unp.nominal_values(sin_i2)
sin_i1_err = unp.std_devs(sin_i1)
sin_i2_err = unp.std_devs(sin_i2)
plt.errorbar(sin_i1_vals, sin_i2_vals, xerr=sin_i1_err, yerr=sin_i2_err, fmt=' g', label='Measurements')
popt, pcov = curve_fit(LinearRegressionModel, sin_i1_vals, sin_i2_vals)
sy = error(sin_i1_vals, popt, pcov)
predictions = LinearRegressionModel(sin_i1_vals, *popt)
plt.plot(sin_i1_vals, predictions)
plt.fill_between(sin_i1_vals, # x
                 predictions+sy, # upper = y + error
                 predictions-sy, # lower = y - error
                 alpha = 0.1, # make transparent
                )
plt.title('sin i₂ = f(sin i₁)')
plt.xlabel('sin i₁')
plt.ylabel('sin i₂')
plt.xlim([0, 1])
plt.ylim([0, 0.7])
plt.xticks(np.arange(0, 1.1, 0.1))
plt.yticks(np.arange(0, 0.8, 0.1))

我向

y

变量注入了一些噪声，以稍微夸大错误，否则它并不真正可见。

您可以像使用

sy

或

sin_i1_err

一样使用

sin_i2_err

。

这是一个值得考虑的好来源。