Vectorized distributions

[Daniel Lee]

Just so we have it on hand, the list of vectorized and tested distributions are:
- beta
- cauchy

[Daniel Lee]

- double_exponential
- normal

[Daniel Lee]

Just so we have it on hand, the list of vectorized and tested distributions are:
- beta
- cauchy
- double_exponential
- normal

- bernoulli
- bernoulli_logit


I'm going to have to extend the parser to allow for int vector types,
but that shouldn't be a problem. At least I don't think so.

[Bob Carpenter]

On 9/22/12 3:34 AM, Daniel Lee wrote:

> I'm going to have to extend the parser to allow for int vector types,
> but that shouldn't be a problem.

You mean a type that would translate to
Eigen's Matrix<int,Dynamic,1>?

What would it be used for?

Totally doable, of course. It would require
adding the new base type and type inference
rules to ast.hpp. And new I/O routines everywhere.

Int types are a bit simpler in that there's
no transform needed for inits or reverse transform
needed for parameters.

We'd talked about having something like

vector<basetype=int>

vector<basetype=real> == vector<> == vector

The other alternative would be to introduce
types int_vector, int_row_vector, int_matrix
and leave vector, row_vector and matrix as is.

- Bob

[Bob Carpenter]

On 9/22/12 3:34 AM, Daniel Lee wrote:

> Just so we have it on hand, the list of vectorized and tested distributions are:

...
> - bernoulli
> - bernoulli_logit

Those last two are new and just what Andrew
wanted. They should speed up logistic
regression the same way they sped up linear
regression.

That must've taken a heroic amount of
derivative calculation and coding to do this
since yesterday afternoon!

I'll change the doc to match ASAP. I introduced
a convention for vectorized distributions to use
"reals" for a vectorized real type (I stole
the natural language plural convention from Prolog,
who I believe stole it from Lisp).

Also, everyone else -- I managed to generate
a Stan function index. The only problem is that
the underscore isn't searchable in Preview in our
docs. Who knows what pdflatex is generating.
PDFs are suboptimal for search, as I've
noted before for Acrobat (Matt -- please go kick some
corporate butt on this matter):

http://lingpipe-blog.com/2007/03/18/a-pdfn-mess/

- Bob

[Ben Goodrich]

On Saturday, September 22, 2012 4:10:49 PM UTC-4, Bob Carpenter wrote:

Also, everyone else -- I managed to generate
a Stan function index.  The only problem is that
the underscore isn't searchable in Preview in our
docs.  Who knows what pdflatex is generating.
PDFs are suboptimal for search, as I've
noted before for Acrobat (Matt -- please go kick some
corporate butt on this matter):

http://lingpipe-blog.com/2007/03/18/a-pdfn-mess/

 
FWIW, I can find "fulfills" in that PDF using okular, which apparently can also be installed on inferior operating systems.

Ben

[Bob Carpenter]

On 9/22/12 4:39 PM, Ben Goodrich wrote:
> On Saturday, September 22, 2012 4:10:49 PM UTC-4, Bob Carpenter wrote:
>
> Also, everyone else -- I managed to generate
> a Stan function index. The only problem is that
> the underscore isn't searchable in Preview in our
> docs. Who knows what pdflatex is generating.
> PDFs are suboptimal for search, as I've
> noted before for Acrobat (Matt -- please go kick some
> corporate butt on this matter):
>

> http://lingpipe-blog.com/2007/03/18/a-pdfn-mess/ <http://lingpipe-blog.com/2007/03/18/a-pdfn-mess/>

>
>
> FWIW, I can find "fulfills" in that PDF using okular, which apparently can also be installed on inferior operating systems.

Works in Preview, too. I guess they've
improved these things in the past five years.

- Bob

[Andrew Gelman]

Just tell me when they're in the current downloadable stan and I'll give them a try.
A

> --
>

[Ben Goodrich]

On Saturday, September 22, 2012 4:10:49 PM UTC-4, Bob Carpenter wrote:

Also, everyone else -- I managed to generate
a Stan function index.  The only problem is that
the underscore isn't searchable in Preview in our
docs.  Who knows what pdflatex is generating.
PDFs are suboptimal for search, as I've
noted before for Acrobat (Matt -- please go kick some
corporate butt on this matter):

http://lingpipe-blog.com/2007/03/18/a-pdfn-mess/

 
Does adding \usepackage[T1]{fontenc} to the preamble help with Preview?

http://www.latex-community.org/forum/viewtopic.php?f=48&t=8891

Ben

[Daniel Lee]

On Sep 22, 2012 5:22 PM, "Andrew Gelman" <gel...@stat.columbia.edu> wrote:
>
> Just tell me when they're in the current downloadable stan and I'll give them a try.
> A

They will be available in 1.0.2. ETA is next Stan meeting +/- a day.

[Daniel Lee]

I might have been thinking about the problem wrong.

What I am trying to do is add the function signatures for
bernoulli_log in stan/gm/function_signatures.hpp.

The function signature of bernoulli_log looks like:

typename return_type<T_prob>::type
bernoulli_log(const T_n& n, const T_prob& theta)

where T_n can be an int, std::vector<int>, Matrix<int,Dynamic,1>, or
Matrix<int,1,Dynamic> and T_prob can be double, std::vector<double>,
Matrix<double,Dynamic,1>, Matrix<double,1,Dynamic> or with vars
instead of doubles. That's all been tested on the C++ side and it
works. I just need to add the vectorized function signatures for
stanc.

Does this cover all the function signatures I need to add?

for (size_t i = 0; i < vector_types.size(); ++i) {
add("bernoulli_log",DOUBLE_T,INT_T,vector_types[i]);
add("bernoulli_log",DOUBLE_T,expr_type(INT_T,1U),vector_types[i]);
}

Or should the 3rd argument for add() also include some representation
equivalent to VECTOR_T and ROW_T (like you suggested above with the
int_vector and int_row_vector)?

I don't think I was really clear how I explained what I'm trying to
do, but hopefully there's enough context to figure it out. If not, we
can talk Monday. I'm going to hammer out more of these vectorized
distributions this weekend.




Daniel

[Bob Carpenter]

On 9/22/12 5:45 PM, Daniel Lee wrote:

> What I am trying to do is add the function signatures for
> bernoulli_log in stan/gm/function_signatures.hpp.
>
> The function signature of bernoulli_log looks like:
>
> typename return_type<T_prob>::type
> bernoulli_log(const T_n& n, const T_prob& theta)
>
> where T_n can be an int, std::vector<int>, Matrix<int,Dynamic,1>, or
> Matrix<int,1,Dynamic> and T_prob can be double, std::vector<double>,
> Matrix<double,Dynamic,1>, Matrix<double,1,Dynamic> or with vars
> instead of doubles.

Nice -- that makes sense on the C++ side.

> That's all been tested on the C++ side and it
> works. I just need to add the vectorized function signatures for
> stanc.
>
> Does this cover all the function signatures I need to add?
>
> for (size_t i = 0; i < vector_types.size(); ++i) {
> add("bernoulli_log",DOUBLE_T,INT_T,vector_types[i]);
> add("bernoulli_log",DOUBLE_T,expr_type(INT_T,1U),vector_types[i]);
> }

Yes, that covers what we have now in the way of basic
int and vectors of ints.

> Or should the 3rd argument for add() also include some representation
> equivalent to VECTOR_T and ROW_T (like you suggested above with the
> int_vector and int_row_vector)?

We'd need to add those types, and unless there are
vector ops we need to do on them, I don't see the
point (other than that Matrix is a bit tighter
in storage than vectors of vectors, so should be
a bit faster in the 2-D case).

> ... I'm going to hammer out more of these vectorized
> distributions this weekend.

That's great. I'm working on math/matrix and agrad/matrix.
There were some holes in what we exposed and different
holes in what we tested, but I'm filling them all in.

