Numpy polyfit: возможна ошибка масштабирования ковариационной матрицы?

Question

Мне сложно определить масштабирование ковариационной матрицы в numpy polyfit.

В документации я читал, что коэффициент масштабирования до go из немасштабированный до масштабированной ковариационной матрицы:

chi2 / sqrt(N - DOF).

В приведенном ниже коде кажется, что коэффициент масштабирования на самом деле

chi2 / DOF

Вот мой код

# Generate synthetically the data
# True parameters
import numpy as np

true_slope = 3
true_intercept = 7

x_data = np.linspace(-5, 5, 30)

# The y-data will have a noise term, to simulate imperfect observations
sigma = 1
y_data = true_slope * np.linspace(-5, 5, 30) + true_intercept
y_obs = y_data + np.random.normal(loc=0.0, scale=sigma, size=x_data.size)

# Here I generate artificially some unequal uncertainties 
# (even if there is no reason for them to be so)
y_uncertainties = sigma * np.random.normal(loc=1.0, scale=0.5*sigma, size=x_data.size)

# Make the fit
popt, pcov = np.polyfit(x_data, y_obs, 1, w=1/y_uncertainties, cov='unscaled')
popt, pcov_scaled = np.polyfit(x_data, y_obs, 1, w=1/y_uncertainties, cov=True)

my_scale_factor = np.sum((y_obs - popt[0] * x_data  - popt[1])**2 / y_uncertainties**2)\
                         / (len(y_obs)-2)

scale_factor =  pcov_scaled[0,0] / pcov[0,0]

Если я запускаю код, я вижу, что фактический коэффициент масштабирования равен chi2 / DOF, а не значение, указанное в документации. Это правда или я что-то упустил?

У меня еще вопрос. Почему предлагается использовать только обратную ошибку y-данных вместо квадрата обратной величины ошибок y-данных для весов в случае, если неопределенности распределены нормально?

Редактировать, чтобы добавить данные, сгенерированные запуском кода

x_data = array([-5.        , -4.65517241, -4.31034483, -3.96551724, -3.62068966,
   -3.27586207, -2.93103448, -2.5862069 , -2.24137931, -1.89655172,
   -1.55172414, -1.20689655, -0.86206897, -0.51724138, -0.17241379,
    0.17241379,  0.51724138,  0.86206897,  1.20689655,  1.55172414,
    1.89655172,  2.24137931,  2.5862069 ,  2.93103448,  3.27586207,
    3.62068966,  3.96551724,  4.31034483,  4.65517241,  5.        ])

y_obs = array([-7.27819725, -8.41939411, -3.9089926 , -5.24622589, -3.78747379,
   -1.92898727, -1.375255  , -1.84388812, -0.37092441,  0.27572306,
    2.57470918,  3.860485  ,  4.62580789,  5.34147103,  6.68231985,
    7.38242258,  8.28346559,  9.46008873, 10.69300274, 12.46051285,
   13.35049975, 13.28279961, 14.31604781, 16.8226239 , 16.81708308,
   18.64342284, 19.37375515, 19.6714002 , 20.13700708, 22.72327533])

y_uncertainties = array([ 0.63543112,  1.07608924,  0.83603265, -0.03442888, -0.07049299,
    1.30864191,  1.36015322,  1.42125414,  1.04099854,  1.20556608,
    0.43749964,  1.635056  ,  1.00627014,  0.40512511,  1.19638787,
    1.26230966,  0.68253139,  0.98055035,  1.01512232,  1.83910276,
    0.96763007,  0.57373151,  1.69358475,  0.62068133,  0.70030971,
    0.34648312,  1.85234844,  1.18687269,  1.23841579,  1.19741206])

С этими данными я получаю, что scale_factor = 1.6534129347542432, my_scale_factor = 1.653412934754234 и что «номинальный» масштабный коэффициент, указанный в документации. , т.е.

nominal_scale_factor = np.sum((y_obs - popt[0] * x_data  - popt[1])**2 /\  
                               y_uncertainties**2) / np.sqrt(len(y_obs) - len(y_obs) + 2)

имеет значение nominal_scale_factor = 32.73590595145554

PS. моя numpy версия - 1.18.5 3.7.7 (default, May 6 2020, 11:45:54) [MSC v.1916 64 bit (AMD64)]

Answer 1

Относительно документации numpy.polyfit:

По умолчанию ковариация масштабируется по chi2 / sqrt (N-dof), т.е. веса считаются ненадежными, за исключением относительного значения и все масштабируется таким образом, что уменьшенный chi2 равен единице.

Это похоже на ошибку в документации. Правильный коэффициент масштабирования для ковариации: chi_square/(N-M), где M - количество подходящих параметров, а N-M - количество степеней свободы. Похоже, что np.polyfit реализован правильно, потому что my_scale_factor и scale_factor согласованы.

Относительно вопроса, почему не «квадрат обратной величины ошибок y-данных»: полиномиальная аппроксимация или, в более общем смысле, аппроксимация методом наименьших квадратов включает решение вектора p в

A @ p = y

, где A - матрица (N, M) для N точек данных в y и M элементов в p, а каждый столбец в A - это полиномиальный член, вычисленный при соответствующих значениях x.

Решение минимизирует

    (SUM_j A[i, j] p[j] - y[i])^2
SUM -----------------------------
 i           sigma_y[i]^2

В вычислительном отношении самый дешевый способ вычисления это происходит путем умножения каждой строки в A и каждого значения y на соответствующее 1/sigma_y с последующим взятием стандартного решения уравнения A@p=y методом наименьших квадратов. Если пользователь предоставит обратные ошибки, это избавит процедуру подгонки от обработки деления на ноль и медленных операций возведения в квадрат - root.

Answer 2

Что касается первой части, я открыл проблему Github

https://github.com/numpy/numpy/issues/16842

Вывод по этому потоку - документация неверна, но функция ведет себя

Документация должна быть обновлена ​​до

По умолчанию ковариация масштабируется на chi2 / dof , т. е. предполагается, что веса ненадежны за исключением относительного значения, и все масштабируется так, что приведенная chi2 равна единице.