# Von_mises_rng with kappa=0

#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.

#3

Thanks Bob.