我曾经在计算表中运行DROITEREG
函数。下面是一个例子:
左上角是数据,底部是函数DROITEREG
的结果,是一个2乘5的表。我写了几个细胞的标签。a和b是线性回归的参数,u(a)和u(b)是a和b上的不确定性。我想从一个numpy函数计算这些不确定性。在
我从curve_fit
函数成功:
import numpy as np
from scipy.stats import linregress
from scipy.optimize import curve_fit
data_o = """
0.42 2.0
0.97 5.0
1.71 10.0
2.49 20.0
3.53 50.0
3.72 100.0
"""
vo, So = np.loadtxt(data_o.split("\n"), unpack=True)
def f_model(x, a, b):
return a * x + b
popt, pcov = curve_fit(
f=f_model, # model function
xdata=1 / So, # x data
ydata=1 / vo, # y data
p0=(1, 1), # initial value of the parameters
)
# parameters
print(popt)
# uncertaintes :
print(np.sqrt(np.diag(pcov)))
输出如下,结果与DROITEREG
得到的结果一致:
但这并不完全令人满意,因为这应该很容易从一个简单的最小二乘函数得到。所以我尝试使用polyfit
。在
(a, b), Mcov = np.polyfit(1 / So, 1 / vo, 1, cov=True)
print("a = ", a, " b = ", b)
print("SSR = ", sum([(y - (a * x + b))**2 for x, y in zip(1 / So, 1 / vo)]))
print("Cov mat\n", Mcov)
print("Cov mat diag ", np.diag(Mcov))
print("sqrt 1/2 cov mat diag ", np.sqrt(0.5 * np.diag(Mcov)))
输出为:
a = 4.35522612104 b = 0.186297716685
SSR = 0.00398117627681
Cov mat
[[ 0.01144455 -0.00167853]
[-0.00167853 0.00057795]]
Cov mat diag [ 0.01144455 0.00057795]
sqrt 1/2 cov mat diag [ 0.07564571 0.01699926]
最后,我注意到来自polyfit
的Mcov矩阵是curve_fit
的pcov矩阵的2倍。如果我尝试用更大次数的多项式进行拟合,我发现因子等于参数的个数。在
我没有成功地使用linregress
,因为我不知道如何获得参数估计的协方差矩阵。同样,我也成功地使用了scipy.odr
,但它比上面的解决方案更复杂,这对于一个简单的线性回归来说也是如此。也许我遗漏了一些东西,因为我对统计学不太在行,而且我不真正理解这个协方差矩阵的意义。在
因此,我想知道的是获得线性回归参数和相关不确定性的最简单方法(相关系数也是一个很好的点,但它更容易计算)。例如,我的主要目标是给化学或物理专业的学生一个简单的方法来做这个线性回归,并计算与这个模型相关的参数。在
目前没有回答
相关问题 更多 >
编程相关推荐