employing a ragged array of cutpoints

[terrance savitsky]

I have an multvariate ordered outcome, y[N,D], where the number of categories depends on dimension d = 1, ..., D.   I'd like to set an ordered_logistic likelihood for each y[i,d], where the length of each row of the cutpoint "matrix", gamma[d,] depends on "d".   I know this is not literally possible, but i'm hoping there is a way to do so without brute force definition of distinct vector objects for each dimension, d, with no dimension index tie between them.   What would make this go easily would be if I could define matrix[D,max(ncuts)] gamma and then use gamma[d,1:ncuts[d]] in the ordered_logistic statement, but I know we may not use partial indexing.

Thanks for your help.

--
Thank you, Terrance Savitsky

[terrance savitsky]

i'd really appreciate any advice on the forwarded topic, below, about dealing with a ragged array in stan.  in my case, it is an array of cutpoints where the number of "rows" depends on the the column index.  

[Ben Goodrich]

Is it so bad to make all variables have the maximum number of cutpoints? I don't know for a fact that Stan would still work, but I hope it would just make the posterior distribution of the extra cutpoints for each variable be the same as the prior distribution.

Ben

[terrance savitsky]

the problem is in use of the ordered_logistic likelihood.  the probability masses assigned to each observed category are keyed off of the length of the cutpoint vector (which we know is the number of categories minus one) so this isn't a case where a square matrix may be used in a fashion that allows the data to ignore some positions (and all will be well under a proper prior).  if is the max number of cutpoints is 5 over all dimensions and dimension "d" has 3 categories, then I need something like, y[i,d] ~ ordered_logistic(c[i,d],cutpoints[d,1:2]), where "cutpoints" is a D x 5 matrix.

--
 
 

[Ben Goodrich]

Maybe the attached patch will make it work? If the patch doesn't apply for whatever reason, just change

STAN_HOME/src/stan/prob/distributions/univariate/discrete/ordered_logistic.hpp

so that the line that says

int K = c.size() + 1;

becomes

int K = fmax(y);

Ben

[Ben Goodrich]

Actually, you would need to use max rather than fmax. Correct patch attached.

Ben

[Matt Hoffman]

I'd personally love to have ragged arrays in Stan. I'm looking at
modeling large numbers of sequences, most of which are very short but
some of which are longish.

Matt

> --
>
>

[Bob Carpenter]

It would be nice. Does R let you specify ragged
arrays some way?

There's always the explicit DB-like coding of a ragged
array using the indices. Suppose x is a ragged array
of size I, where the entry for index i in 1:I is
an array of size N[i].

K | ii nn word
------------------
1 | 1 1 17
2 | 1 2 192
3 | 1 3 42
...
92| 1 92 598
93| 2 1 17
94| 2 2 3
95| 2 3 971

The matrix x is determined by:

x[ ii[k], nn[k] ] = word[k]

If we assume the ii and nn are sorted as above, we can
use an index structure to indicate where the entries
for x[i,] start:

I| kk
-------
1| 1
2| 93
...

So kk[i] is the entry where the entries for x[i,] start
and size(x[i,]) in the ragged array is (kk[i+1] - kk[i])
(assuming you pad one on the end).

With this secondary kk structure, you can do away with the
i and n indices and just keep the word array.

For a start, it'd be easy to write a Stan operation

ragged_array_of_arrays(word,k)

that would take those word and start arrays and create a
ragged array of arrays. The same kind of thing would work
for a vector.

But there's no good type for the result in Stan. We
could always add something like:

int x[I,];

for a ragged array of arrays of ints.

I also could imagine specifying a ragged array in Stan
by providing its indexes.

int<lower=0> N; // size of ragged array
int<lower=0> x_sizes[N]; // size of each entry of ragged array
ragged_array x[x_sizes]; // ragged array

That's more than a bit clunky, but it would be
workable given what we have now.

Essentially x_sizes[N] is just the diff of kk with itself
offset by one.

- Bob

[W Robert Long]

R has stack() and unstack() which I believe can be used for ragged arrays

[Ben Goodrich]

On Monday, October 29, 2012 3:55:28 PM UTC-4, Bob Carpenter wrote:

It would be nice.  Does R let you specify ragged
arrays some way?

I think the most common way that an R user would approach this problem would be to have a list of arrays, like

my_list <- list(A = array_1, B = array_2)

For a start, it'd be easy to write a Stan operation 

   ragged_array_of_arrays(word,k)

that would take those word and start arrays and create a
ragged array of arrays.  The same kind of thing would work
for a vector.

But there's no good type for the result in Stan.

