CDF of a Von Mises distribution

General
1

Hello
I need to use the CDF of a Von Mises distribution (like scipy.stats.vonmises.cdf would do in Scipy). Stan does not seem to have a cdf for the Von Mises, only the lpdf. Am I wrong ? Is there a nice way to compute the cdf using the von_mises_lpdf already existing in stan (in order to avoid computing the infinite series).
Thank you

1 Like
2

hey i dont think Stan has a Von Mises CDF function but it should be straightforward to write your own by porting the scipy python function available here to a Stan user defined function.

2 Likes
3

indeed, thanks !

4

I opened a math library issue to request someone add it to our library so that we can have it in the language:

We also need the Cauchy circular distribution!

1 Like
5

this is my copy of the scipy one, I hope it’s right :

real von_mises_cdf_series(real k,real x,int p) {
		real s=sin(x);
		real c=cos(x);
		real sn=sin(p*x);
		real cn=cos(p*x);
		real R=0;
		real V=0;
		real new_R ;
		real new_V ;
		real new_sn ;
		real new_cn ;
		int n ;
		for (index in 1:p-1) {
			n=p-index;
			new_sn = sn*c-cn*s;
			new_cn = cn*c+sn*s;
			sn=new_sn;
			cn=new_cn;
			new_R = 1.0/(2.0*n/k+R);
			R=new_R;
			new_V=R*(sn/n+V);
			V=new_V;
		}
		return 0.5+x/(2.0*pi())+V/pi();
	}

	real von_mises_cdf_normalapprox(real k, real x) {
		real b=sqrt(2.0/pi())*exp(k)/modified_bessel_first_kind(0,k);
		real z=b*sin(x/2.0);
		return normal_cdf(z,0,1);
	}

	real von_mises_cdf(real k, real x) {
		real x2;
		real x3 ;
		int ck;
		real a[4] ;
		real F;
		real value_prov2 ;
		x3 = 2.0*pi()*round(x/(2.0*pi())) ;
		x2 = x-x3;
		ck=50;
		a[1] = 28.0;
		a[2] = 0.5;
		a[3] = 100.0;
		a[4]= 5.0;

		if (k<ck) {
			real value_prov = a[1]+k*a[2]-a[3]/(k+a[4]);
			real p=ceil(value_prov);
			int p2 = 1; 
			while (p2 < p) {p2+=1;}
			value_prov2 = von_mises_cdf_series(k,x2,p2);
			F=fmax(0, fmin(value_prov2,1)); // equivalent to clip in numpy
		}
		else {
			F = von_mises_cdf_normalapprox(k, x2);
		}
		return F+x3;
	}

	real von_mises_cdf_with_loc(real k, real x, real loc) {
		real val ;
		real addition_cdf ;
		addition_cdf = 0 ;
		val = x-loc ;
		while (val>pi()) {
			addition_cdf+=1;
			val-=2*pi();
		}
		while (val<-pi()) {
			val+=2*pi();
			addition_cdf-=1;
		}
		return von_mises_cdf(k,val) + addition_cdf ;
	}

[edit: escaped code]

2 Likes
6

It’s generally much easier to understand code if you (a) add some space around operators, (b) declare variables near where they’re used, and © eliminate unneeded temps. Here’s how I’d have coded that first function:

real von_mises_cdf_series(real k, real x, int p) {
  real s = sin(x);
  real c = cos(x);
  real sn = sin(p * x);
  real cn = cos(p * x);
  real R = 0;
  real V = 0;
  for (index in 1:p - 1) {
    int n = p - index;
    real new_sn = sn * c - cn * s;
    cn  = cn * c + sn * s;
    sn = new_sn;
    R = 1.0 / (2.0 * n / k + R);
    V = R * (sn / n + V);
  }
  return 0.5 + x / (2.0 * pi()) + V / pi();
}
2 Likes
7

thank you

8

Just a suggestion Bob, but it may be worth doing a diff between distributions available in scipy and those in Stan more generally and porting then across en masse.

It may help shortcut past whatever journey the scipy folks had to go on to decide they should add the distributions we don’t have.

9

also I got an error because round(0.5)=0 in python and round(0.5)=1 in stan. Just to keep in mind

3 Likes
10

Hello yes indeed would be great to have the CDF of the von mises in Stan. My solution works but is extremely slow…