I am testing calibration in R of a model in Julia. It works, but having to use juliaGet() for each model call makes it extremely slow. Is there a faster way to do it, where the Julia functions act like native R functions without the overhead?
# linmod_module.jl
module linmod_module
export linmod, linmodvec, loglikelihood
function linmod(x, a, b)
y = a + b * x
return y
end
function linmodvec(datax, a, b)
y = linmod(datax, a, b)
return y
end
function loglikelihood(datay, datax, a, b, error)
datan = size(datax, 1)
y = linmod(datax, a, b)
z = (datay - y) / error
ll = - 0.5 * datan * log(2 * pi * error * error) - 0.5 * sum(z .* z);
return ll
end
end
# try BayesianTools with Julia model
library(XRJulia)
library(BayesianTools)
# library(coda)
# make some data
x_data <- seq(0, 1, length.out=1000)
y_data <- 0.5+0.5*x_data+rnorm(length(x_data),0,0.1)
model_se <- 0.1
# find Julia functions
findJulia(test=TRUE)
ev <- RJulia() # create Julia evaluator
ev$AddToPath(getwd()) # find module linmod
# test Julia model
linmod <- JuliaFunction("linmod", module="linmod_module")
test <- linmod(0.5, 0.5, 0.5)
# test vectorised Julia model
linmodvec <- JuliaFunction("linmodvec", module="linmod_module")
test <- linmodvec(x_data, 0.5, 0.5)
juliaGet(test) # but this is very slow
# test Julia loglikelihood
loglikelihood <- JuliaFunction("loglikelihood", module="linmod_module")
test <- loglikelihood(y_data, x_data, 0.5, 0.5, model_se)
# BayesianTools setup
bt_names <- c("a", "b")
bt_prior <- createUniformPrior(c(0,0), c(1,1))
bt_likelihood <- function(params){
loglikelihood(y_data, x_data, params[1], params[2], model_se)
}
bt_parallel <- FALSE
bt_setup <- createBayesianSetup(likelihood=bt_likelihood,
prior=bt_prior,
parallel=bt_parallel,
names=bt_names)
nBurnin <- 900
nSample <- 900
nTotal <- nSample + nBurnin
nChains <- 3
nInternal <- 3
bt_settings <- list(startValue=nInternal, # DREAMzs internal chains
iterations=nTotal/nChains, # iterations per external chain
nrChains=nChains,
burnin=nBurnin/nChains+1, # discard these ones
parallel=bt_parallel,
message=TRUE)
bt_out <- runMCMC(bayesianSetup=bt_setup,
sampler = "DREAMzs",
settings=bt_settings)
suppressWarnings(summary(bt_out))
tracePlot(bt_out)
correlationPlot(bt_out)
bt_predict <- function(params){
juliaGet(linmodvec(x_data, params[1], params[2]))
}
bt_error <- function(mean, params){
return(rnorm(length(mean), mean=mean, sd=model_se))
}
bt_samples <- getSample(bt_out, start=2)
summary(bt_samples)
# this throws an error
suppressMessages({
plotTimeSeriesResults(sampler=bt_samples,
model=bt_predict,
observed=y_data,
error=bt_error,
plotResiduals=TRUE,
main="Residual Analysis")
})
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent