Re: [stan-users] Von Mises from new bessel function

[Michael Betancourt]

Not quite.  The problem is that there's also a Jacobian floating around from the transformation from x to cos(x - mu).  In fact, there are an infinite number of Jacobians because the transformation has to be evaluated at every zero, in this case cos(x - mu + 2pi), so you have to increment lp an infinite number of times.  This is actually a known problem in the circular statistics literature, and it arises because we sample from a non-compact space (R^{n}) where as the von Mises lives on a compact space (S^{1}).  Any transformation between the two topologies causes problems that either end up requiring things like complex numbers.

Fortunately Ben and Marcus realized that you could sample from S^{1} uniformly by using a trick that involves sampling two gaussians (actually, it can be two of any variate defined on all of R) and then normalizing.  You should be able to just use the code they were developing for the unit_vector type (did that ever get implemented?) in 2D and augment with a the Radon-Nikodym derivative between the von-Mises and uniform-on-a-circle distribution.  Which in this case I think just reduces to the density you quoted above.

On Thu, Sep 26, 2013 at 5:31 PM, Mike Lawrence <mike....@gmail.com> wrote:

Hi folks,

In anticipation of stan 2.0, I've been exploring the git repo. I *think* I've figured out how to implement the Von Mises density function using the new modified_bessel_first_kind function, but would appreciate input from others if not. 

So, to start, I had to figure out what the Von Mises density function is in the first place. The source for the "circular" package indicates that the density function is (in R code):

    dens <- 1 / ( 2 * pi * besselI( x = kappa, nu = 0, expon.scaled = TRUE ) ) * ( exp( cos( x - mu ) - 1 ) ) ^ kappa

where x and mu are in radians. Now, the above uses the expon.scaled=TRUE argument to bessell(), which is absent from the stan modified_bessel_first_kind() function, but I think that (again, in R):

    besselI(x=kappa,nu=0,expon.scaled=T) 

is equal to:

    besselI(x=kappa,nu=0,expon.scaled=F) * exp(-kappa)

(excepting a small precision difference). So, in my stan code I should be able to do

    lp__ <- lp__+ log( 1 / ( 2 * pi() * modified_bessel_first_kind(0,kappa) * exp(-kappa) ) * ( exp( cos( x - mu ) - 1 ) ) ^ kappa )

to achieve an increment of the lp__ by the log-density of the Von Mises at datum x, correct?

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Marcus Brubaker]

The unit_vector type hasn't been revised yet based on our discussions.  I'm actually working on a paper right now which tests out exactly the best way to do it in terms of sampling efficiency.  Once that's done and submitted I'll go back and fix the unit_vector parameterization.

Cheers,

Marcus

[Mike Lawrence]