Я столкнулся со странной ситуацией в моем коде регрессии. Во-первых, позвольте мне показать мои модели (сумма экспонент)
$ f(t, \vec{p}) = \sum_n A_n\exp{-B_n t} $
import functools as ft
import numpy as np
import scipy.linalg as scilin
import scipy.optimize as opt
def my_func(t, *p):
val = 0
for j in range(0, len(p), 2):
val += p[j]*np.exp(-p[j+1]*t)
return val
def my_jac(t, *p):
jac = np.zeros(shape=(len(p), len(t)))
for j in range(0, len(p), 2):
jac[j] = np.exp(-p[j+1]*t)
jac[j+1] = -p[j]*t*np.exp(-p[j+1]*t)
return jac.T
def my_hess(t, *p):
hess = np.zeros((len(p), len(p), len(t)))
for j in range(0, len(p), 2):
hess[j+1, j] = -t*np.exp(-p[j+1]*t)
hess[j, j+1] = -t*np.exp(-p[j+1]*t)
hess[j+1, j+1] = p[j]*(t**2)*np.exp(-p[j+1]*t)
return hess
Ниже приведены следующие функции «стоимости», которые я хочу минимизировать
def chisquare(xo, xdata, ydata, covar, func):
# covar is covariance matrix
# residual is func(t, *xo) - ydata
res = func(xdata, *xo) - ydata
L = scilin.cholesky(covar, lower=True)
xres = scilin.solve_triangular(L, res, lower=True)
chisq = xres.dot(xres)
return chisq
def chi_jac(jac, xo, xdata, ydata, covar, func):
# I introduce a wrapper that takes in the jacobian of func, jac
# as a parameter. The wrapped function has the same arguments as chisquare
res = func(xdata, *xo) - ydata
L = scilin.cholesky(covar, lower=True)
# (J * Linv) * (Linv.H * res)
# lhs * rhs
lhs = scilin.solve_triangular(L, jac(xdata, *xo), lower=True)
rhs = scilin_solve_triangular(L, res)
chiJ = lhs.T.dot(rhs)
return 2*chiJ
def chi_hess(jac, hess, xo, xdata, ydata, covar, func):
# I introduce a wrapper that takes in the jacobian and hessian of func, jac and hess,
# as a parameter. The wrapped function has the same arguments as chisquare
res = func(xdata, *xo) - ydata
L = scilin.cholesky(covar, lower=True)
# (a) rhs: Solve L*H = H_l; L*res = xres
xres = scilin.solve_triangular(L, res, lower=True)
Hl = scilin.solve_triangular(L, hess(xdata, *xo).T, lower=True)
rhs = 2*Hl.T.dot(xres)
# (b) lhs: solve L*J = J_l
Jl = scilin.solve_triangular(L, jac(xdata, *xo), lower=True)
lhs = 2*Jl.T.dot(Jl)
return lhs + rhs
Ниже у меня есть две реализации изгиба кривой. В одном из них используется модуль curve_fit scipy.optimize, а в другом - модуль минимизации scipy.optimize.
def minsolve(xdata, ydata, func, xo, jac, hess):
my_args = (xdata, ydata, covar, func)
wrap_jac = ft.partial(chi_jac, jac)
wrap_hess = ft.partial(chi_hess, jac, hess)
res = opt.minimize(chisquare, xo, args=my_args, jac=wrap_jac,
hess=wrap_hess, method='trust-exact')
if not res.success:
res_message = f"Status: {res.status};\nMessage: {res.message}"
raise Exception(res_message)
return res.x
# for curve_fit I just use the module directly)
popt, pcov = opt.curve_fit(func, xdata, ydata, p0=xo, sigma=covar,
method='trf', tr_solver='exact')
При использовании кривой_fit мой алгоритм завершается, и я получаю как popt, так и pcov. Тем не менее, когда я запускаю мою minsolve реализацию optimize.minimize, успех не удался. В частности, я получаю следующее сообщение:
Предупреждение. Неправильная аппроксимация вызвала сбой в прогнозировании улучшения.
Это странно, потому что я использую один и тот же алгоритм Trust-Region для обоих модулей в одном и том же набор данных.
Причина, по которой я хочу использовать модуль минимизации, заключается в том, что я хочу изменить свою функцию $ \ chi ^ 2 $, чтобы она могла принимать такие вещи, как байесовские приоры для моего анализа. Пока что я застрял в этой проблеме.
Есть ли ошибка в моей реализации моих функций эффективной стоимости, якобиан, гессиан?
Ребята, есть ли у вас какие-либо рекомендации относительно того, что я должен сделать, чтобы все заработало?
Кстати, есть ли способ заставить латекс рендериться здесь? Я хотел бы показать функции, которые я использую в форме LaTeX, чтобы сделать вещи немного яснее.
РЕДАКТИРОВАТЬ: Дополнительная информация
Источник сообщения об ошибке scipy .optimize.minimize происходит от алгоритма поиска строки в точной области доверия.
С другой стороны, критерий сходимости 'gtol' соответствует scipy.optimize.curve_fit