Re: [stan-users] Re: employing a ragged array of cutpoints

[Matt Hoffman]

(Moving to stan-dev.)

I think this also relates to the older question of being able to
inspect the length of an array and not having to specify it as a
separate variable. It'd be really nice to be able to write something
like
data {
int[int[]] x;
}
model {
for (i in 1:x.size)
for (j in 1:x[i].size)
x ~ normal(mu, 1);
}
to load an Rdump that looks like
x <-
list(c(1, 2, 3), 1, c(2, 8))

Not sure how high this goes on the to-do list, though.

Matt

On Mon, Oct 29, 2012 at 12:55 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> 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
>
>
> On 10/29/12 2:56 PM, Matt Hoffman wrote:
>>
>> 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
>>
>> On Sun, Oct 28, 2012 at 6:39 PM, Ben Goodrich <goodri...@gmail.com>
>> wrote:
>>>
>>> Actually, you would need to use max rather than fmax. Correct patch
>>> attached.
>>>
>>> Ben
>>>
>>>
>>> On Sunday, October 28, 2012 9:36:47 PM UTC-4, Ben Goodrich wrote:
>>>>
>>>>
>>>> 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
>>>>
>>>> On Sunday, October 28, 2012 9:12:33 PM UTC-4, tds151 wrote:
>>>>>
>>>>>
>>>>> 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.
>>>>>
>>>>> On Sun, Oct 28, 2012 at 1:10 PM, Ben Goodrich <goodri...@gmail.com>
>>>>> wrote:
>>>>>>
>>>>>>
>>>>>> 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
>>>>>>
>>>>>>
>>>>>> 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
>>>>>>> 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
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> Thank you, Terrance Savitsky
>>>>>>
>>>>>>
>>>>>> --
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Thank you, Terrance Savitsky
>>>
>>>
>>> --
>>>
>>>
>>
>
> --
>
>

[Bob Carpenter]

On 10/29/12 4:09 PM, Matt Hoffman wrote:
> (Moving to stan-dev.)
>
> I think this also relates to the older question of being able to
> inspect the length of an array and not having to specify it as a
> separate variable. It'd be really nice to be able to write something
> like
> data {
> int[int[]] x;
> }
> model {
> for (i in 1:x.size)
> for (j in 1:x[i].size)
> x ~ normal(mu, 1);
> }
> to load an Rdump that looks like
> x <-
> list(c(1, 2, 3), 1, c(2, 8))

It's much easier to see how to do this in the data { ... } block
than elsewhere. The data gets read out of a structured file
that declares dimensions (albeit inconveniently postifx).

The model transforms a vector of parameters into structured
parameters. I don't know how to do that without declared
sizes for structures.

For data, we'd need to rewrite the io/reader.hpp code to
just buffer all the input until it got to the dimension
declarations, then construct an object based on dims
or total number of objects read into which to write.

You can already do the following in Stan:

for (i in 1:size(x))
for (j in 1:size(x[i]))
...


I'd rather see a declaration that matched our
regular syntax like:

int x[,];

Now we COULD start generalizing all of our syntax so
it doesn't mix the syntax into the variable names like C
does. So for Stan 2.0, we could support a new format of:

int[5,3] x;

instead of the current:

int x[5,3];

or we could support both formats and even a mix of both
with

int[3] x[5]

or even

vector[2][3,4] x[5,6,7];

- Bob

[Matt Hoffman]

I think just having it in the data block would be fine. I can imagine
a couple of approaches to dealing with ragged parameter arrays. One
looks something like this:

data {
int x[,];
}
parameters {
float z[size(x),];

for (i in 1:size(x))

z[i].resize(size(x[i]));
}

and another looks something like this:

data {
int x[,];
}
parameters {
float z[x];
}

The second is pithier, but maybe more confusing. The first puts the
responsibility of defining array shapes on the user. Because the code
is in the parameters block, we only have to run it once, just after
we've read in the data and know what the shapes should be.

Matt

> --
>
>

[Bob Carpenter]

On 11/2/12 11:32 AM, Matt Hoffman wrote:
> I think just having it in the data block would be fine. I can imagine
> a couple of approaches to dealing with ragged parameter arrays. One
> looks something like this:
>
> data {
> int x[,];
> }
> parameters {
> float z[size(x),];
> for (i in 1:size(x))
> z[i].resize(size(x[i]));
> }

I don't see how to make this work given that we need to
know how to convert the parameter vector from HMC/NUTS
to the structured parameters of the model.

It probably could be finagled somehow. But I think your
next suggestion's more promising.

> and another looks something like this:
>
> data {
> int x[,];
> }
> parameters {
> float z[x];
> }
>
> The second is pithier, but maybe more confusing. The first puts the
> responsibility of defining array shapes on the user. Because the code
> is in the parameters block, we only have to run it once, just after
> we've read in the data and know what the shapes should be.

It's related to my proposal to have declarations be based
on sizes. So I'd prefer something like this in general:

data {
int Ns[];
}
parameters {
real z[Ns];
...
}

So size(z) = size(N), an size(z[i]) = Ns[i];

Then if we had an operator like sizes(), you could
generate what you want in a more direct way:

data {
int x[,];
}
parameters {

real z[sizes(x)];
}

where sizes(x) returns an array of sizes of x. It generalizes,
too, so you could have

int z[,,];

and sizes(z) would return an int[,].

Then we'd just think of

int w[I,J];

as being shorthand for the fully specified thing that
took an int[] indicating all the ragged sizes.

- Bob