I'm thinking it'll also be worth vectorizing the rest of
the functions over reals, which mostly means things that
would then return vectors. This should still be more efficient
if we get around to writing custom auto-dif for the ops.

- Bob

[Daniel Lee]

On Sat, Sep 22, 2012 at 6:32 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> On 9/22/12 5:45 PM, Daniel Lee wrote:
>
>> What I am trying to do is add the function signatures for
>> bernoulli_log in stan/gm/function_signatures.hpp.
>>
>> The function signature of bernoulli_log looks like:
>>
>> typename return_type<T_prob>::type
>> bernoulli_log(const T_n& n, const T_prob& theta)
>>
>> where T_n can be an int, std::vector<int>, Matrix<int,Dynamic,1>, or
>> Matrix<int,1,Dynamic> and T_prob can be double, std::vector<double>,
>> Matrix<double,Dynamic,1>, Matrix<double,1,Dynamic> or with vars
>> instead of doubles.
>
>
> Nice -- that makes sense on the C++ side.
>
>
>> That's all been tested on the C++ side and it
>> works. I just need to add the vectorized function signatures for
>> stanc.
>>
>> Does this cover all the function signatures I need to add?
>>
>> for (size_t i = 0; i < vector_types.size(); ++i) {
>> add("bernoulli_log",DOUBLE_T,INT_T,vector_types[i]);
>> add("bernoulli_log",DOUBLE_T,expr_type(INT_T,1U),vector_types[i]);
>> }
>
>
> Yes, that covers what we have now in the way of basic
> int and vectors of ints.

Awesome! That makes life easy.

>> Or should the 3rd argument for add() also include some representation
>> equivalent to VECTOR_T and ROW_T (like you suggested above with the
>> int_vector and int_row_vector)?
>
>
> We'd need to add those types, and unless there are
> vector ops we need to do on them, I don't see the
> point (other than that Matrix is a bit tighter
> in storage than vectors of vectors, so should be
> a bit faster in the 2-D case).

That makes sense. We'll hold off on that and implement sometime in the
future, if necessary.

>> ... I'm going to hammer out more of these vectorized
>> distributions this weekend.
>
>
> That's great. I'm working on math/matrix and agrad/matrix.
> There were some holes in what we exposed and different
> holes in what we tested, but I'm filling them all in.
>
> I'm thinking it'll also be worth vectorizing the rest of
> the functions over reals, which mostly means things that
> would then return vectors. This should still be more efficient
> if we get around to writing custom auto-dif for the ops.
>

Sounds like a good thing to do.

[Daniel Lee]

On Sat, Sep 22, 2012 at 3:34 AM, Daniel Lee <bea...@alum.mit.edu> wrote:
> Vectorized and tested distributions are:

> - beta
> - cauchy
> - double_exponential
> - normal
> - bernoulli
> - bernoulli_logit

- gamma
- inv_gamma

[Bob Carpenter]

I just installed Xcode 4.5 and haven't had
any problem running our code or tests.

- Bob

[Daniel Lee]

Vectorized and tested distributions are:
- beta
- cauchy
- double_exponential
- normal
- bernoulli
- bernoulli_logit
- gamma
- inv_gamma

new:
- logistic

[Daniel Lee]

Vectorized and tested distributions are:
- beta
- cauchy
- double_exponential
- normal
- bernoulli
- bernoulli_logit
- gamma
- inv_gamma

- logistic

new:
- lognormal

[Daniel Lee]

Vectorized and tested distributions are:
- beta
- cauchy
- double_exponential
- normal
- bernoulli
- bernoulli_logit
- gamma
- inv_gamma
- logistic

- lognormal
- pareto

I'm going to shift focus on vectorizing all the distributions without
doing anything clever, then walk through the distributions evaluating
the partials.


Daniel

[Ben Goodrich]

 
Daniel, how difficult would it be to make

y ~ normal(scalar, scalar);

