Bessel function of fractional order

[ppernot]

Summary:

Extension of Bessel and modified Bessel functions to fractional order

Description:

Extension to fractional order would conform with their implementation in base R.

Current Version:

v2.9.0

[syclik]

Do you have a specific use-case for fractional order Bessel functions and
modified Bessel functions?

On Mon, Jun 27, 2016 at 11:11 AM, Pascal PERNOT notifications@github.com
wrote:

Summary:

Extension of Bessel and modified Bessel functions to fractional order
Description:

Extension to fractional order would conform with their implementation in
base R.
Current Version:

v2.9.0

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
#1939, or mute the thread
https://github.com/notifications/unsubscribe/AAZ_F6TGvZwHW599wLXJM3XW_towNPwCks5qP-gSgaJpZM4I_MkO
.

[ppernot]

Hi,

yes, I am using them to compute second virial coefficients.

Here is the R implementation of my function.

virial2 = function(Temp,eps,sig) {
  avogad=6.0221353e23
  Tstar = Temp/eps
  x = 1/(2*Tstar)
  Bstar = 2^0.5*pi*x*exp(x)*(besselI(x,-3/4)+ besselI(x, 3/4)-besselI(x, 1/4)-besselI(x,-1/4))
  virial = 2*pi*sig^3*Bstar/3
  # Experimental results given in m^3/kmol
  virial = virial*(1000.*avogad)*1e-30
  return(virial)
}

Extension to stan would enable me to use it for bayesian inference on the parameters eps and sig...

Thanks,

Pascal

[bgoodri]

I have a unpushed branch that does this for the Bessel K function only.

[ppernot]

I am not familiar with stan development. Would it be difficult to extend
this to the I type ?

[bgoodri]

The particularity of Stan that we are dealing with here is that if we allow
the order to be a real number, then we have to be able to differentiate the
Bessel function(s) with respect to the order parameter and do so with a
high degree of numerical accuracy throughout the (nu, x) plane. The Bessel
K function is where I started because it is pretty well behaved numerically
since it is almost like an exponential PDF. The Bessel I function would be
harder to deal with, although it would help if you had any links to papers
that talk about differentiating it.

On Tue, Jun 28, 2016 at 2:48 AM, Pascal PERNOT notifications@github.com
wrote:

I am not familiar with stan development. Would it be difficult to extend
this to the I type ?

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#1939 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqqMiEUYlwgxL5f3d0CMOmAv6S_nDks5qQMOpgaJpZM4I_MkO
.

[ppernot]

Isn't a recurrence relation accurate enough ?
Ref. Wolfram

dBesselI = function(x,nu) {
  besselI(x,nu-1) - nu/x*besselI(x,nu)
}

[bgoodri]

Differentiating with respect to x is not so bad. Differentiating with
respect to nu is harder. We would need to implement (with numerical
accuracy) things like

http://functions.wolfram.com/Bessel-TypeFunctions/BesselI/20/01/01/

On Tue, Jun 28, 2016 at 3:30 AM, Pascal PERNOT notifications@github.com
wrote:

Isn't a recurrence relation accurate enough ?

dBesselI = function(x,nu) {
besselI(x,nu-1) - nu/x*besselI(x,nu)
}

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#1939 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqqHyZ23dfsUZ3DAx9CaTtPte45q4ks5qQM2WgaJpZM4I_MkO
.

[ppernot]

OK, I see the problem...
I assume we cannot skip differentiation wrt nu ?

[bgoodri]

We skip it now but only because we restrict nu to be an integer.

On Tue, Jun 28, 2016 at 3:45 AM, Pascal PERNOT notifications@github.com
wrote:

OK, I see the problem...
I assume we cannot skip differentiation wrt nu ?

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#1939 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqhkNL7-Bc-mDNHidNYBWzwlBHNdGks5qQNEGgaJpZM4I_MkO
.

[ppernot]

That would still hold if we constrain nu to be a multiple of 1/4...
An ugly solution, I know, but it would solve MY problem ;-)

[bgoodri]

Yeah, I would just define a function in the Stan language

real my_bessel(int nu4, real x) {
real nu;
nu <- 0.25 * nu4;
...
}

and implement whatever special case formula you need in your situation.

On Tue, Jun 28, 2016 at 4:55 AM, Pascal PERNOT notifications@github.com
wrote:

That would still hold if we constrain nu to be a multiple of 1/4...
An ugly solution, I know, but it would solve MY problem ;-)

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#1939 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqm2TXhIvj6wyLuOhtX2DU4DLjMLPks5qQOGcgaJpZM4I_MkO
.