注册用户享全站资源 并成为粉丝 不定时福利发放
 
非盈利性学习社区
分享网站

分享网站

让更多人可以获取免费教程

合作广告

合作广告

最大优惠

人事面试

人事面试

套路太深

部分企业人事部收到求职者简历后,会预先进行一轮电话面试,来考察求职者的综合素质,因此了解懂得求职面试的基本技巧,将...

1、着装要适当  穿着不一定要名贵,但一定要合体、干净、整洁,而且颜色和图案的搭配一定要协调。鞋子应该是舒服而又引...

一、隐瞒真实个人资料的不诚实者简历是求职的第一步,只有面试官对你的简历有兴趣才会通知你面试。在简历中适当地突出个人...

一、首先你得已经成为公司里“最好”的程序员,或者你已经找不到可作为老师和导师的人关于这一点,很多人都会过度自信,所...

面试是求职的必经阶段,这个过程中,人的性格不能说决定一切,但是却还是有一些影响的。性格活泼的人,能迅速的融入环境,...

程序员段子

程序员段子

乐一乐

公司高层公司副总A:咱们开个会研究一下这个事情怎么处理。公司副总B:如果老板没有救成功,下任是谁呢?会不会影响公司...

FBI WARNING:以下内容纯属恶搞,如有雷同,那我也没办法。(其实程序猿还是很有前途啊!!!)尼玛,扔个肥皂...

放心吧,程序猿是不会有外遇的产品经理是大家的公敌……Chrome 浏览器原来是这么来的TMD,瞧我这火爆脾气,发起...

笑话秘籍入门篇大学时候,刚学会上网,于是在网上搜A片,然后搜出个PHP大全.rmvb。以为是拍黄片大全,谁知道竟然...

人们常说程序员的生活枯燥为人刻板,其实这是你不懂程序员,代码之外,这些高智商的人幽默有趣,论坛常常是他们展现才华的...

优化过程 PK :Julia 能打败 Python 和 R 笑到最后吗?
 二维码
发表时间:2018-10-30 17:51

在这篇文章中,作者通过一个简单的似然函数优化(Maximum Likelihood Optimization)问题来对比 Julia,R 和 Python。这是一个比较小的优化问题,性能上的差异表现可能不太明显,但解决问题的过程能很好地反应三者各自的优劣势。

作者在撰写本文时,对这三种语言的熟悉程度如下:

语言实战经验
R9 年
Julia6 个月
Python新手

Julia 布道者 ChrisRackauckas 曾经说过:

如果你用 Julia 处理一个 10 秒内的问题,它的优势并不能体现出来。 而一旦处理的问题变复杂,需要花费比较长的时间,这时 Julia 的优势就会慢慢体现了。

有人用 Python 和 Julia 做过对比实验。以 10⁵ 为界点进行计算,当数值比 10⁵ 更小时 Python 比 Julia 快的。但数值大于 10⁵ 后,Julia 的速度就比 Python 快很多了。

优化问题

观察序列 Q1,Q2,...,Qn,我们需要找到优化该似然函数的参数 μ 和 σ:

通常我们会尝试优化对数似然:

在统计学上,这是截断的正态分布的最大似然估计(MLE)。

Julia 的测试情况

以下是作者使用 Julia 进行测试的情况。使用 Julia 中的 Optim.jl,可以直接使用特殊符号(symbols)作为变量名称,按照使用习惯,此处作者使用了希腊字母 μσ。Julia 还有一个 JuMP.jl 包用于优化问题。但 JuMP.jl 更适合用于更高级的优化问题,用在此处有点小题大做。

Julia 第一次优化

Julia 在执行第一次优化用了 7 秒,比 R 和 Python 都慢。对此,ChrisRackauckas 指出:

如果你需要解决 100 个 10 秒的优化问题,第一次执行需要花费 17 秒,接下来的优化不需要编译,大约只需要 10 秒。因此,总运行时常为 1007 秒。所以,当用 Julia 处理一个 10⁵ 秒的问题时,这 7 秒基本可以忽略不记;但如果用 Julia 处理 5 秒甚至更小的问题时,这 7 秒的差异就特别明显。

作者在下方硬编码了在 MLE 估计中使用的 Q_t 的值:

using Distributions, Optim# hard coded data\observationsodr=[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 parameterfunction 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)]) # 7.5 seconds@time res = optimize(neglikfn, [mean(Q_t), std(Q_t)]) # 0.000137 seconds# the \mu and \sigma estimatesOptim.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]

输出效果如下,排版看起来很舒服,也支持数学公示显示:

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

由此看出 Julia 的优势

指数描述
 Truncated(DN, lower, upper) 是定义截断分布的非常简单的方法
 logpdf 函数适用于任何分布式函数
 输出结果条理清晰,可读性强

Julia 的不足:

指数描述
 如果只是处理 10 秒内的简单问题,7.5 秒的编译时间会很烦人

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 的优势:

指数描述
 又处理截断正态的专用包
 马上输出结果,编译比 Julia 快

R 的不足:

指数描述
 截断正态没有对数密度; 没有简单的方法来定义任意分布的截断分布; 稀疏输出

Python 的测试情况

作者利用已有的 Python 学习经验想出如下方案,输入代码:

import numpy as npfrom scipy.optimize import minimizefrom scipy.stats import norm# generate the dataodr=[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 datadef 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 optimizeneglik = 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 的优势:

指数描述
 易于学习,各种支持非常好
 能很快输出结果,比 Julia 编译快

Python 的不足:

指数描述
 输出的可读性有待提高

综上所述,三种的综合对比如下:

语言优势不足
Julia

易于使用;完美支持截断正态分布;可读性强

第一次运行编译时间长
R易于使用可读性对比 Julia 较差
Python易于使用可读性对比 Julia 较差


全部评论(0条)
亲~快来评论噢!
Java帮帮公众号生态

Java帮帮公众号生态

总有一款适合你

Java帮帮-微信公众号

Java帮帮-微信公众号

将分享做到极致

Python帮帮-公众号

Python帮帮-公众号

人工智能,爬虫,学习教程

大数据驿站-微信公众号

大数据驿站-微信公众号

一起在数据中成长

九点编程-公众号

九点编程-公众号

深夜九点学编程

程序员服务区-公众号

程序员服务区-公众号

吃喝玩乐,听学吐画

Java帮帮学习群生态

Java帮帮学习群生态

总有一款能帮到你

Java学习群

Java学习群

与大牛一起交流

大数据学习群

大数据学习群

在数据中成长

九点编程学习群

九点编程学习群

深夜九点学编程

python学习群

python学习群

人工智能,爬虫

测试学习群

测试学习群

感受测试的魅力

Java帮帮生态承诺

Java帮帮生态承诺

一直坚守,不负重望

初心
勤俭
诚信
正义
分享
战略合作
关于我们
友链申请
友链交换:加帮主QQ2524138991 留言即可 24小时内答复  
快速换友链:在你的网站设置好Java帮帮-IT免费资源网友链,截图发给帮主即可
会员登录
获取验证码
登录
登录
我的资料
留言
回到顶部