work for a multidimensional y, like a matrix or array? It may already work in principle (but isn't exposed to the parser) if VectorView works for whatever y is.

Thanks,
Ben

[Daniel Lee]

That works and should be in v1.0.2. If you are having any difficulties, can you attach an example and I will figure out what wasn't implemented correctly.

We should have tests in test/gm/model_specs/function_signatures*.stan that flex every combination of arguments. 

Thanks,
Ben

--
 
 

[Ben Goodrich]

goodrich@CYBERPOWERPC:/tmp/2D$ cat twoD.stan
parameters {
  real y[3,2]; // or matrix[3,2] y;
}
model {
  y ~ normal(0.0,1.0);
}

goodrich@CYBERPOWERPC:/opt/stan$ make O=0 /tmp/2D/twoD

--- Translating Stan graphical model to C++ code ---
bin/stanc /tmp/2D/twoD.stan --o=/tmp/2D/twoD.cpp
Model name=twoD_model
Input file=/tmp/2D/twoD.stan
Output file=/tmp/2D/twoD.cpp

INVALID COMMAND-LINE ARGUMENT
LOCATION:  file=/tmp/2D/twoD.stan; line=5, column=3

  y ~ normal(0.0,1.0);
  ^-- here

DIAGNOSTIC(S) FROM PARSER:

no matches for function name="normal_log"
    arg 0 type=real[2]
    arg 1 type=real
    arg 2 type=real
unknown distribution=normal
Parser expecting: "}"

 
Very tangentially related. Is it possible in the makefiles to just ignore the optimization level for the bin/stanc{.exe} target only? Am I correct in thinking that using -O3 is essentially useless for the parser and just takes more time and RAM to compile? The only thing that makes me concerned is that make bin/stanc also compiles libstan where I guess that O3 is important.

Thanks,
Ben

[Daniel Lee]

I haven't thought about how to generalize to multidimensional matrices
and 2d arrays yet. Can we define a consistent behavior for the
distributions? I'm guessing that you want:

matrix ~ normal(scalar, scalar)

to mean that each element in the matrix should be treated as y in the
distribution. That makes sense, but then what should be the
appropriate behaviors for things like:
matrix ~ normal(matrix, vector)
or
matrix ~ normal(vector, row_vector)
or
matrix[3,2] ~ normal(matrix[3,1], matrix[1,2]) (where dim is in the
square brackets)

If we can define consistent behavior, we can try to beef up the
VectorView stuff to do what we want. I really haven't thought about
these things too much, so maybe there's a simple answer to all of
this.

[Ben Goodrich]

On Friday, October 12, 2012 12:16:36 PM UTC-4, Daniel Lee wrote:

I haven't thought about how to generalize to multidimensional matrices
and 2d arrays yet. Can we define a consistent behavior for the
distributions? I'm guessing that you want:

matrix ~ normal(scalar, scalar)

to mean that each element in the matrix should be treated as y in the
distribution.

Yes. Just flatten whatever is on the left-hand side.
 

That makes sense, but then what should be the
appropriate behaviors for things like:
matrix ~ normal(matrix, vector)
or
matrix ~ normal(vector, row_vector)
or

I would say throw exception unless everything is the same size or is a scalar.
 

matrix[3,2] ~ normal(matrix[3,1], matrix[1,2]) (where dim is in the
square brackets)

I would have guessed that that was already legal, but if not it should be.
 

If we can define consistent behavior, we can try to beef up the
VectorView stuff to do what we want. I really haven't thought about
these things too much, so maybe there's a simple answer to all of
this.

I don't think there is a simple answer for the

multidimensional ~ normal(multidimensional, multidimensional)

case unless all the multidimensional objects are the same dimension. R recycles and throws a warning if recycling doesn't come out right, but I think Stan should only recycle scalars.

Ben

[Daniel Lee]

On Oct 12, 2012, at 9:38 AM, Ben Goodrich <goodri...@gmail.com> wrote:

On Friday, October 12, 2012 12:16:36 PM UTC-4, Daniel Lee wrote:

I haven't thought about how to generalize to multidimensional matrices
and 2d arrays yet. Can we define a consistent behavior for the
distributions? I'm guessing that you want:

matrix ~ normal(scalar, scalar)

to mean that each element in the matrix should be treated as y in the
distribution.

Yes. Just flatten whatever is on the left-hand side.

Would we also flatten the others parameters as well? We currently make no distinction between the arguments in the c++ distribution function. I just want to be clear as to what happens when the lhs is a 3x2 matrix and a parameter is a 2x3 matrix. Or perhaps a 3 vector or a 6 vector. 

That makes sense, but then what should be the
appropriate behaviors for things like:
matrix ~ normal(matrix, vector)
or
matrix ~ normal(vector, row_vector)
or

I would say throw exception unless everything is the same size or is a scalar.
 

matrix[3,2] ~ normal(matrix[3,1], matrix[1,2]) (where dim is in the
square brackets)

I would have guessed that that was already legal, but if not it should be.
 

We don't allow the lhs, but this is currently illegal. The lengths of the parameters do not match. One is length 3 and the other is length 2. 

If we can define consistent behavior, we can try to beef up the
VectorView stuff to do what we want. I really haven't thought about
these things too much, so maybe there's a simple answer to all of
this.

I don't think there is a simple answer for the

multidimensional ~ normal(multidimensional, multidimensional)

case unless all the multidimensional objects are the same dimension. R recycles and throws a warning if recycling doesn't come out right, but I think Stan should only recycle scalars.

 

What does "recycling" mean?

[Ben Goodrich]

On Saturday, October 13, 2012 7:19:11 PM UTC-4, Daniel Lee wrote:

On Oct 12, 2012, at 9:38 AM, Ben Goodrich <goodri...@gmail.com> wrote:

On Friday, October 12, 2012 12:16:36 PM UTC-4, Daniel Lee wrote:

I haven't thought about how to generalize to multidimensional matrices
and 2d arrays yet. Can we define a consistent behavior for the
distributions? I'm guessing that you want:

matrix ~ normal(scalar, scalar)

to mean that each element in the matrix should be treated as y in the
distribution.

Yes. Just flatten whatever is on the left-hand side.

Would we also flatten the others parameters as well? We currently make no distinction between the arguments in the c++ distribution function. I just want to be clear as to what happens when the lhs is a 3x2 matrix and a parameter is a 2x3 matrix.

I would vote for throwing an exception.
 

Or perhaps a 3 vector or a 6 vector.

I could vote for throwing an exception here too, but that seems like arguably valid input if it is not too hard to implement. The 3 vector could be interpreted as applying to both columns and the 6 vector as applying sequentially to the six matrix entries. But if that is hard, then don't worry about it. It is not too much to ask users to make the input match or be scalars.
 

 

That makes sense, but then what should be the
appropriate behaviors for things like:
matrix ~ normal(matrix, vector)
or
matrix ~ normal(vector, row_vector)
or

I would say throw exception unless everything is the same size or is a scalar.
 

matrix[3,2] ~ normal(matrix[3,1], matrix[1,2]) (where dim is in the
square brackets)

I would have guessed that that was already legal, but if not it should be.
 

We don't allow the lhs, but this is currently illegal. The lengths of the parameters do not match. One is length 3 and the other is length 2. 

Incompatible lengths should definitely result in a fatal exception. I mistakenly thought the above was a weird way to say scalar ~ normal(scalar, scalar).
 

If we can define consistent behavior, we can try to beef up the
VectorView stuff to do what we want. I really haven't thought about
these things too much, so maybe there's a simple answer to all of
this.

I don't think there is a simple answer for the

multidimensional ~ normal(multidimensional, multidimensional)

case unless all the multidimensional objects are the same dimension. R recycles and throws a warning if recycling doesn't come out right, but I think Stan should only recycle scalars.

 

What does "recycling" mean?

Recycling means repeating the elements of a vector until they match the length of a longer vector. In other words, this is technically legal in R

> set.seed(12345)
> rnorm(6, mean = 0:2, sd = 4:5) # implicit recycling
[1]  2.342115  4.547330  1.562787 -2.267486  3.423550 -7.089780
> set.seed(12345)
> rnorm(6, mean = c(0:2, 0:2), sd = c(4:5, 4:5, 4:5)) # explicit recycling
[1]  2.342115  4.547330  1.562787 -2.267486  3.423550 -7.089780

where the colon means a sequence of integers from the left to the right. But I think implicit recycling is too dangerous for Stan, except recycling scalars to match the length of a vector / matrix / array.

Ben

 

[Bob Carpenter]

Option 0. The way I would think about this is in terms
of coercing the various inputs to sequences,
throwing an exception if they don't match in
sequence size, and repeating the function as many times
as they do.

scalar -- recycle to max length
std::vector -- natural sequence
Eigen column vec -- natural sequence
Eigen row vec -- natural sequence

This is what we have now

Option 1. Extend this in the obvious way to matrices:

Eigen matrix -- column major order

This would match what we do now and also match how
Eigen defines operator()(int) for matrices.

If we do add them, I think it'll be simplest to just
define the natural sequentialization and allow
anything with matching sequence size. This matches
Eigen's unary operator()(int) behavior.

I agree that recycling is a bad idea for anything
other than a scalar. The only case I see it working
is recycling a row or column vector up to a matrix.
Or subarrays up to full arrays, like allowing
a CxD array to recycle up to an AxBxCxD array.
But that's starting to get awfully subtle.

Option 2. The only other thing that seems simple is to require
all the types to match other than scalar. But I
think everyone wants to allow 1D arrays, vectors
and row vectors to mix.

The major downside of this proposal is that you
could put in a 3 x 2 and a 2 x 3 matrix, and they'd
get converted to sequences as follows:

a[1,1], a[2,1], a[3,1], a[2,1], a[2,2], a[3,2]

b[1,1], b[2,1], b[1,2], b[2,2], b[1,3], b[2,3]

I can't imagine anywhere a user would do this and
it wouldn't be an error.

Option 3. The in between spec would be to require the sizes
to all match, and if there are matrices and/or arrays,
their dimensions all have to match.

I'm OK with any of these options, including not extending
to matrices.

Also, what about 2D, 3D or greater arrays?

- Bob

[Ben Goodrich]

On Saturday, October 13, 2012 8:11:52 PM UTC-4, Bob Carpenter wrote:

Option 3.  The in between spec would be to require the sizes
to all match, and if there are matrices and/or arrays,
their dimensions all have to match. 

I think I would be in favor of option 3. Almost anything should be allowed on the LHS, scalars should always be allowed on the RHS, but if a parameter on the RHS is not a scalar, it should match the dimensionality of the LHS and not just the size.

Ben

[Bob Carpenter]

I think what's allowed on the left and right
should be symmetric (not yet considering when
there are multivariate quantities on the left
or right).

So this would mean:

OK: normal(matrix[m,n], real, real)
normal(real, matrix[m,n], real) (*)
normal(real, real, matrix[m,n])

BAD: normal(matrix[m,n], vector[m*n], real)
normal(matrix[m,n], matrix[n,m], real)

OK: normal(matrix[m,n], matrix[m,n], real);
normal(matrix[m,n], matrix[m,n], matrix[m,n]);

etc.

I understand that it's not going to be likely to see
the case (*) above, but I think it's easier to
keep all the implementations allowing it.

What about arrays?

??: normal(real[m,n], matrix[m,n], real)

- Bob

[Daniel Lee]

Vectorized and tested distributions are:
- beta
- cauchy
- double_exponential
- normal
- bernoulli
- bernoulli_logit
- gamma
- inv_gamma
- logistic
- lognormal
- pareto

slow implementations:
- chi_square
- exponential
- inv_chi_square
- scaled_inv_chi
- student_t
- uniform
- weibull

CDFs (slow):
- normal_cdf