Как вычислить ошибки стандартного отклонения с помощью scipy.optimize.least_squares

Question

Я сравниваю установку с optimize.curve_fit и optimize.least_squares. С curve_fit я получаю ковариационная матрица pcov в качестве выхода, и я могу вычислить стандартные ошибки отклонения для моих подогнанных переменных, что:

perr = np.sqrt(np.diag(pcov)) 

Если я фитинг с least_squares, я не получаю никакого вывода ковариационной матрицы и Я не могу вычислить ошибки стандартного отклонения для моих переменных.

Вот мой пример:

#import modules 
import matplotlib 
import numpy as np 
import matplotlib.pyplot as plt 
from scipy.optimize import curve_fit 
from scipy.optimize import least_squares 

noise = 0.5 
N = 100 
t = np.linspace(0, 4*np.pi, N) 

# generate data 
def generate_data(t, freq, amplitude, phase, offset, noise=0, n_outliers=0, random_state=0): 
    #formula for data generation with noise and outliers 
    y = np.sin(t * freq + phase) * amplitude + offset 
    rnd = np.random.RandomState(random_state) 
    error = noise * rnd.randn(t.size) 
    outliers = rnd.randint(0, t.size, n_outliers) 
    error[outliers] *= 10 
    return y + error 

#generate data 
data = generate_data(t, 1, 3, 0.001, 0.5, noise, n_outliers=10) 

#initial guesses 
p0=np.ones(4) 
x0=np.ones(4) 

# create the function we want to fit 
def my_sin(x, freq, amplitude, phase, offset): 
    return np.sin(x * freq + phase) * amplitude + offset 

# create the function we want to fit for least-square 
def my_sin_lsq(x, t, y): 
    # freq=x[0] 
    # phase=x[1] 
    # amplitude=x[2] 
    # offset=x[3] 
    return (np.sin(t*x[0]+x[2])*x[1]+ x[3]) - y 

# now do the fit for curve_fit 
fit = curve_fit(my_sin, t, data, p0=p0) 
print 'Curve fit output:'+str(fit[0]) 

#now do the fit for least_square 
res_lsq = least_squares(my_sin_lsq, x0, args=(t, data)) 
print 'Least_squares output:'+str(res_lsq.x) 


# we'll use this to plot our first estimate. This might already be good enough for you 
data_first_guess = my_sin(t, *p0) 

#data_first_guess_lsq = x0[2]*np.sin(t*x0[0]+x0[1])+x0[3] 
data_first_guess_lsq = my_sin(t, *x0) 

# recreate the fitted curve using the optimized parameters 
data_fit = my_sin(t, *fit[0]) 
data_fit_lsq = my_sin(t, *res_lsq.x) 

#calculation of residuals 
residuals = data - data_fit 
residuals_lsq = data - data_fit_lsq 
ss_res = np.sum(residuals**2) 
ss_tot = np.sum((data-np.mean(data))**2) 
ss_res_lsq = np.sum(residuals_lsq**2) 
ss_tot_lsq = np.sum((data-np.mean(data))**2) 

#R squared 
r_squared = 1 - (ss_res/ss_tot) 
r_squared_lsq = 1 - (ss_res_lsq/ss_tot_lsq) 
print 'R squared curve_fit is:'+str(r_squared) 
print 'R squared least_squares is:'+str(r_squared_lsq) 

plt.figure() 
plt.plot(t, data) 
plt.title('curve_fit') 
plt.plot(t, data_first_guess) 
plt.plot(t, data_fit) 
plt.plot(t, residuals) 

plt.figure() 
plt.plot(t, data) 
plt.title('lsq') 
plt.plot(t, data_first_guess_lsq) 
plt.plot(t, data_fit_lsq) 
plt.plot(t, residuals_lsq) 

#error 
perr = np.sqrt(np.diag(fit[1])) 
print 'The standard deviation errors for curve_fit are:' +str(perr) 

Я был бы очень благодарен за любую помощь, наилучшими пожеланиями

пса: Я много входа получил из этого источника и использовал часть кода Robust regression

Answer 1

Результат optimize.least_squares имеет параметр внутри него, называемый jac. Из documentation:

JAC: ndarray, разреженной матрицей или LinearOperator, форме (т, п)

Модифицированный матрицу Якоби при решении, в том смысле, что J^т является приближение Гаусса-Ньютона из гессиан функции стоимости. Тип такой же, как тот, который используется алгоритмом.

Это может быть использовано для оценки матрицы ковариации параметров по следующей формуле: Sigma = (J'J)^- 1. Из исходного кода optimize.curve_fit они также масштабируют эту оценку MSE остатков. В коде это будет выглядеть как

J = res_lsq.jac 
cov = np.linalg.inv(J.T.dot(J)) * (residuals_lsq**2).mean() 

Надеюсь, это поможет.