Limits on positive ordered

[Sachinthaka Abeywardana]

I am currently running rstan 2.8.2 and trying to put limits on a postive_ordered variable such that:

parameters{

positive_ordered[k]<lower=1/max(y), upper=1/min(y)> lambda;

}

However, it doesnt like the limits in <> braces. If I were to take them out it works perfectly. What is the correct way of implementing these limits for positive ordered variables?

Cheers,

Sachin

[Daniel Lee]

Hi Sachin,

1. The positive_ordered type doesn't have range constraints. If you look at the Stan reference guide, you can find it in the BNF grammars (in the Modeling Language Syntax section).

2. If it were to be allowed, it would be before the brackets: positive_ordered<lower=1/max(y), upper=1/min(y)>[k] (but this won't work)

3. There are lots of ways to accomplish this in Stan. One of the ways is to build the transform yourself. Another could be to use a simplex and each value could represent the difference between two lambdas.

Hopefully that helps.

Daniel

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

[Andrew Gelman]

Programming the transform oneself could be a lot of work and it would probably slow down the computation.  A simplex would work (if the simplex variable is x, then then ordered variable is just “lower” + (“upper”- “lower”) * cumulative sum of x, and I assume Stan has a cumulative sum function.  Or it’s possible that the upper and lower constraints should really be soft in any case (I’m thinking this because it seems unusual to define a variable with data-based constraints).

A

[Sachinthaka Abeywardana]

I think I should have mentioned that I am hoping to put exponential priors on each lambda such that:

lambda[1]~exponetial(1) 1(lambda[1]\in[a,b])

lambda[2]~exponetial(1) 1(lambda[2]\in[a,lambda[1]]) ...

lambda[k]~exponetial(1) 1(lambda[k]\in[a,lambda[k-1]])

where [a,b] are the higher and lower limits. 

So I'm not entirely sure where or how a simplex would fit into this framework. 

S

[Andrew Gelman]

You can throw priors on to any model in addition to whatever else is there.

[Bob Carpenter]

One way to use a simplex to generate
a positive ordered K-vector lambda that falls in
the interval (a, b) is as follows:

parameters {
simplex[K + 1] theta;
...
}
transformed parameters {
positive_ordered[K] lambda;
...

lambda <- a + head(cumulative_sum(theta), K) * (b - a);
...
}

You can then use a uniform prior on kappa and a Dirichlet(alpha)
prior on kappa, which is nice because it'll control spacing.
The prior alpha can be used to control how uniform the spacing will
be.

To see how the cumulative sum works, take K=3, so we
generate a (K+1) simplex such as

(0.2, 0.1, 0.3, 0.4)

Its cumulative sum is

(0.2, 0.3, 0.6, 1.0)

The first K elements are picked out by the head (0.2, 0.3, 0.6),
and then those are multiplied by (b - a), and then you add
a.

(a + (b-a) * 0.2, a + (b-a) * 0.3, a + (b - a) * 0.6)

It's all nice and symmetric and you don't need any Jacobians anywhere.

I have no idea how well it'll work in practice, though, so please
report back if you try it!

We don't have good prior recommendations for these cases yet, either.

- Bob

[Sachinthaka Abeywardana]

So I made the model as Bob suggested and looks like:

parameters {
 simplex[k] pi;
 simplex[k+1] theta;
}
transformed parameters{
    positive_ordered[k] lambda;
    lambda <- a + head(cumulative_sum(theta), k) * (b - a);
}
model {
 lambda~exponential(1);
        ....
}

(The entire model is attached in the files below).

So I have a few questions on this even though it works:

1. Is it valid to do lambda~exponential(1). The reason I ask this is because lambda depends on theta (a simplex). Also isn't the model going to try to put a uniform distribution on theta, even though I specified the prior on lambda? Basically how does stan consolidate the fact that it has an exponential prior and the restrictions placed by the transformed parameters block (a rough explanation would be fine).

2. I got the following warning and not sure what it means:

Warning (non-fatal):

Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.

If so, you need to call increment_log_prob() with the log absolute determinant of the Jacobian of the transform.

Left-hand-side of sampling statement:

    lambda ~ exponential(...)

3. Numerical problems:

validate transformed params: lambda is not a valid ordered vector. The element at 2 is 0.00089928, b     2

validate transformed params: lambda is not a valid ordered vector. The element at 2 is 3.72184, but      1

Does the above mean that there were 3 numerical problems. Apologies if I was not meant to ask this on the

[Bob Carpenter]

> On Dec 1, 2015, at 8:16 PM, Sachinthaka Abeywardana <sachin.ab...@gmail.com> wrote:
>
> So I made the model as Bob suggested and looks like:
> parameters {
> simplex[k] pi;
> simplex[k+1] theta;
> }
> transformed parameters{
> positive_ordered[k] lambda;
> lambda <- a + head(cumulative_sum(theta), k) * (b - a);
> }
> model {
> lambda~exponential(1);
> ....
> }
>
> (The entire model is attached in the files below).
>
> So I have a few questions on this even though it works:

I'm glad it works.

> 1. Is it valid to do lambda~exponential(1). The reason I ask this is because lambda depends on theta (a simplex). Also isn't the model going to try to put a uniform distribution on theta, even though I specified the prior on lambda?

What matters if the final log density. Theta is
involved in both a default uniform prior over the simplex,
plus the distribution on lambda, which is defined in terms
of theta. Changing theta affects the value of lambda, which
in turn affects the exponential distribution's contribution
to the log density.

> Basically how does stan consolidate the fact that it has an exponential prior and the restrictions placed by the transformed parameters block (a rough explanation would be fine).

The way it really works is through HMC. There's an algorithm
section in the manual.

> 2. I got the following warning and not sure what it means:
> Warning (non-fatal):
> Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
> If so, you need to call increment_log_prob() with the log absolute determinant of the Jacobian of the transform.
> Left-hand-side of sampling statement:
> lambda ~ exponential(...)

This you have to worry about. Hence the warning.

Because you apply a change of variables, you have to apply
a Jacobian adjustment. That's also described thoroughly in
the manual.

In this case, you're lucky, because a and b are constants
and the cumulative sum has a unit Jacobian determinant, so it's
a constant linear transform of a unit Jacobian, so there's no
affect on the posterior from the change of variables.

So you don't need to do anything.

I still think the Dirichlet on the bins and maybe an exponential
on the overall scale makes a lot of sense here, especially if
these are cutpoints for an ordinal regression.

> 3. Numerical problems:
> validate transformed params: lambda is not a valid ordered vector. The element at 2 is 0.00089928, b 2
> validate transformed params: lambda is not a valid ordered vector. The element at 2 is 3.72184, but 1
> Does the above mean that there were 3 numerical problems. Apologies if I was not meant to ask this on the

Looks like your response got cut off here.

How often does this turn up? I don't see how lambda can fail
to be an ordered vector. Is a < 0 perhaps?

- Bob

> <model2.stan><lambdaStan1.R><y250.csv>