Cauchy equivalent of Von Mises distribution in Stan

[Ben Lambert]

Hello all,

I am working on a problem where we track animal paths over time. In particular we track animal's angle of motion relative to their previous path, and want to build a model for angle of motion. Our initial thoughts were to use the Von Mises distribution. 

However, the distribution of angles moved by the animals has much fatter tails than the Von Mises (itself basically a circular normal.) Data that looks similar to our real data can be generated via the following in R:

library(CircStats)

phi <- rwrpcauchy(1000000,pi,rho = 0.7)-pi

hist(phi,100)

Basically what we would ideally like to do is to use a more robust version of the Von Mises; equivalent to replacing the normal by a Student t or Cauchy in non-circular statistics. There is a circular version of the distribution I describe: it is called the wrapped Cauchy distribution (see https://en.wikipedia.org/wiki/Wrapped_Cauchy_distribution) - the same distribution I used to generate the above.

My question is - how can I fit a similar distribution in Stan? I was thinking that basically I could do something like:

phi ~ von_mises(0,kappa);

kappa ~ gamma(a,a);

So basically use a mixture of Von Mises to yield something like a wrapped Cauchy. (This is through analogy of representing a Student T as a mixture or normals).

I have tried the above however, and it doesn't seem to work; the posterior predictive distribution for phi looks pretty much the same as for the vanilla Von Mises.

Does anyone have an idea as to how I might be able to do this type of workaround here?

Best,

Ben

PS Here is the full code:

data{

  int N;

  real phi[N];

}

parameters {

  real<lower=0> kappa;

  real<lower=0> a;

}

model{

  if (kappa < 100){

    phi ~ von_mises(0,kappa);

  } else{

    phi ~ normal(0,sqrt(1/kappa));

  }

  

  kappa ~ gamma(a,a);

  

}

generated quantities{

  real phi_gen;

  if (kappa < 100){

    phi_gen <- von_mises_rng(0,kappa);

  } else{

    phi_gen <- normal_rng(0,sqrt(1/kappa));

  }

  

}

[Mike Lawrence]

You'll want kappa to be a vector with same length as phi

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

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

~ Certainty is (possibly) folly ~

[Mike Lawrence]

You'll also probably want to put an explicit prior on "a". Right now it's uniform bounded at zero and infinity, which can be trouble.

[Ben Lambert]

Thanks very much! Really appreciated. In case anyone else needs it, the following worked perfectly for us!

data{

  int N;

  real phi[N];

}

parameters {

  real<lower=0> kappa[N];

  real<lower=0> a;

  real<lower=0> b;

}

model{

    for (i in 1:N){

      if (kappa[i] < 100){

        phi[i] ~ von_mises(0,kappa[i]);

      } else{

        phi[i] ~ normal(0,sqrt(1/kappa[i]));

      }

        kappa[i] ~ gamma(a,b);

    }

    a ~ gamma(2,1);

    b ~ gamma(2,1);

}