circular parameter

[Hugh]

Hi Stan users,

I would like to ask about estimating a parameter that for a circular quantity (e.g. time).

I have included a toy dataset/models in the attached file that runs in r.

I am trying to model some oscillating time series data - e.g. y[i,j] ~ normal(sin(time[I,[j]+pos[j]), sigma)

Different subjects have different phase shifts (pos[j]) in the range -pi to pi

I would like to estimate the phase shift parameter.

Initially I used a uniform prior over the range -pi to pi for offset (model 1 in attached file). This resulted in the posterior distributions of pos[j] near -pi or pi being poorly estimated, presumably as the optimization method has difficulties moving between -pi and pi.

Then I tried a von meis prior with a very small kappa (0 is not allowed) for offset (model 2 in attached file). The results were similar to model 1.

Then I tried using 2 parameters with uniform prior over the range -1 to 1 (model 3, xpos and ypos). I then took the atan of ypos/xpos to get the angle. Since there are multiple values of xpos and ypos that give the same angle, I set parameter sigma to be equal to sqrt(xpos^2+ypos^2)+0.1 , to make xpos and ypos identifiable. The seemed to work better than models 1 or 2, but had the disadvantage that each individual now had their own sigma.

By chance I noticed that if I neglected to index sigma in the likleyhood function (line 241, model 4), then the estimates for the angle were more precise . Whilst the print of the model indicated there were  still as many sigma as individuals, they were now all of the same value.

My questions are

- the way I used to estimate the offset/angle in model 3 or 4 seemed better than in model 1 or 2, but I could not find examples of others doing this, so I was unsure if this is a reasonable way to do it?

- What is going on in model 4 vs model 3? I would  have expected not indexing sigma to result in an error?

Thanks for any advice you may have

Hugh

Please see the attached modell/data  file.

[Bob Carpenter]

There's a large literature on circular statistics.

We'll finally have a working version of unit vector
parameters in Stan 2.10 (out soon). And I wrote up a new
chapter for the manual (warning---I had to learn as I
went, never having done circular stats). I'm attaching
it prior to release (comments welcome, as always).

Your diagnosis is accurate --- the sampler can't move
between -pi and pi if you treat it as a line. Stan uses
an auxiliary variable parameterization somewhat like yours
(but unconstrained), and a unit normal regulaization
described in the attached chapter; there aren't any free
parameters in that, though --- you just need to fix the
scale.

If you did this

for (n in 1:N)
y[n] ~ normal(mu[n], sigma);

and sigma is an N-vector, then it's equivalent to this

for (n in 1:N)
for (m in 1:N)
y[n] ~ normal(mu[n], sigma[m]);

which is not what you want to do. It artificially multiplies
data size by N and ties all the parameters sigma[n] together.

- Bob

> --
> 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.
> <circular example.txt>

[Hugh]

Thanks Bob - exciting to hear unit vector parameters are on the way.

I think you may have forgotten to attach the chapter - would be most interested to see this,

Best,

Hugh

[Bob Carpenter]

I indeed forgot to attach it last time. Here it is.

- Bob

[Hugh]

Thanks Bob,

I read the manual bit you sent which was most helpful. I changed the parameterization I was using based on what is in the manual, and used unit normal priors. It worked much better, but I got some divergent transitions (model 7 in the attached). I tried reducing the standard deviation of the unit normal priors and this reduced the number of divergent transitions (model 8 in attached). Reducing the SD further just made the posterior distributions more uniform (model 9).

Any thoughts on why I get the divergent transitions, and why there are less when I reduce the SD of the normal priors?

Best wishes

Hugh

[Bob Carpenter]

The program looks OK, but maybe Michael will have some
clue.

You could start by trying to put a prior on sigma.

If you're in R, there are ways to see what the
divergences are graphically (but I can never remember what
they are).

Some efficiency tips:

You can replace ((x^2)+(y^2))^0.5 with hypot(x, y).

I don't think you have to worry about

if (x_t==0&&y_t==1) {ans <- 0.5;} // to avoid div by 0
if (x_t==0&&y_t==-1) {ans <- 1;} // to avoid div by 0

unless you'r initializing at 0, because those are measure-0
points. I'd also recommend some space around operators and in between
curly braces, as well as consistent indentation. You can
replace this:

for(h in 1:H) {
xpos[h] ~ normal(0,1);
ypos[h] ~ normal(0,1);
}

with

xpos ~ normal(0, 1);
ypos ~ normal(0, 1);

It's much more efficient.

We'll soon have the unit_vector data structure built in and
also have sin() vectorized so the likelihood can be easily
vectorized.

- Bob

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

> <circular example_v2.txt>

[Hugh]

Hi Bob,

Thanks for your suggestions. I put a prior on sigma, and made the efficiency changes. Still getting divergent transitions, even if I increase adapt_delta to 0.99.

Below is the pairs plot with the first and second subjects parameters and sigma. I note the x and y variables are correlated as one would expect, since the ratio stays the same as they increase multiplicatively. Do you think that could be causing the sampling problems?

Best wishes

Hugh

[Bob Carpenter]

I retagged the subject and maybe Michael (who knows most about
the sampler) or Ben (who suggested this parameterization) can jump
in.

- Bob

> On Apr 24, 2016, at 10:22 PM, 'Hugh' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Hi Bob,
>
> Thanks for your suggestions. I put a prior on sigma, and made the efficiency changes. Still getting divergent transitions, even if I increase adapt_delta to 0.99.
> Below is the pairs plot with the first and second subjects parameters and sigma. I note the x and y variables are correlated as one would expect, since the ratio stays the same as they increase multiplicatively. Do you think that could be causing the sampling problems?
>
> Best wishes
>
> Hugh

> <Auto Generated Inline Image 1.png>
> <Auto Generated Inline Image 2.png>

> <r output.txt><circular example_v3.txt><Auto Generated Inline Image 1.png><Auto Generated Inline Image 2.png>

[Michael Betancourt]

I think you’re going to need to dump to a diagnostic file to get
at the unconstrained parameters to have any idea from where
the divergences arise.

[Bob Carpenter]

They're the bright blue dots in the pairs() plots.

They look like they're all over the posterior here, which
hasn't been the case with other examples I've seen, where
they tend to clump near boundaries.

- Bob

[Michael Betancourt]

Yeah, I know. I think the constraining transformation is
obfuscating the source of the divergences.

[Hugh]

Hi,

Thanks for having a look. I simplified the example so that all individuals are of the same phase and moved the transformation from the functions to  the transformed parameters block so as to be able to see all of the parameters on the pairs plot. Still the divergent transitions seem well spread. I attach  the model (circular_examplev5.txt), pairs plot,  diagnostic (mod12diag1->3) and sample  (mod12samp1->3) files and would be most grateful for any advice.

Best,

Hugh