2019-04-05 Обновление: в предыдущей версии этого сообщения были некоторые серьезные опасения по поводу проблемы задержки компиляции с пакетом Julia Optim.jl. Проблема была решена, и пакет Julia на самом деле довольно эффективен.
В этом посте я попытаюсь сравнить и противопоставить Julia, R и Python с помощью простой задачи оптимизации максимального правдоподобия, которая мотивирована проблемой из области кредитного риска и более подробно обсуждается в этом посте .
Примечание – уровень опыта автора на момент написания статьи
R | 9 лет |
Юля | 6 месяцев |
Питон | начинающий |
TL;DR
Для такой простой задачи оптимизации R, Julia и Python/SciPy будут все делать грамотную работу , так что явного победителя нет. Тем не менее, у Джулии есть преимущество, поскольку она имеет лучший результат и имеет лучшие способы работы с усеченными распределениями
Sub 1s ждет первого запуска из – за задержки компиляции | Простота в использовании; Замечательная поддержка усеченных распределений; Легкий для понимания вывод | Юля |
Вывод модели не так читаем как у Джулии | Простота в использовании | R |
Вывод модели не так читаем как у Джулии | Простота в использовании | Питон |
Приведенные наблюдения
L = ∏ ( ϕ ( Q i , μ , σ ) Φ ( максимум Q t , μ , σ ) ) L = \prod\left(\frac{\phi(Q_i,\mu,\sigma)}{\Phi(\max Q_t,\mu,\sigma)}\right) L = ∏ ( Φ ( максимум Q t , μ , σ ) ϕ ( Q i , μ , σ ) )
часто вместо этого мы пытаемся оптимизировать логарифмическую вероятность
бревно L = l = ( ∑ i бревно ϕ ( Q i , μ , σ ) ) − n бревно Φ ( максимум Q t , μ , σ ) \log = \left(\sum_i \log\phi(Q_i,\mu,\sigma)\right) – n\log\Phi(\max Q_t,\mu,\sigma) Ло g L = l = ( i ∑ Ло g ϕ ( Q i , μ , σ ) ) − n Ло g Φ ( максимум Q t , μ , σ )
Говоря статистическим языком, это оценка максимального правдоподобия (MLE) для усеченных нормальных распределений.
В этой статье Ergashev et al., 2016 (за платной стеной) выведено необходимое условие существования решения и выведено уравнение в одном параметре, которое при решении (с использованием простых методов uniroot) может восстанавливать как μ\muμ, так и σ\sigmaσ. Из моего тестирования применение формулы Эргашева дает примерно 50-кратную скорость до решения R. Однако в этом посте в блоге я стремлюсь сравнить и противопоставить функцию оптимизации в Julia vs. R против Python, и поэтому я решил не реализовывать методы Эргашева. |
Джулия решение
Ниже приведена моя реализация Julia с использованием Optim.jl
В Julia можно использовать символы в именах переменных, поэтому я использовал
Я отметил, что время, необходимое для выполнения первого запуска задачи оптимизации, составляет менее 1 секунды. Это значительное улучшение по сравнению с 7,5 секундами в более ранней версии этого поста!
Если вы хотите, чтобы использование Optim вообще не имело задержки и чтобы вызов optim был эффективным, вы можете попробовать PackageCompiler.jl предварительно скомпилировать пакет Optim.jl с помощью using PackageCompiler; compile_incremental(:Optim,). После завершения в ОТВЕТЕ должны быть инструкции о том, как использовать скомпилированный образ для сверхскоростного использования вызовов Optim и optim |
Ниже я жестко закодировал значения Q_t
, которые я хотел бы использовать в оценке MLE.
using Distributions, Optim # hard-coded data\observations odr=[0.10,0.20,0.15,0.22,0.15,0.10,0.08,0.09,0.12] Q_t = quantile.(Normal(0,1), odr) # return a function that accepts `[mu, sigma]` as parameter function neglik_tn(Q_t) maxx = maximum(Q_t) f(μσ) = -sum(logpdf.(Truncated(Normal(μσ[1],μσ[2]), -Inf, maxx), Q_t)) f end neglikfn = neglik_tn(Q_t) # optimize! # start searching @time res = optimize(neglikfn, [mean(Q_t), std(Q_t)]) # 0.8 seconds @time res = optimize(neglikfn, [mean(Q_t), std(Q_t)]) # 0.000137 seconds # the \mu and \sigma estimates Optim.minimizer(res) # [-1.0733250637041452,0.2537450497038758] # or # use `fieldnames(res)` to see the list of field names that can be referenced via . (dot) res.minimizer # [-1.0733250637041452,0.2537450497038758]
что дало мне следующий вывод, который, безусловно, является наиболее описательным из стиха Julia/R/Python.
Results of Optimization Algorithm * Algorithm: Nelder-Mead * Starting Point: [-1.1300664159893685,0.22269345618402703] * Minimizer: [-1.0733250637041452,0.2537450497038758] * Minimum: -1.893080e+00 * Iterations: 28 * Convergence: true * √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: true * Reached Maximum Number of Iterations: false * Objective Calls: 59
Что хорошего было в Джулии?
Усеченный(DN, нижний, верхний) – это удивительно простой способ определения усеченного распределения | |
Функция logpdf, которая работает для любого дистрибутива | |
Очень четкий и информативный вывод |
Что было не так хорошо в Джулии?
Приемлемый | Задержка компиляции менее 1 секунды |
R решение
R имеет пакет truncnorm
для работы с усеченными нормалями.
odr=c(0.10,0.20,0.15,0.22,0.15,0.10,0.08,0.09,0.12) x = qnorm(odr) library(truncnorm) neglik_tn = function(x) { maxx = max(x) resfn = function(musigma) { -sum(log(dtruncnorm(x, a = -Inf, b= maxx, musigma[1], musigma[2]))) } resfn } neglikfn = neglik_tn(x) system.time(res <- optim(c(mean(x), sd(x)), neglikfn)) res
что и дало выход
$par [1] -1.0733354 0.2537339 $value [1] -1.89308 $counts function gradient 55 NA $convergence [1] 0 $message NULL
Что хорошего было в Р?
Имеет пакет для усеченного нормального |
Что было не так хорошо в R?
Пара незначительных вещей: нет плотности журнала для усеченного нормального; и нет простого способа определить усеченное распределение для произвольных распределений; разреженный вывод (печать) |
Несмотря на то, что у меня нет опыта работы с Python, простые поиски в Google позволили мне придумать это решение. Я использовал дистрибутив Anaconda, который избавил меня от многих хлопот с точки зрения установки пакетов, так как все пакеты, которые мне нужны, предварительно установлены.
import numpy as np from scipy.optimize import minimize from scipy.stats import norm # generate the data odr=[0.10,0.20,0.15,0.22,0.15,0.10,0.08,0.09,0.12] Q_t = norm.ppf(odr) maxQ_t = max(Q_t) # define a function that will return a return to optimize based on the input data def neglik_trunc_tn(Q_t): maxQ_t = max(Q_t) def neglik_trunc_fn(musigma): return -sum(norm.logpdf(Q_t, musigma[0], musigma[1])) + len(Q_t)*norm.logcdf(maxQ_t, musigma[0], musigma[1]) return neglik_trunc_fn # the likelihood function to optimize neglik = neglik_trunc_tn(Q_t) # optimize! res = minimize(neglik, [np.mean(Q_t), np.std(Q_t)]) res
и вот такие результаты
fun: -1.8930804441641282 hess_inv: array([[ 0.01759589, 0.00818596], [ 0.00818596, 0.00937868]]) jac: array([ -3.87430191e-07, 3.33786011e-06]) message: 'Optimization terminated successfully.' nfev: 40 nit: 6 njev: 10 status: 0 success: True x: array([-1.07334252, 0.25373624])
Что хорошего было в Питоне?
Легко гуглить даже для начинающих |
Что было не так хорошо в Python?
Результат мог бы быть лучше, но это незначительный момент |
Для такой простой задачи оптимизации вы не можете ошибиться, выбрав любой из трех языков/экосистем. Все три языка/экосистемы позволяют вам определять замыкания, которые я нахожу хорошим способом генерировать функции, которые встраивают необходимые данные, оставляя вам только функцию с параметрами распределения в качестве единственных аргументов. Для меня Джулия явно лучший выбор, хотя и не с большим отрывом.