Would it make sense to utilize an STL container that accepts keys? Another alternative that would have been acceptable in Terrance's case would be to create an Eigen::Array whose "scalar type" is a vector_v or matrix_v or whatever. I think for dynamic sized matrices, they wouldn't have to all be the same size. I don't know if either of these would be less clunky than what you suggested below.

 

 We
could always add something like:

   int  x[I,];

for a ragged array of arrays of ints.

I also could imagine specifying a ragged array in Stan
by providing its indexes.

   int<lower=0> N;            // size of ragged array
   int<lower=0> x_sizes[N];   // size of each entry of ragged array
   ragged_array x[x_sizes];   // ragged array

That's more than a bit clunky, but it would be
workable given what we have now.

Essentially x_sizes[N] is just the diff of kk with itself
offset by one.

 Ben

[terrance savitsky]

Ben, You're likely not going to remember this issue where I have a matrix of ordered scores, y[i,d], where the number of categories depends on "d".   the ordered_logistic() distribution cannot address such a case.  You suggested a patch to the templated function ordered_logistic.hpp of the form:

int K = c.size() + 1;

becomes

int K = max(y);

This won't work because ordered_logistic is univariate (so that 'y' is an integer).   Setting int K = y won't work because y may not hold the maximum possible value for y[,d].   The immediate approach I think of, at present, is to edit my copy of ordered_logistic.hpp to take in a (new) integer indicator for "K" (set to max(y[,d]).  Before I made the change on my end, I would appreciate to have your input about whether you think such is okay.

terrance. 

[Ben Goodrich]

Argh. That might get vectorized soon, but apparently not yet. Passing in K is fine in principle, but you would have to doctor the code in several places to do it. Honestly, it is probably easier for you to roll your own ordered logistic distribution by manipulating lp__ in the .stan file.

Ben

[terrance savitsky]

thx, Ben. I would need to use int_step() to deal with multiple possible lp__ update nuggets based on whether 'y' == 1 or 'y==K' or a value in between.  I assume that's fine to do as a hack to implement a conditional statement.   i suppose it is, however, inefficient each of the 3 alternatives multiplied by an int_step() command are computed.

As an aside, I did doctor the ordered_logistic.hpp file in several places.  it was obvious where to do so as the code was well laid out.  that didn't work as I received the following error:

no matches for function name="ordered_logistic_log"

    arg 0 type=int

    arg 1 type=real

    arg 2 type=vector

    arg 3 type=int

unknown distribution=ordered_logistic

Parser expecting: <statement>

The "arg 3" input is my addition.   I thought the .hpp performed the prototyping/declaration??

terrance.

[Daniel Lee]

Can you clarify what you're trying to do with ordered_logistic? The distribution function expects 3 arguments, not 4.

--
 
 

[terrance savitsky]

absolutely correct.  i performed a direct edit of my copy of the ordered_logistic.hpp file (containing the templated function) to add another input for the number of categories so that the function could allow the number of categories to depend on dimension of the data.   

[Daniel Lee]

That gets a little tricky.

1) You will need to update all 4 definitions of the distribution function. I'm assuming you've done that successfully. If not, the policy type will match the fourth argument and the behavior won't be as expected.

2) You need to update stanc. The place you'll need to touch is src/stan/gm/function_signatures.hpp. You'll need to update the function signature accepted by the parser for ordered_logistic_log.

Once all that's done, you should be able to get away with it.

If you're just altering generated code, you won't need to do step 2.

--
 
 

[Ben Goodrich]

On Tuesday, November 13, 2012 10:42:10 PM UTC-5, tds151 wrote:

thx, Ben. I would need to use int_step() to deal with multiple possible lp__ update nuggets based on whether 'y' == 1 or 'y==K' or a value in between.  I assume that's fine to do as a hack to implement a conditional statement.   i suppose it is, however, inefficient each of the 3 alternatives multiplied by an int_step() command are computed.

Yes, it would be ugly. But at least it would be ugly .stan code as opposed to ugly hacks to the C++ code along the lines of what Daniel mentioned to get past the parser error. Hacking the generated C++ code might be a bit better than hacking the Stan library if you only need to do a few models like this.

Ben

[Bob Carpenter]

I'm still not sure I understand the intention, so please
bear with my trying to rephrase it with more concrete types:

int<lower=0> N; // number of data items
int<lower=2> K; // max number of cats
int<lower=1,upper=K> y[N]; // outcome
int<lower=2,upper=K> ncats[N]; // number of possible cats

The constraint on y is tighter than can be declared in Stan,

y[n] in 1:ncats[n].

and here's the part we don't yet support:

