Julia vs R vs Python: простая оптимизация

Автор оригинала: ZJ.

2019-04-05 Обновление: в предыдущей версии этого сообщения были некоторые серьезные опасения по поводу проблемы задержки компиляции с пакетом Julia Optim.jl. Проблема была решена, и пакет Julia на самом деле довольно эффективен.

В этом посте я попытаюсь сравнить и противопоставить Julia, R и Python с помощью простой задачи оптимизации максимального правдоподобия, которая мотивирована проблемой из области кредитного риска и более подробно обсуждается в этом посте .

Примечание – уровень опыта автора на момент написания статьи

TL;DR

Для такой простой задачи оптимизации R, Julia и Python/SciPy будут все делать грамотную работу , так что явного победителя нет. Тем не менее, у Джулии есть преимущество, поскольку она имеет лучший результат и имеет лучшие способы работы с усеченными распределениями

Приведенные наблюдения Q 1 , Q 2 , . . . , Q n Q_1,\, Q_2,\, …,\, Q_n Q 1 , Q 2 , . . . , Q |/n , мы стремимся найти параметры μ \mu μ и σ \sigma σ , которые оптимизируют эту функцию правдоподобия

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) для усеченных нормальных распределений.

Джулия решение

Ниже приведена моя реализация Julia с использованием Optim.jl В Julia можно использовать символы в именах переменных, поэтому я использовал μ σ \mu\sigma μ σ в качестве имени переменной. У Джулии также есть популярный пакет под названием JuMP.jl для задач оптимизации. Однако пакет Optim.jl более чем адекватен для такой простой задачи, и я буду рассматривать JuMP.jl только в будущем, когда буду иметь дело с более сложными задачами оптимизации.

Я отметил, что время, необходимое для выполнения первого запуска задачи оптимизации, составляет менее 1 секунды. Это значительное улучшение по сравнению с 7,5 секундами в более ранней версии этого поста!

Ниже я жестко закодировал значения 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

Что хорошего было в Джулии?

Что было не так хорошо в Джулии?

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?

Для такой простой задачи оптимизации вы не можете ошибиться, выбрав любой из трех языков/экосистем. Все три языка/экосистемы позволяют вам определять замыкания, которые я нахожу хорошим способом генерировать функции, которые встраивают необходимые данные, оставляя вам только функцию с параметрами распределения в качестве единственных аргументов. Для меня Джулия явно лучший выбор, хотя и не с большим отрывом.