Пусть есть параметрическая модель $M\left( \theta \right)$, где $\theta$ - параметры.
Функция правдоподобия $L\left( X | M\left( \theta \right) \right)$ определят достоверность получения набора данных $X$ при заданном наборе параметров и данной модели.
Задача: определить такой набор параметров $\theta$, для которого функция принимает наибольшее значение.
По порядку производной:
Не использует производных $L$
Использует первую производную $\frac{\partial L}{\partial \theta_i}$ (градиент)
Использует вторые прозиводные $\frac{\partial^2 L}{\partial \theta_i \partial \theta_j}$ (гессиан)
(brute force)
import numpy as np
from scipy.optimize import minimize
def rosen(x):
"""The Rosenbrock function"""
return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
minimize(rosen, x0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True})
Optimization terminated successfully. Current function value: 0.000000 Iterations: 339 Function evaluations: 571
final_simplex: (array([[ 1. , 1. , 1. , 1. , 1. ], [ 1. , 1. , 1. , 1. , 1. ], [ 1. , 1. , 1. , 1.00000001, 1.00000001], [ 1. , 1. , 1. , 1. , 1. ], [ 1. , 1. , 1. , 1. , 1. ], [ 1. , 1. , 1. , 1. , 0.99999999]]), array([ 4.86115343e-17, 7.65182843e-17, 8.11395684e-17, 8.63263255e-17, 8.64080682e-17, 2.17927418e-16])) fun: 4.8611534334221152e-17 message: 'Optimization terminated successfully.' nfev: 571 nit: 339 status: 0 success: True x: array([ 1., 1., 1., 1., 1.])
Направление на максимум всегда в направлении градиента функции:
$$ \theta_{k+1} = \theta_k + \beta_k \nabla L $$
Модификация метода - метод сопряженных градиентов на самом деле требует второй производной.
Главная формула:
$$ L(\theta) = L(\theta_0) + \nabla L( \theta - \theta_0) + \frac{1}{2} (\theta-\theta_0)^T H (\theta-\theta_0) + o(\theta-\theta_0)$$
$$\nabla f(\theta_k) + H(\theta_k)(\theta_{k+1} - \theta_k) = 0$$
$$ \theta_{k+1} = \theta_k - H^{-1}(\theta_k)\nabla L(\theta_k) $$
Можно добавить выбор шага:
$$ \theta_{k+1} = \theta_k - \lambda_i H^{-1}(\theta_k)\nabla L(\theta_k) $$
from scipy import optimize
optimize.newton(lambda x: x**3 - 1, 1.5)
1.0000000000000016
Примеры:
minimize(method=’L-BFGS-B’)
В случае анализа спектров имеем:
$$ L(X | \theta) = \prod p_i (x_i | \theta)$$
Или:
$$\ln{ L(X | \theta)} = \sum \ln{ p_i (x_i | \theta)}$$
В случае нормальных распределений:
$$\ln{ L(X | \theta)} \sim \sum{ \left( \frac{y_i - \mu(x_i, \theta)}{\sigma_i} \right)^2 }$$
Пусть: $$r_i = \frac{y_i - \mu(x_i, \theta)}{\sigma_i}$$ $$J_{ij} = \frac{\partial r_i}{\partial \theta_j} = - \frac{\partial \mu(x_i, \theta)}{\sigma_i \partial \theta_j}$$
Тогда:
$$ \theta_{(k+1)} = \theta_{(k)} - \left( J^TJ \right)^{-1}J^Tr(\theta)$$
$$ \theta_{(k+1)} = \theta_{(k)} + \delta$$
$$ (J^TJ + \lambda I)\delta = J^Tr(\theta)$$
При этом $\lambda$ - фактор регуляризации, выбирается произвольным образом.
Идея: Есть некоторая статистика (функция данных) $f(x)$. Для оптимального решения среднее от этой функции по экспериментальным данным и по модели должны совпадать: $$ E_\theta(f(x)) = \sum_i{f(x_i)} $$
Можно показать, что оптимальная эффективность получается когда
$$ f = \frac{\partial \ln L}{\partial \theta} $$
В этом случае и если ошибки распределены по Гауссу или Пуассону, решение для оптмального $\theta$ можно получить как:
$$ \sum_{i}{\frac{\mu_{i}\left( \mathbf{\theta},E_{i} \right) - x_{i}}{\sigma_{i}^{2}}\left. \ \frac{\partial\mu_{i}\left( \mathbf{\theta},E_{i} \right)}{\partial\mathbf{\theta}} \right|_{\mathbf{\theta}_{\mathbf{0}}}} = 0. $$