Общий случай

Постановка задачи

Пусть есть параметрическая модель $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)

• Строим сетку и ищем на ней максимум.
• Возможен только для одномерных, максимум двумерных задач.
• Точность ограничена размером ячкйки сетки.

Симплекс методы

1. Строим многоугольник в пространстве параметров с $n+1$ вершинами, где $n$ - размерность пространства.
2. Орпделеляем значения функции в каждой вершине.
3. Находим вершину с наименьшим значением и двигаем ее к центру масс многоугольника.

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 $$

• Не понятно, как определять $\beta$
• Не понятно, когда останавливаться.

Модификация метода - метод сопряженных градиентов на самом деле требует второй производной.

Вторые производные

Главная формула:

$$ 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

Методы с переменной метрикой

• Вычислять $\nabla L$ и $H$ очень дорого
• Давайте вычислять их итеративно.

Примеры:

• MINUIT
• scipy 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. $$