Hi!

I am trying to program up the model for cylindrical variables from the paper by Mardia and Sutton which can be found here mardia_sutton_cylindrical_dist.pdf

This is a combination of values on a cylinder, of which one is cyclic and described by a von Mises distribution and the other is described by a normal distribution.

My final goal is to cluster/classify my data points, I have programmed it up in a similar way for only Gaussians.

Unfortunately I can’t write formula in here, but basically I want to program the 4 equations which can be found in (1.1) (1.2) and (1.3).

I generate the 4 parameters, mu and kappa for the von Mises distribution and mu and sigma for the normal distribution in the parameters block.

In the transformed parameters block I apply the equations in (1.3) and in the model block I apply the equation (1.2) before I finally use the equation (1.1).

Now I get the error even though afterward the warmup and sampling starts:

```
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can't start sampling from this initial value.
```

I don’t know where this comes from, maybe someone knows an answer?

Also do you see anything wrong with my approach to this code for the model?

Here is my code:

```
data {
int K; // number classes
int N; // number of all data points
vector[N] angle; //theta
vector[N] dist; //x
}
parameters{
// variables for von mises distribution, mu and kappa
vector [K] a_mu;
vector<lower=0, upper=90>[K] kappa;
// variables for normal distribution, mu and sigma
vector [K] d_mu;
vector<lower=0, upper=6> [K] d_sigma;
simplex [K] weights;
}
transformed parameters{
real p_1 = 0.5;
real p_2 = 0.5;
real p = sqrt(pow(p_1,2) + pow(p_2,2));
vector[K] d_sigma_c;
for (k in 1:K){
d_sigma_c[k] = d_sigma[k] * (1 - pow(p,2));
}
model{
for (n in 1:N) {
vector [K] pb;
vector[K] d_mu_c;
for (k in 1:K) {
// calculate d_mu_c with the helper function
d_mu_c[k] = d_mu[k] + d_sigma[k] * sqrt(kappa[k])
* (p_1 *(cos(angle[n]) - cos(a_mu[k])) + p_2*(sin(angle[n]) - sin(a_mu[k])));
//calculate real function
pb[k] = log(weights[k]) + (1/(2 * pi() * bessel_first_kind(0,kappa[k]))
* exp(kappa[k] * (cos(angle[n] - a_mu[k]))) * 1/sqrt(2 * pi() * d_sigma_c[k])
* exp(pow((dist[n] - d_mu_c[k]),2) / (2 * d_sigma_c[k])));
}
target += log_sum_exp (pb);
}
}
```