Von_mises_rng with kappa=0

Modeling
1

Hi all,
I am trying to fit a mixture of a von Mises distribution and a uniform distribution on the circle to my data (responses from a memory experiment, where subjects select the colour of remembered items on a 360 degrees colour wheel).
I wanted to start out by generating some data from the model. Taking advantage of the fact that when kappa=0, the von Mises distribution is the uniform distribution on the circle (similarly to the implementation of my colleagues in JAGS, see Oberauer, Stoneking, Wabersich, & Lin, 2017, Journal of Vision), i wrote the following stan code:

data { 
  real<lower=0> p_mem; 
  real<lower=0> kappa;
  int n_sim; 
} 
model { } 
generated quantities { 
  real y[n_sim];
  real z[n_sim];
  for (n in 1:n_sim) {
    z[n] = bernoulli_rng(p_mem);
    y[n] = von_mises_rng(0, z[n]*kappa); 
  } 
}

and ran it using:

mixture_dat <- list(kappa = 20, p_mem = 0.7, n_sim = 1000) #p_mem is the mixing coefficient indicating the probability of belonging to the von Mises
gen_mixture_one_level <- stan(file = "generate_mixture_one_level.stan", data = mixture_dat,
                              iter=1, chains = 1, algorithm = "Fixed_param")

the error I got was:

[2] "  Exception thrown at line 12: von_mises_rng: inverse of variance is 0, but must be > 0!"

I also stumbled upon this potentially related discussion:

Any corrections or suggestions?
Thanks!

Dan

2

Judging from that issue, it’s a bug in our implementation.

First, you usually only want to do one simulation in the generated quantities and then just let Stan run multiple iterations. Otherwise, you have n_sim copies of the same variables here. So this should probably look like this once you include constraints:

generated quantities {
  int<lower = 0, upper = 1> z = bernoulli_rng(p_mem);
  real<lower = -pi(), upper= pi()> y = von_mises_rng(0, z * kappa);
}

You can code around the bug in von_mises_rng with:

y = z ? von_mises_rng(0, kappa) : uniform_rng(-pi(), pi());

But the problem you’re going to run into is that a real value running from -\pi to \pi isn’t the right topology, because those two values should be next to each other. It’s possible to use unit_vector to do proper circular statistics and I wrote some discussion up in the manual, but we don’t even have a simple case study for it.

1 Like
3

Thanks Bob.