I would like to use an implementation of BesselK that accepts the order (v) as real (parameter).
I found this thread Vignette for interfacing with external C++ code via rstan, that includes exactly the BesselK example.
I was just wandering if from 2016 there was an alternative method to include it in rstan (from R C++ code for example)
(The goal is to enable coding of distributions such as SICHEL https://www.rdocumentation.org/packages/gamlss.dist/versions/5.1-1/topics/SICHEL)
That thread also put forward concerns about it’s numerical stability for differentiation, any reccomendations?
I think you have to wait until
integral_1d is exposed.
are you able to give a rough time estimate? (so I can plan some work for the next weeks/months).
It is in math. I think we just need to merge the PR for the language.
Any update on that topic? I am trying to build up spatial regression model with Matérn covariance function that require modified Bessel of second kind with order as real parameter.
The PR was closed
@lionel68 The PR was closed for now, but there is hope :-)
It is currently out of my short-term plan to complete this, but if there is need and will to cooperate on a final push to make it happen, then a lot of the work has been done. If you read the comments on the PR (start at the end) you’ll notice that there is a remaining part of the parameter space where the code has large error. This needs to be resolved if the code is to be included in the main codebase and I ran out of ideas. The computation is also very much slow.
If you don’t care about this part of the parameter space I have adapted, a slightly older version of the code to work as a custom C++ function in Stan, so you can use it right away: https://github.com/stemangiola/RNAseq-noise-model/blob/SICHEL-reparametrised/dev/besselk.hpp (attribution would be nice)
To use it, have
real log_besselk_frac(real v, real z);
in your stan file and compile the model with
includes = '\n#include "besselk.hpp" \n'
We might see if we can extend the approach given in
to the BesselK case.
Thanks @martinmodrak for the update and the infos on where things are stuck. However given my very limited math and C++ skills I do not feel comfortable to make the push. Also it seems that some Matern covariance functions are part of the new Stan release: https://github.com/stan-dev/stan/pull/2758, I’ll wait until this trickles down to RStan to have a go at using it for the models I have in mind.