search

alashworth / test-issue-import

0 stars 0 forks source link

Bessel function of fractional order #101

Open alashworth opened 5 years ago

alashworth commented 5 years ago

Issue by ppernot Monday Jun 27, 2016 at 15:11 GMT Originally opened as https://github.com/stan-dev/stan/issues/1939

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

alashworth commented 5 years ago

Comment by syclik Monday Jun 27, 2016 at 15:16 GMT

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 https://github.com/stan-dev/stan/issues/1939, or mute the thread https://github.com/notifications/unsubscribe/AAZ_F6TGvZwHW599wLXJM3XW_towNPwCks5qP-gSgaJpZM4I_MkO .

alashworth commented 5 years ago

Comment by ppernot Monday Jun 27, 2016 at 15:33 GMT

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

alashworth commented 5 years ago

Comment by bgoodri Tuesday Jun 28, 2016 at 02:16 GMT

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

alashworth commented 5 years ago

Comment by ppernot Tuesday Jun 28, 2016 at 06:48 GMT

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

alashworth commented 5 years ago

Comment by bgoodri Tuesday Jun 28, 2016 at 07:19 GMT

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 https://github.com/stan-dev/stan/issues/1939#issuecomment-228964883, or mute the thread https://github.com/notifications/unsubscribe/ADOrqqMiEUYlwgxL5f3d0CMOmAv6S_nDks5qQMOpgaJpZM4I_MkO .

alashworth commented 5 years ago

Comment by ppernot Tuesday Jun 28, 2016 at 07:30 GMT

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

dBesselI = function(x,nu) {
  besselI(x,nu-1) - nu/x*besselI(x,nu)
}
alashworth commented 5 years ago

Comment by bgoodri Tuesday Jun 28, 2016 at 07:40 GMT

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 https://github.com/stan-dev/stan/issues/1939#issuecomment-228973106, or mute the thread https://github.com/notifications/unsubscribe/ADOrqqHyZ23dfsUZ3DAx9CaTtPte45q4ks5qQM2WgaJpZM4I_MkO .

alashworth commented 5 years ago

Comment by ppernot Tuesday Jun 28, 2016 at 07:45 GMT

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

alashworth commented 5 years ago

Comment by bgoodri Tuesday Jun 28, 2016 at 07:46 GMT

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 https://github.com/stan-dev/stan/issues/1939#issuecomment-228976008, or mute the thread https://github.com/notifications/unsubscribe/ADOrqhkNL7-Bc-mDNHidNYBWzwlBHNdGks5qQNEGgaJpZM4I_MkO .

alashworth commented 5 years ago

Comment by ppernot Tuesday Jun 28, 2016 at 08:55 GMT

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 ;-)

alashworth commented 5 years ago

Comment by bgoodri Tuesday Jun 28, 2016 at 09:14 GMT

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 https://github.com/stan-dev/stan/issues/1939#issuecomment-228991354, or mute the thread https://github.com/notifications/unsubscribe/ADOrqm2TXhIvj6wyLuOhtX2DU4DLjMLPks5qQOGcgaJpZM4I_MkO .

alashworth commented 5 years ago

Comment by bob-carpenter Friday Mar 30, 2018 at 20:15 GMT

See discussion at: http://discourse.mc-stan.org/t/log-of-modified-bessel-first-kind/842/6