mirror of
https://github.com/fjosw/pyerrors.git
synced 2025-03-15 06:40:24 +01:00
Merge pull request #121 from fjosw/feat/num_diff_fit
Least square fit with numerical differentiation
This commit is contained in:
Fabian Joswig
GitHub
commit
889e24367d
2 changed files with 98 additions and 12 deletions
|
|
@ -9,8 +9,11 @@ import matplotlib.pyplot as plt
|
|||
from matplotlib import gridspec
|
||||
from scipy.odr import ODR, Model, RealData
|
||||
import iminuit
|
||||
from autograd import jacobian
|
||||
from autograd import jacobian as auto_jacobian
|
||||
from autograd import hessian as auto_hessian
|
||||
from autograd import elementwise_grad as egrad
|
||||
from numdifftools import Jacobian as num_jacobian
|
||||
from numdifftools import Hessian as num_hessian
|
||||
from .obs import Obs, derived_observable, covariance, cov_Obs
|
||||
|
||||
|
||||
|
|
@ -114,6 +117,8 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
|
|||
If True, a plot which displays fit, data and residuals is generated (default False).
|
||||
qqplot : bool
|
||||
If True, a quantile-quantile plot of the fit result is generated (default False).
|
||||
num_grad : bool
|
||||
Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
|
||||
'''
|
||||
if priors is not None:
|
||||
return _prior_fit(x, y, func, priors, silent=silent, **kwargs)
|
||||
|
|
@ -160,6 +165,8 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
|
|||
corrected by effects caused by correlated input data.
|
||||
This can take a while as the full correlation matrix
|
||||
has to be calculated (default False).
|
||||
num_grad : bool
|
||||
Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
|
@ -174,6 +181,13 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
|
|||
|
||||
x_shape = x.shape
|
||||
|
||||
if kwargs.get('num_grad') is True:
|
||||
jacobian = num_jacobian
|
||||
hessian = num_hessian
|
||||
else:
|
||||
jacobian = auto_jacobian
|
||||
hessian = auto_hessian
|
||||
|
||||
if not callable(func):
|
||||
raise TypeError('func has to be a function.')
|
||||
|
||||
|
|
@ -268,7 +282,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
|
|||
|
||||
fitp = out.beta
|
||||
try:
|
||||
hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel())))
|
||||
hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel())))
|
||||
except TypeError:
|
||||
raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
|
||||
|
||||
|
|
@ -277,7 +291,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
|
|||
chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((d[n_parms + m:].reshape(x_shape) - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2)
|
||||
return chisq
|
||||
|
||||
jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
|
||||
jac_jac_x = hessian(odr_chisquare_compact_x)(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
|
||||
|
||||
# Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
|
||||
try:
|
||||
|
|
@ -290,7 +304,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
|
|||
chisq = anp.sum(((d[n_parms + m:] - model) / dy_f) ** 2) + anp.sum(((x_f - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2)
|
||||
return chisq
|
||||
|
||||
jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f)))
|
||||
jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f)))
|
||||
|
||||
# Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
|
||||
try:
|
||||
|
|
@ -318,6 +332,11 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs):
|
|||
|
||||
x = np.asarray(x)
|
||||
|
||||
if kwargs.get('num_grad') is True:
|
||||
hessian = num_hessian
|
||||
else:
|
||||
hessian = auto_hessian
|
||||
|
||||
if not callable(func):
|
||||
raise TypeError('func has to be a function.')
|
||||
|
||||
|
|
@ -406,14 +425,15 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs):
|
|||
if not m.fmin.is_valid:
|
||||
raise Exception('The minimization procedure did not converge.')
|
||||
|
||||
hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(params))
|
||||
hess = hessian(chisqfunc)(params)
|
||||
hess_inv = np.linalg.pinv(hess)
|
||||
|
||||
def chisqfunc_compact(d):
|
||||
model = func(d[:n_parms], x)
|
||||
chisq = anp.sum(((d[n_parms: n_parms + len(x)] - model) / dy_f) ** 2) + anp.sum(((d[n_parms + len(x):] - d[:n_parms]) / dp_f) ** 2)
|
||||
return chisq
|
||||
|
||||
jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((params, y_f, p_f)))
|
||||
jac_jac = hessian(chisqfunc_compact)(np.concatenate((params, y_f, p_f)))
|
||||
|
||||
deriv = -hess_inv @ jac_jac[:n_parms, n_parms:]
|
||||
|
||||
|
|
@ -441,6 +461,13 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
|
|||
|
||||
x = np.asarray(x)
|
||||
|
||||
if kwargs.get('num_grad') is True:
|
||||
jacobian = num_jacobian
|
||||
hessian = num_hessian
|
||||
else:
|
||||
jacobian = auto_jacobian
|
||||
hessian = auto_hessian
|
||||
|
||||
if x.shape[-1] != len(y):
|
||||
raise Exception('x and y input have to have the same length')
|
||||
|
||||
|
|
@ -571,9 +598,9 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
|
|||
fitp = fit_result.x
|
||||
try:
|
||||
if kwargs.get('correlated_fit') is True:
|
||||
hess = jacobian(jacobian(chisqfunc_corr))(fitp)
|
||||
hess = hessian(chisqfunc_corr)(fitp)
|
||||
else:
|
||||
hess = jacobian(jacobian(chisqfunc))(fitp)
|
||||
hess = hessian(chisqfunc)(fitp)
|
||||
except TypeError:
|
||||
raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
|
||||
|
||||
|
|
@ -589,7 +616,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
|
|||
chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
|
||||
return chisq
|
||||
|
||||
jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f)))
|
||||
jac_jac = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))
|
||||
|
||||
# Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv
|
||||
try:
|
||||
|
|
|
|||
|
|
@ -83,8 +83,65 @@ def test_least_squares():
|
|||
assert math.isclose(pcov[i, i], betac[i].dvalue ** 2, abs_tol=1e-3)
|
||||
|
||||
|
||||
def test_least_squares_num_grad():
|
||||
x = []
|
||||
y = []
|
||||
for i in range(2, 5):
|
||||
x.append(i * 0.01)
|
||||
y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens"))
|
||||
|
||||
num = pe.fits.least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True)
|
||||
auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False)
|
||||
|
||||
assert(num[0] == auto[0])
|
||||
assert(num[1] == auto[1])
|
||||
|
||||
|
||||
def test_prior_fit_num_grad():
|
||||
x = []
|
||||
y = []
|
||||
for i in range(2, 5):
|
||||
x.append(i * 0.01)
|
||||
y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens"))
|
||||
|
||||
num = pe.fits.least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True, priors=y[:2])
|
||||
auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False, piors=y[:2])
|
||||
|
||||
|
||||
def test_least_squares_num_grad():
|
||||
x = []
|
||||
y = []
|
||||
for i in range(2, 5):
|
||||
x.append(i * 0.01)
|
||||
y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens"))
|
||||
|
||||
num = pe.fits.least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True)
|
||||
auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False)
|
||||
|
||||
assert(num[0] == auto[0])
|
||||
assert(num[1] == auto[1])
|
||||
|
||||
|
||||
assert(num[0] == auto[0])
|
||||
assert(num[1] == auto[1])
|
||||
|
||||
|
||||
def test_total_least_squares_num_grad():
|
||||
x = []
|
||||
y = []
|
||||
for i in range(2, 5):
|
||||
x.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens"))
|
||||
y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens"))
|
||||
|
||||
num = pe.fits.total_least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True)
|
||||
auto = pe.fits.total_least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False)
|
||||
|
||||
assert(num[0] == auto[0])
|
||||
assert(num[1] == auto[1])
|
||||
|
||||
|
||||
def test_alternative_solvers():
|
||||
dim = 192
|
||||
dim = 92
|
||||
x = np.arange(dim)
|
||||
y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim)
|
||||
yerr = 0.1 + 0.1 * np.random.rand(dim)
|
||||
|
|
@ -158,7 +215,7 @@ def test_correlated_fit():
|
|||
|
||||
|
||||
def test_fit_corr_independent():
|
||||
dim = 50
|
||||
dim = 30
|
||||
x = np.arange(dim)
|
||||
y = 0.84 * np.exp(-0.12 * x) + np.random.normal(0.0, 0.1, dim)
|
||||
yerr = [0.1] * dim
|
||||
|
|
@ -470,7 +527,7 @@ def test_correlated_fit_vs_jackknife():
|
|||
|
||||
|
||||
def test_fit_no_autograd():
|
||||
dim = 10
|
||||
dim = 3
|
||||
x = np.arange(dim)
|
||||
y = 2 * np.exp(-0.08 * x) + np.random.normal(0.0, 0.15, dim)
|
||||
yerr = 0.1 + 0.1 * np.random.rand(dim)
|
||||
|
|
@ -486,6 +543,8 @@ def test_fit_no_autograd():
|
|||
with pytest.raises(Exception):
|
||||
pe.least_squares(x, oy, func)
|
||||
|
||||
pe.least_squares(x, oy, func, num_grad=True)
|
||||
|
||||
with pytest.raises(Exception):
|
||||
pe.total_least_squares(oy, oy, func)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue