matrix vector multiplication with diag_matrix
Marcus Brubaker
It's a simple AR model and the specification has the line:
lambda[t] ~ multi_normal( diag_matrix(1-phi)*mu + diag_matrix(phi)*lambda[t-1], U );
where lambda, phi and mu are vectors and U is a cov_matrix.
The parser chokes with:
lambda[t] ~ multi_normal( diag_matrix(1-phi)*mu + diag_matrix(phi)*lambda[t-1], U );
^-- here
DIAGNOSTIC(S) FROM PARSER:
expression is ill formed
Parser expecting: <eps>
Any suggestions? I've tried a number of different variations on the expression, including defining transformed variables to represent the mean to no avail.
Also, on a related note, has anyone created a time-series demo for stan? If not I'd be more than happy to contribute this one once I get it running.
Cheers,
Marcus
Bob Carpenter
stan-users, not that it's had any traffic yet.
I really need to write some doc! And work on some
better parser error messages.
The basic problem is that Stan's not very sophisticated
about the syntax for matrices. All the operations need
to be spelled out Lisp-like (I plan to add better parsing
for them at some point).
The expression for matrix product is "multiply(x,y)".
You can see the available ops in src/stan/math/matrix.hpp and
src/stan/agrad/matrix.hpp.
> lambda[t] ~ multi_normal( diag_matrix(1-phi)*mu +
> diag_matrix(phi)*lambda[t-1], U );
>
> where lambda, phi and mu are vectors and U is a cov_matrix.
To get (phi - 1), you need to either create a unit
vector one and use "minus(phi,one)", or write it
out longhand (which is just as fast as everything's
automatically "vectorized" when the generated code gets
compiled). For instance, you can do phi-1 this way if
phi is a parameter:
transformed parameters {
vector(K) phi_m_1; // K should be sample size
for (k in 1:K) phi_m_1[k] <- phi[k] - 1;
}
The same goes for the other matrix ops.
We've been working on some time-series data on
the side. But nothing's been contributed. We'll
set up some way to let people contribute models and
data so other people can run them. I think it'd
be a great resource.
- Bob
Marcus Brubaker
At any rate, thanks for the tips!
Bob Carpenter
I'll go check the settings for stan-dev and
make sure we require approval there, too.
We don't want to discourage anyone from contributing,
it's just a matter of where the discussions
should go.
I'll work on getting some more natural notation
into matrix arithmetic. op_name(x,y) looks very awkward
compared to (x op y).
- Bob
Jack Tanner
>
> I'll work on getting some more natural notation
> into matrix arithmetic. op_name(x,y) looks very awkward
> compared to (x op y).
Indeed. And x op y would be consistent with JAGS. See the manual
http://sourceforge.net/projects/mcmc-jags/files/Manuals/3.x/jags_user_manual.pdf/download
In general, is there an aim to keep Stan consistent with JAGS?
Bob Carpenter
or with JAGS (which aren't themselves consistent
with each other).
Our variables are typed. We have many more special
fucntions. We don't have the same parameterizations
of densities. We allow improper priors. We allow
direct manipulation of the total log probability
function so you can define your own densities or
transformations of existing densities (with Jacobian
determinants).
To be consistent with JAGS, we'd also need to use the funky
"%*%" notation for matrix operators, whereas I'm just
going to overload "*", etc.
We're behind in not having all the cumulative
densities (yet) and not having all the random generation
functions (again, as of yet). Those, we plan to add
as we go along.
- Bob
Jack Tanner
>
> To be consistent with JAGS, we'd also need to use the funky
> "%*%" notation for matrix operators, whereas I'm just
> going to overload "*", etc.
A plea: don't reinvent the wheel for this. Pick some existing language
convention (Matlab, R, JAGS, whatever) and use that for all operators. I
don't care that %*% is funky, because I only have to learn it once, and
I need to distinguish between element-wise ops and matrix-wise ops
somehow in any case.
Matt Hoffman
distinguish between elementwise multiplication and inner products, and
also that we should adopt some existing convention for that. As a
Matlab user, I find R's %*% notation a little obnoxious (it smells
like something that someone who doesn't do much matrix multiplication
might invent), but if we anticipate having more R users than Matlab
users then it might make sense to use that. Then there's Python's
approach, which is to do matrix multiplication using the dot()
function. Which is also a little ridiculous, but only because "dot()"
isn't a very accurate name.
Matt
Marcus Brubaker
Bob Carpenter
Python-like (or is that really numpy-like?)
and be done with it.
Seriously, though, I haven't at all been concerned
with compatibility with anything existing. I've tried to
be internally consistent, but we could use some review.
I went with full versions of distribution names,
like normal(0,1), instead of norm(0,1) or R/BUGS/JAGS
dnorm(0,1). And our parameterization uses scale, not
precision or squared scale.
I suppose we could add a precision-parameterized
normal(mu,tau)
and then add a scale-parameterized version,
normal_scale(mu,sigma);
and one with squared scale (aka variance), as in:
normal_squared_scale(mu,sigma_sq);
We also have the Cholesky-factor parameterization in the
multivariate case.
There are no element-wise ops, at least not yet. I've
reduced results where possible. So there are ops like:
double multiply(row_vector,vector);
matrix multiply(vector,row_vector);
But there are NO element-wise ops like:
vector multiply(vector,vector);
Should we add elementwise ops? And what should they be
called?
Does JAGS have element-wise ops? If I do (v1 * v2),
is that interpreted element-wise? Is that why they use the
funny %*% notation for matrix ops? I didn't see it in their
doc (I could've missed it or it could just not be in the doc).
I haven't implemented dim(matrix), so no:
vector<int> dim(matrix);
On the other hand, there are generic dot products,
double dot_product(vector,vector);
double dot_product(row_vector,vector);
but I spelled the name out in full.
We're also not allowing the left-hand-side assignment
algebra, as in:
logit(theta) <- x
and instead require:
theta <- inv_logit(x)
We just have to make up our minds before I start the doc
process.
- Bob
Andrew Gelman
A
On Mar 3, 2012, at 7:22 PM, Bob Carpenter wrote:
> Maybe we can call our currently function syntax
> Python-like (or is that really numpy-like?)
> and be done with it.
Saying "python"-like or matlab-like is fine with me. I have no desire to be bugs or jags compliant. We will have 50+ examples of Bugs/Jags code converted to Stan code so that shouldn't be an issue.
My current inclination is to use R-like syntax, but if everyone thinks that's too ugly we don't have to.
>
> Seriously, though, I haven't at all been concerned
> with compatibility with anything existing. I've tried to
> be internally consistent, but we could use some review.
>
> I went with full versions of distribution names,
> like normal(0,1), instead of norm(0,1) or R/BUGS/JAGS
> dnorm(0,1). And our parameterization uses scale, not
> precision or squared scale.
>
I like normal(), but I would not object to having norm() as an alias. As we discussed before, I _don't_ want to use dnorm(). The "d" there adds nothing but confusion.
> I suppose we could add a precision-parameterized
>
> normal(mu,tau)
>
I think the default should be the scale-parameterized version (as in R). I say that even though I use the variance-parameterized version in BDA. I do not not not not not want the so-called precision-parameterized version to be the default. But I agree that it would be ok to have alternative forms.
The alternative forms will be particularly important for inverse-chi^s dist and negative-binomial ("overdispersed Poisson").
>
> We also have the Cholesky-factor parameterization in the
> multivariate case.
For multivariate, I think the default parameterization should use the covariance matrix.
>
> There are no element-wise ops, at least not yet. I've
> reduced results where possible. So there are ops like:
>
> double multiply(row_vector,vector);
>
> matrix multiply(vector,row_vector);
>
> But there are NO element-wise ops like:
>
> vector multiply(vector,vector);
>
> Should we add elementwise ops? And what should they be
> called?
Personally, I'm happy with the R convention that +, -, *, / operate component-wise. I find that to be very convenient. But if it's confusing, we don't need it.
>
> Does JAGS have element-wise ops? If I do (v1 * v2),
> is that interpreted element-wise? Is that why they use the
> funny %*% notation for matrix ops? I didn't see it in their
> doc (I could've missed it or it could just not be in the doc).
>
Again, I don't see the point of being Jags-compliant, I'd prefer to be consistent with R or Matlab or Python or whatever.
> I haven't implemented dim(matrix), so no:
>
> vector<int> dim(matrix);
>
> On the other hand, there are generic dot products,
>
> double dot_product(vector,vector);
> double dot_product(row_vector,vector);
>
> but I spelled the name out in full.
>
> We're also not allowing the left-hand-side assignment
> algebra, as in:
>
> logit(theta) <- x
>
> and instead require:
>
> theta <- inv_logit(x)
That's fine with me. I've never liked the "logit(theta) <- x" notation. This relates to my desire for our code to be executable code rather than just an interpreted model. One thing that I found difficult about Bugs is that its model is _not_ actually code that runs in order. Also, in my books I used the invlogit notation. (I actually call it invlogit() in R, but if you prefer inv_logit, that's fine. Some people write expit, which makes logical sense but I've never used it, as it seems like just one more thing to learn.
Bob Carpenter
> Hi--see below for my comments.
>
>
>> Maybe we can call our currently function syntax
>> Python-like (or is that really numpy-like?)
>> and be done with it.
>
> Saying "python"-like or matlab-like is fine with me. I have no desire to be bugs or jags compliant. We will have 50+ examples of Bugs/Jags code converted to Stan code so that shouldn't be an issue.
I agree, but all else being equal, doing what R does
seems like a not unreasonable thing to do.
Does R support matrix multiplication in some way with
syntax like (m1 %*% m2) a la JAGS?
> My current inclination is to use R-like syntax, but if everyone thinks that's too ugly we don't have to.
We're basically using R-like syntax with "<-" for assignments
and simple loops being "(for n in 1:N)". But you forgot the
semicolons!
I do not want to give up the semicolons. I find R, BUGS, and JAGS's
context-dependent use of whitespace really really annoying.
By that I mean that it's OK to write:
a <- b +
c
but not OK to write:
a <- b
+ c
The latter is the standard convention in mathematical typesetting
and the convention adopted by almost all coding conventions I've ever
seen. Here I want to follow C, not R/BUGS/JAGS.
>> Seriously, though, I haven't at all been concerned
>> with compatibility with anything existing. I've tried to
>> be internally consistent, but we could use some review.
>>
>> I went with full versions of distribution names,
>> like normal(0,1), instead of norm(0,1) or R/BUGS/JAGS
>> dnorm(0,1). And our parameterization uses scale, not
>> precision or squared scale.
>>
>
> I like normal(), but I would not object to having norm() as an alias. As we discussed before, I _don't_ want to use dnorm(). The "d" there adds nothing but confusion.
I'd rather not have things aliased. I think that leads
to confusion when you look at someone else's code and
wonder if it means the same as you.
>> I suppose we could add a precision-parameterized
>>
>> normal(mu,tau)
>>
>
> I think the default should be the scale-parameterized version (as in R). I say that even though I use the variance-parameterized version in BDA. I do not not not not not want the so-called precision-parameterized version to be the default. But I agree that it would be ok to have alternative forms.
Right -- you've said that before, which is why we went with it.
> The alternative forms will be particularly important for inverse-chi^s dist and negative-binomial ("overdispersed Poisson").
We have a couple alternative parameterizations now, but will
add more as we go. Ben wants to add a bunch of useful multivariate
parameterizations.
And I think it'd be really nice from a speed perspective to
use some unconstrained outcome parameterizations. That would
avoid having to transform all the variables and apply the
absolute derivative of the transform independently -- we can
roll it into the density and reduce it.
We'd also like to add more compound operations like
log_sum_exp() and log1p() and multiply_log() [the latter
of which we needed for a * log(b) when both a and b were 0].
>> We also have the Cholesky-factor parameterization in the
>> multivariate case.
>
> For multivariate, I think the default parameterization should use the covariance matrix.
That's the default for now. Other than for the normal,
I followed BDA.
>> There are no element-wise ops, at least not yet. I've
>> reduced results where possible. So there are ops like:
>>
>> double multiply(row_vector,vector);
>>
>> matrix multiply(vector,row_vector);
>>
>> But there are NO element-wise ops like:
>>
>> vector multiply(vector,vector);
>>
>> Should we add elementwise ops? And what should they be
>> called?
>
> Personally, I'm happy with the R convention that +, -, *, / operate component-wise. I find that to be very convenient. But if it's confusing, we don't need it.
So that's a yes, we should add componentwise ops? I'll put
it on the to-do list.
I can overload the basic arithmetic to do component-wise ops
pretty easily.
But what about operator notation for matrices? %*%, etc.?
I find those incredibly ugly. Not considering precedent,
I'd prefer to make the matrix ops the defaults and write something
else for component-wise ops. Sounds like that'll confuse
everyone, though.
>> Does JAGS have element-wise ops? If I do (v1 * v2),
>> is that interpreted element-wise? Is that why they use the
>> funny %*% notation for matrix ops? I didn't see it in their
>> doc (I could've missed it or it could just not be in the doc).
>>
>
> Again, I don't see the point of being Jags-compliant, I'd prefer to be consistent with R or Matlab or Python or whatever.
The only point is to make it easy for people to move.
It's a double-edged sword, though. If we make it too
much like, but not exactly like something else, it
could be even more confusing.
>> I haven't implemented dim(matrix), so no:
>>
>> vector<int> dim(matrix);
>>
>> On the other hand, there are generic dot products,
>>
>> double dot_product(vector,vector);
>> double dot_product(row_vector,vector);
>>
>> but I spelled the name out in full.
>>
>> We're also not allowing the left-hand-side assignment
>> algebra, as in:
>>
>> logit(theta)<- x
>>
>> and instead require:
>>
>> theta<- inv_logit(x)
>
> That's fine with me. I've never liked the "logit(theta)<- x" notation. This relates to my desire for our code to be executable code rather than just an interpreted model. One thing that I found difficult about Bugs is that its model is _not_ actually code that runs in order. Also, in my books I used the invlogit notation. (I actually call it invlogit() in R, but if you prefer inv_logit, that's fine. Some people write expit, which makes logical sense but I've never used it, as it seems like just one more thing to learn.
Again, it'd be easy enough to add to the parser and it'd
just be short hand for applying the inverse op on the RHS.
But I always found it confusing in use.
- Bob
Andrew Gelman
>
> I agree, but all else being equal, doing what R does
> seems like a not unreasonable thing to do.
Yes, I agree. We could take R or Matlab or Python as our model. I don't think it makes sense to use Bugs and Jags. Once Stan is working, people won't need Bugs or Jags, they'll be forgotten. But R/Matlab/Python will still be around.
>
> Does R support matrix multiplication in some way with
> syntax like (m1 %*% m2) a la JAGS?
R does matrix mulitplication using %*%. I assume that's where Jags got it from.
>
> We're basically using R-like syntax with "<-" for assignments
> and simple loops being "(for n in 1:N)". But you forgot the
> semicolons!
I can live with the semicolons!
>
> I'd rather not have things aliased. I think that leads
> to confusion when you look at someone else's code and
> wonder if it means the same as you.
>
No prob. The dist names should be clear enough from the examples, anyway.
>>>
>>> Should we add elementwise ops? And what should they be
>>> called?
>>
>> Personally, I'm happy with the R convention that +, -, *, / operate component-wise. I find that to be very convenient. But if it's confusing, we don't need it.
>
> So that's a yes, we should add componentwise ops? I'll put
> it on the to-do list.
>
> I can overload the basic arithmetic to do component-wise ops
> pretty easily.
>
> But what about operator notation for matrices? %*%, etc.?
> I find those incredibly ugly. Not considering precedent,
> I'd prefer to make the matrix ops the defaults and write something
> else for component-wise ops. Sounds like that'll confuse
> everyone, though.
For addition and subtraction, there's no ambiguity. So the only real issue is multiplication.
>>>
>>> We're also not allowing the left-hand-side assignment
>>> algebra, as in:
>>>
>>> logit(theta)<- x
>>>
>>> and instead require:
>>>
>>> theta<- inv_logit(x)
>>
>> That's fine with me. I've never liked the "logit(theta)<- x" notation. This relates to my desire for our code to be executable code rather than just an interpreted model. One thing that I found difficult about Bugs is that its model is _not_ actually code that runs in order. Also, in my books I used the invlogit notation. (I actually call it invlogit() in R, but if you prefer inv_logit, that's fine. Some people write expit, which makes logical sense but I've never used it, as it seems like just one more thing to learn.
>
> Again, it'd be easy enough to add to the parser and it'd
> just be short hand for applying the inverse op on the RHS.
> But I always found it confusing in use.
I don't like it because it leads people to read the lines of the model as mathematical expressions rather than computable expressions. So if you're in agreement on this, let's not do it...