python - 如何用Python确定拟合参数的不确定性？

Question

我有以下 x 和 y 数据:

x   y
1.71    0.0
1.76    5.0
1.81    10.0
1.86    15.0
1.93    20.0
2.01    25.0
2.09    30.0
2.20    35.0
2.32    40.0
2.47    45.0
2.65    50.0
2.87    55.0
3.16    60.0
3.53    65.0
4.02    70.0
4.69    75.0
5.64    80.0
7.07    85.0
9.35    90.0
13.34   95.0
21.43   100.0

对于上面的数据，我试图将数据拟合成以下形式:

但是，存在与 x 和 y 相关的某些不确定性，其中 x 的不确定性为 x 的 50%，而 y 具有固定的不确定性。我试图用这个 uncertainties package 确定拟合参数的不确定性.但是，我在使用 scipy optimize 的曲线拟合函数进行曲线拟合时遇到了问题。我收到以下错误:

minpack.error: Result from function call is not a proper array of floats.

如何修复跟随错误并确定拟合参数(a、b 和 n)的不确定性？

MWE

from __future__ import division
import numpy as np
import re
from scipy import optimize, interpolate, spatial
from scipy.interpolate import UnivariateSpline
from uncertainties import unumpy


def linear_fit(x, a, b):
    return a * x + b


uncertainty = 0.5
y_error = 1.2
x = np.array([1.71, 1.76, 1.81, 1.86, 1.93, 2.01, 2.09, 2.20, 2.32, 2.47, 2.65, 2.87, 3.16, 3.53, 4.02, 4.69, 5.64, 7.07, 9.35, 13.34, 21.43])
x_uncertainty = x * uncertainty
x = unumpy.uarray(x, x_uncertainty)
y = np.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0])
y = unumpy.uarray(y, y_error)


n = np.arange(0, 5, 0.005)
coefficient_determination_on = np.empty(shape = (len(n),))
for j in range(len(n)):
    n_correlation = n[j]
    x_fit = 1 / ((x) ** n_correlation)
    y_fit = y
    fit_a_raw, fit_b_raw = optimize.curve_fit(linear_fit, x_fit, y_fit)[0]
    x_prediction = (fit_a_raw / ((x) ** n_correlation)) + fit_b_raw
    y_residual_squares = np.sum((x_prediction - y) ** 2)
    y_total_squares = np.sum((y - np.mean(y)) ** 2)
    coefficient_determination_on[j] = 1 - (y_residual_squares / y_total_squares)

Answer 1

让我首先说明这个问题不可能“很好地”解决，因为你想解决 a、b 和 n。这是因为对于一个固定的 n，你的问题承认一个封闭形式的解决方案，而如果你让 n 是自由的，它就不会，事实上这个问题可能有多个解决方案。因此，经典错误分析(例如 uncertanities 使用的错误分析)失效了，您必须求助于其他方法。

案例n已修复

如果 n 是固定的，你的问题是你调用的库不支持 uarray，所以你必须做一个解决方法。值得庆幸的是，线性拟合(在 l2 距离下)很简单 Linear least squares它承认一个封闭形式的解决方案，只需要用一个填充值然后解决 normal equations .

地点:

你可以这样做:

import numpy as np
from uncertainties import unumpy

uncertainty = 0.5
y_error = 1.2
n = 1.0

# Define x and y
x = np.array([1.71, 1.76, 1.81, 1.86, 1.93, 2.01, 2.09, 2.20, 2.32, 2.47, 2.65,
              2.87, 3.16, 3.53, 4.02, 4.69, 5.64, 7.07, 9.35, 13.34, 21.43])
# Take power of x values according to n
x_pow = x ** n
x_uncertainty = x_pow * uncertainty
x_fit = unumpy.uarray(np.c_[x_pow, np.ones_like(x)],
                      np.c_[x_uncertainty, np.zeros_like(x_uncertainty)])

y = np.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0,
              55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0])
y_fit = unumpy.uarray(y, y_error)

# Use normal equations to find coefficients
inv_mat = unumpy.ulinalg.pinv(x_fit.T.dot(x_fit))
fit_a, fit_b = inv_mat.dot(x_fit.T.dot(y_fit))

print('fit_a={}, fit_b={}'.format(fit_a, fit_b))

结果:

fit_a=4.8+/-2.6, fit_b=28+/-10

案例n未知

n 未知，你真的有麻烦了，因为问题是非凸的。在这里，线性误差分析(由 uncertainties 执行)将不起作用。

一种解决方案是执行 Bayesian inference , 使用像 pymc 这样的包.如果你对此感兴趣，我可以尝试写一篇文章，但它不会像上面那样干净。