BesselK with the order (v) as real (parameter)

Modeling
1

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?

Thanks

2

I think you have to wait until integral_1d is exposed.

3

Thanks Ben,

are you able to give a rough time estimate? (so I can plan some work for the next weeks/months).

1 Like
4

It is in math. I think we just need to merge the PR for the language.

5

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.

1 Like
6

The PR was closed

7

@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

functions {
  real log_besselk_frac(real v, real z);
}

in your stan file and compile the model with

includes = '\n#include "besselk.hpp" \n'
2 Likes
8

We might see if we can extend the approach given in

to the BesselK case.

1 Like
9

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.