ragged_matrix[K-1] c; // cutpoint vecs in rows

For instance, with K=5, we might have

c = [c1 ;
c2 c3 ;
c4 c5 c6 ;
c7 c8 c9 c10]

so that c[1] is a 1-row-vector [c1],
c[2] is the 2-row-vector [c2 c3], and so on.

And then what you want to do model-wise is something like this:

for (n in 1:N) {
eta[n] <- ...; // calculate predictor
y[n] ~ ordered_logistic(eta[n], c[ncuts[n]]');
}

Before getting into how Stan can deal with this, can you let
us know if this is your intent?

- Bob

> On Sunday, October 28, 2012 4:04:47 PM UTC-4, tds151 wrote:
>
> i'd really appreciate any advice on the forwarded topic, below, about dealing with a ragged array in
> stan. in my case, it is an array of cutpoints where the number of "rows" depends on the the column index.
>
> ---------- Forwarded message ----------
> From: *terrance savitsky* <tds...@gmail.com>
> Date: Sat, Oct 27, 2012 at 7:43 PM
> Subject: employing a ragged array of cutpoints
> To: stan-...@googlegroups.com
>
>
> I have an multvariate ordered outcome, y[N,D], where the number of categories depends on dimension d =

> 1, ..., D.I'd like to set an ordered_logistic likelihood for each y[i,d], where the length of each

[terrance savitsky]

yes, this outline expresses my situation.  you have expressed it in a more general and elegant manner to provide an example of the utility for a ragged array.   it would be wonderful to have the specification for a ragged_array in the fashion you outline.  this particular issue may also be resolved in a less elegant if partial indexing were available.  in this case, i might have a square matrix, but only use a portion of the entries.  

--

[Bob Carpenter]

It is too bad there's not an elegant way to do this
within Stan. We could handle this case with a triangular
matrix data type -- we wouldn't even need full ragged arrays.
We'll probably add one in the near future (months, not weeks
or days, that is).

As to some of the other proposed solutions, you don't want
unnecessary "junk" parameters as they'll distort the
posterior geometry and what the U-turn condition in
NUTS does (there are ways to mitigate their effect,
but it's still work by the algorithm).

Here's one solution that stays entirely within
Stan's modeling language and doesn't add any
unnecessary parameters. But it does make you recode things
like ordered vectors from scratch in the model.

Suppose you have K maximum possible outcome categories, which
means anywhere from 1 to (K-1) cut points per instance.
Suppose K=5, then I'm proposing a representation like this:

outcomes num_cuts lowest increments
2 1 first[1]
3 2 first[2] inc[1]
4 3 first[3] inc[2] inc[3]
5 4 first[4] inc[4] inc[5] inc[6]

Then the 3 cutpoints for the case of 4 outcomes
are (lowest[3], lowest[3]+inc[2], lowest[3]+inc[2]+inc[3]).
The rest is just a bunch of by-hand triangular array
indexing, which is majorly ugly.

There's some indexing arithmetic below to compute the
starting indexes 1, 2, and 4 for he increments.
Note that each is the sum of the number of entries in
the previous row, so takes the form of the sum of integers
from 1 to some number (plus 1 because we went with
R-style matrix indexing).

Anyway, here goes.


parameters {
real first_cuts[K-1];
real<lower=0> incr_cuts[((K-1)*(K-2))/ 2];
}
model {

// priors on cutpoints
first_cuts ~ ...;
incr_cuts ~ ...;

for (n in 1:N) {

// build up cutpoint vector
vector[num_cuts[n]] c;
c[1] <- first_cuts[num_cuts[n]];
for (k in 2:num_cuts[n])
c[k] <- c[k-1]
+ incr_cuts[((num_cuts[n]-1)*(num_cuts[n]-2))/2 + k - 1];

// linear predictor (local var or transformed eta[n] if saving)
eta <- ...;

// likelihood
y[n] ~ ordered_logistic(eta,c);

}

I haven't compiled this, so no guarantees on perfect
well-formedness.

The amount of work done to create new vectors (it requires
a malloc then free under the hood by our matrix lib) and all the
arithmetic and copying will probably be dwarfed by the
other operations, which include inverse logits (hence exponentials)
and their derivatives. (The allocations could be reduced
by sorting the data based on number of cut points and then
indexing into it.)

The next step would be to do what Ben suggests and just
rewrite the terms in the ordered logistic, which would require
a bit more fancy indexing footwork, but could be implemented
pretty directly.

An integer binomial coefficient function would make things
cleaner. All of this (n*(n+1))/2 stuff is just (n choose 2),
the sum of the numbers from 1 to n.

- Bob

> --
>
>