【发布时间】：2018-06-08 01:12:11
【问题描述】：
我用 R 和 Julia 语言构建了两个几乎相同的程序。
我事先知道有一种方法可以提高 Julia 代码的性能，因为我没有声明类型，而且通常非向量化代码往往更有效。但是，Julia 代码和 R 代码位于同一地形，因此可以进行比较。
注意：我开始学习 Julia 一周，所以你可以看出我对 Julia 的编程不是很好。
我注意到 Optim 包提供了optimize() 函数。与与 R 一起默认安装的 stats 包的 R 语言的 optim() 函数相比，它非常慢语言。
以下分别是代码R和Julia：
R 代码：
rm(list=ls(all=TRUE))
Gexp <- function(par,x){
lambda <- par[1]
pexp(q = x, rate = lambda, lower.tail = TRUE, log.p = FALSE)
}
gexp <- function(par,x){
lambda <- par[1]
dexp(x = x, rate = lambda, log = FALSE)
}
QGexp <- function(p,...){
qexp(p,...)
}
Gweibull <- function(par,x){
alpha <- par[1]
beta <- par[2]
pweibull(q = x, shape = alpha, scale = beta, lower.tail = TRUE, log.p = FALSE)
}
QGweibull <- function(p,...){
qweibull(p,...)
}
# Função de distribuição acumulada Exponentiated Kw-G class (EKw-G)
cdf_ekwg <- function(cdf,par,x,...){
a <- par[1]
b <- par[2]
c <- par[3]
npar <- c(...)
(1 - (1 - cdf(par = npar, x = x)^a)^b)^c
}
# cdf_ekwg(cdf = G, par = c(0.2,0.4,0.21), x = 1, alpha = 1.1, beta = 1.2, lambda = 1)
# Função densidade de probabilidade Exponentiated Kw-G class (EKw-G)
pdf_ekwg <- function(cdf, pdf, par, x, ...){
a <- par[1]
b <- par[2]
c <- par[3]
#cdf_ekwg_locale <- function(x){
# cdf_ekwg(cdf = cdf, par = par, x, ...)
#}
npar <- c(...)
g = pdf(par = npar, x = x)
G = cdf(par = npar, x = x)
a * b * c * g * G^(a-1) * (1-G^a)^(b-1) * (1 - (1-G^a)^b)^(c-1)
#numDeriv::grad(func = cdf_ekwg_locale, x = x, method = "simple")
}
#integrate(f = pdf_ekwg, par = c(1,1,1.5), lower = 0, upper = Inf, cdf = Gexp, pdf = gexp,
# lambda = 1.5)
# Será fixado os parâmetros de G. Serão estimados os parâmetros a, b e c do
# modelo EKwG.
sample_ekwg <- function(QG, n, par, ...){
a <- par[1]
b <- par[2]
c <- par[3]
u <- runif(n = n, min = 0, max = 1)
p <- (1 - (1 - u^(1/c))^(1/b))^(1/a)
QG(p = p, ...)
}
# Função de distribuição acumulada Exponentiated Kw-G class (EKw-G)
cdf_ekwg <- function(cdf,par,x,...){
a <- par[1]
b <- par[2]
c <- par[3]
npar <- c(...)
(1 - (1 - cdf(par = npar, x = x)^a)^b)^c
}
# cdf_ekwg(cdf = G, par = c(0.2,0.4,0.21), x = 1, alpha = 1.1, beta = 1.2, lambda = 1)
# Função densidade de probabilidade Exponentiated Kw-G class (EKw-G)
pdf_ekwg <- function(cdf, pdf, par, x, ...){
a <- par[1]
b <- par[2]
c <- par[3]
#cdf_ekwg_locale <- function(x){
# cdf_ekwg(cdf = cdf, par = par, x, ...)
#}
npar <- c(...)
g = pdf(par = npar, x = x)
G = cdf(par = npar, x = x)
a * b * c * g * G^(a-1) * (1-G^a)^(b-1) * (1 - (1-G^a)^b)^(c-1)
#numDeriv::grad(func = cdf_ekwg_locale, x = x, method = "simple")
}
# integrate(f = pdf_ekwg, par = c(1.4,1.3,0.5), lower = 0, upper = Inf, cdf = G, alpha = 1.1, beta = 1.2)
loglikelihood <- function(cdf, pdf, par, x, ...){
-sum(log(pdf_ekwg(cdf = cdf, pdf = pdf, par = par, x = x, ...)))
}
myoptim <- function(...) tryCatch(optim(...), error = function(e) NA)
G = Gexp
g = gexp
data = sample_ekwg(QG = QGexp, n = 550, par = c(1,1,1), rate = 1.5)
starts = c(1,1,1)
set.seed(0)
start = Sys.time()
for(i in 1:5){
result <- myoptim(par = starts, fn = loglikelihood, x = data, cdf = G,
pdf = g, method = "Nelder-Mead",rate = 1.5)
}
Sys.time() - start
Julia 代码：
using Distributions
#using Cubature # Calculo de integrais numéricas.
#using Plots
using Optim
#using JuMP
#using NLopt
function gexp(x,par)
λ = par[1]
λ * exp(-λ * x)
end
# valor = hquadrature(x -> gexp(x,1), 0, 100)[1]
function Gexp(x,par)
λ = par[1]
1- exp(-λ * x)
end
function QGexp(x,par)
λ = par[1]
# A função Exponential no pacote Distributions é reparametrizada
# como 1/lambda. Dessa forma, para trabalhar com densidade na forma
# λ * exp(-λ*x) é preciso tomar 1/λ.
quantile.(Exponential(1/λ),x)
end
function sample_ekwg(QG, n, par0, par1...)
a = par0[1]
b = par0[2]
c = par0[3]
u = rand(n)
p = (1 - (1 - u.^(1/c)).^(1/b)).^(1/a)
QG(p, par1...)
end
# Função de distribuição acumulada Exponentiated Kw-G class (EKw-G)
function cdf_ekwg(cdf, x, par0, par1...)
a = par0[1]
b = par0[2]
c = par0[3]
(1 - (1 - cdf.(x,par1...).^a).^b).^c
end
# Função densidade de probabilidade Exponentiated Kw-G class (EKw-G)
function pdf_ekwg(cdf, pdf, x, par0, par1...)
a = par0[1]
b = par0[2]
c = par0[3]
g = pdf(x, par1...)
G = cdf(x, par1...)
a * b * c * g * G.^(a-1) * (1-G.^a).^(b-1) * (1 - (1-G.^a).^b).^(c-1)
end
# valor = hquadrature(x ->
# pdf_ekwg(Gexp,gexp, x, [1,1,1], 1), 0, 100)[1]
function loglike(cdf, pdf, x, par0, par1...)
n = length(x)
soma = 0
for i = 1:n
soma += log(pdf_ekwg(cdf, pdf, x[i], par0, par1...))
end
return -soma # Queremos minimizar loglike.
end
G = Gexp
g = gexp
data = sample_ekwg(QGexp, 550, [1,1,1],1.5)
starts = [1,1,1]
par0 = [1,1,1]
par1 = [1.5]
srand(0) # set seed.
@time for i = 1:5
optimize(par0 -> loglike(G, g, data, par0, par1), [1.3,1.2,2.1])
end
执行这两个代码的机器配置如下：
注意：这些代码不适用于相同的样本。但是，我相信这不能证明计算时间的巨大差异是合理的。
在我的硬件中，R 代码花费 0.6508956 秒，Julia 代码花费 27.180257 秒。
有人知道如何让 Julia 的代码比 R 代码运行得更快吗？ 我想要一个简单的解决方案，因为它承诺在 Julia 中具有出色的计算性能，而无需对编程有太多了解。看，在 R 中，并没有做太多的事情来证明 Julia 代码中的主要修复是合理的。
最好的问候。
【问题讨论】：
-
R 代码是开源的。大概（或者我应该假设）Julia 语言有一个到 C 的接口？
-
别吐这么多。这对数组的性能真的很糟糕，因为数组的长度不是编译时间信息。只需传递数组就可以了。
-
慢的不是 Optim.jl，而是你的函数。您的
loglike被调用了大约 8500 次，优化在我的计算机上需要 16 秒，这与运行loglike8500 次所需的持续时间相同。 -
谢谢。我会尝试将考虑的迭代次数放在R的优化算法中。
-
使用 BenchmarkTools 中的一些宏来获得更高的精度，而不是
@time。