new density notation

[Bob Carpenter]

I want to run this by everyone on stan-dev before
creating the issue and starting to work on it.

1. New density expression with normalization and vertical bar:

cauchy<norm=true>(y | mu, sigma)

defaults to norm=false:

cauchy(y | mu, sigma)


2. New increment log density statement:

ld << cauchy(tau | 0, 2.5),
normal(beta | 0, 2),
bernoulli_logit(y | x * (tau * beta));


comma-separated values all added to log density accumulator.


3.a. Deprecate the use of any function ending in "_log" and
go with the above form.

b. Maintain support for user-defined functions ending in "_log"
to be used this way and keep our density library in C++ ending in "_log".

c. Make a local variable available to functions such as norm__
with the value of the <norm=...> specification (0 or 1 as int).


4. Deprecate use of the sampling statement:

y ~ foo(theta);


5. Rewrite all the examples in the manual!


NOTE ON BACKWARD COMPATIBILITY: The proposal is to maintain backward
compatibility but raise deprecation warnings. Alternatively, we could
just capture old syntax and print specialized update warnings.

- Bob

[Aki Vehtari]

I like | and <<, but have a question:
Would
a <- cauchy(y|mu,sigma);
return density or log-density? And what form is used if you want the other one?

Aki

[Krzysztof Sakrejda]

I like it, especially getting rid of sampling statements. I would like it even more if the notation was just = for assignment and += for increment and "ld" was just special because it was reserved. Seems overspecified to have a special increment syntax for a special variable... Or is the proposal to allow << to increment anything?

[Sebastian Weber]

Hi Bob!

This sounds like a huge thing - is it worth investing this effort? Also, the deprecation period is hopefully very long, i.e. Stan programs will look very different wrt to BUGS programs with this.

My comments:

1. Why introduce another logic with "<norm=true/false>"? Do you plan to introduce further options here, then it makes sense. If only for discriminating between normed/unnormed I would prefer a _norm naming convention thing.

2. What is the rationale to introduce two new operators "<<" and ","? Maybe I haven´t understood fully, but writing

lp << density1 << density2;

makes more sense to me.

In any case, I suggest to try out this new scheme on people who are more statisticians than programmers, i.e. I think someone having a programming background will absorb this easily, but I am not so sure about those people without a programming background.

Best,
Sebastian

[Bob Carpenter]

I was thinking about just that.

* introduce "=" as a synonym for "<-" in assignment statements
and deprecate "<-"

* use "+=" for increment statements

* introduce operators "++" and "--" for expressions
- definitely prefix
- optionally postfix

Stan already includes a deprecated variable "lp__" which holds
the log density accumulator. We don't want people to assign
to that variable --- too easy to lose terms, including the Jacobian.
So we just want an increment. And we want to access its value,
and I thought we'd added a function, but I don't see it; so we need

* function log_density() that returns current value of accumulator

And then we can get rid of lp__ altogether in Stan 3.

There's the remaining question of whether we

* reserve ld as a variable or function name

It won't conflict in the streaming increment notation.
The "<<" is C++ streaming notation and we're thinking of it as streaming
additive terms to the log density accumulator.

- Bob

> On Dec 31, 2015, at 7:57 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> I like it, especially getting rid of sampling statements. I would like it even more if the notation was just = for assignment and += for increment and "ld" was just special because it was reserved. Seems overspecified to have a special increment syntax for a special variable... Or is the proposal to allow << to increment anything?
>

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

It will return the log density. Everything is on the log scale.

If you want the density, it's just exp(cauchy(y | mu, sigma)).

We don't currently have a function for the linear density, just
cauchy_log(). So the proposal is to replace "cauchy_log" with
"cauchy". That's the most confusing part of this, but also
the thing that cleans up the notation the most.

- Bob

[Michael Betancourt]

>
> 3.a. Deprecate the use of any function ending in "_log" and
> go with the above form.

I still prefer moving to an _ld suffix to emphasize that these are
log densities, but I won’t fight no suffixal all.

[Bob Carpenter]

I'm curious what other people think about this, as I think it's
the most contentious part of the proposal.

There are three locations (now is what it is now; proposed is what I proposed):

* Stan language usage (now "_log"; proposed remove)

ld << cauchy_ld(y | mu, sigma);

* user-defined functions to be used as densities (now "_log"; proposed keep)

* Stan math lib suffixes (now "_log"; proposed keep)

- Bob

[Krzysztof Sakrejda]

I'm confused what this "keep" adds. This is probably because I lost some other piece of context. If user-defined densities have a <propto=true> or <propto=false>, then they need special treatment of some sort, otherwise they're just something you would otherwise pass to increment_log_prob()... no? I guess they could increment ld as a side effect rather than returning the increment... but does it limit anything to disallow that? That's just a specialized version of increment_log_prob()...

[Michael Betancourt]

The _log _anywhere_ is misleading and I’d rather
work on deprecating it everywhere if possible,
especially if we’re keeping the GLM-link-function
as-suffix convention.

[Krzysztof Sakrejda]

This is probably the last time that Stan gets to rearrange the API and get away with it. Might as well go to town...

[Avraham Adler]

Speaking as the uninformed peanut gallery here, I agree with Michael that it will cause less confusion if the densities are explicitly denoted as being on the log scale.

Avi

[Jeffrey Arnold]

I agree that _log is misleading and would prefer _ld. Although these proposals are a big change, it's nice to remove the BUGS "original sin".

[Joshua N Pritikin]

On Thursday, December 31, 2015 at 12:57:10 AM UTC+2, Bob Carpenter wrote:
> I want to run this by everyone on stan-dev before
> creating the issue and starting to work on it.
>
> 1. New density expression with normalization and vertical bar:
>
> cauchy<norm=true>(y | mu, sigma)

I know the vertical bar signifies conditional probability, but is the
vertical bar anything more than part of the name of the density?

In other words, the vertical bar can only appear at one position in the
argument list, correct?

This seems a bit idiosyncratic from the point of view of computer
language design. Typically, syntactic elements are introduced to
signify some particular semantic. Can I write cauchy(y, mu | sigma) or
cauchy(y, sigma | mu) or cauchy(y, mu, sigma)? If not, then I wonder
whether the vertical bar serves to clarify or complicate the language.
Does it profit authors to be required to remember where the vertical
bar goes in the argument list?

--
Joshua N. Pritikin
Department of Psychology
University of Virginia
485 McCormick Rd, Gilmer Hall Room 102
Charlottesville, VA 22904
http://people.virginia.edu/~jnp3bc

[Krzysztof Sakrejda]

I agree about the vertical bar, seems mostly like another way to make syntax errors.

[Bob Carpenter]

> On Dec 31, 2015, at 2:58 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> On Thursday, December 31, 2015 at 2:47:10 PM UTC-5, Bob Carpenter wrote:
>>> On Dec 31, 2015, at 2:29 PM, Michael Betancourt <betan...@gmail.com> wrote:
>>>
>>>>
>>>> 3.a. Deprecate the use of any function ending in "_log" and
>>>> go with the above form.
>>>
>>> I still prefer moving to an _ld suffix to emphasize that these are
>>> log densities, but I won’t fight no suffixal all.
>>
>> I'm curious what other people think about this, as I think it's
>> the most contentious part of the proposal.
>>
>> There are three locations (now is what it is now; proposed is what I proposed):
>>
>> * Stan language usage (now "_log"; proposed remove)
>>
>> ld << cauchy_ld(y | mu, sigma);
>>
>> * user-defined functions to be used as densities (now "_log"; proposed keep)
>>

>> * Stan math lib suffixes (now "_log"; proposed keep)
>

> I'm confused what this "keep" adds.

Right now, if you define a function "foo" with a "_log" suffix,
as in "foo_log", as we have done for our density functions,
then we allow

y ~ foo(theta);

to be equal (mod dropped constants) to

increment_log_prob(foo_log(y, theta));

If we change that, we break all the user-defined functions. But
we're going to break ~ anyway (or at least deprecate it).

> This is probably because I lost some other piece of context. If user-defined densities have a <propto=true> or <propto=false>, then they need special treatment of some sort, otherwise they're just something you would otherwise pass to increment_log_prob()... no? I guess they could increment ld as a side effect rather than returning the increment... but does it limit anything to disallow that? That's just a specialized version of increment_log_prob()...

increment_log_prob() is going away --- it'll be deprecated
or eliminated in Stan 3.

- Bob

>>
>> - Bob

[Bob Carpenter]

Do you think

ld << normal_ld(y | mu, sigma);

is a strong enough indicate that "normal_ld" is a log density?

- Bob

[Bob Carpenter]

> On Dec 31, 2015, at 9:50 PM, Joshua N Pritikin <jpri...@pobox.com> wrote:
>
> On Thursday, December 31, 2015 at 12:57:10 AM UTC+2, Bob Carpenter wrote:
>> I want to run this by everyone on stan-dev before
>> creating the issue and starting to work on it.
>>
>> 1. New density expression with normalization and vertical bar:
>>
>> cauchy<norm=true>(y | mu, sigma)
>
> I know the vertical bar signifies conditional probability, but is the
> vertical bar anything more than part of the name of the density?
>
> In other words, the vertical bar can only appear at one position in the
> argument list, correct?
>
> This seems a bit idiosyncratic from the point of view of computer
> language design. Typically, syntactic elements are introduced to
> signify some particular semantic. Can I write cauchy(y, mu | sigma) or
> cauchy(y, sigma | mu) or cauchy(y, mu, sigma)? If not, then I wonder
> whether the vertical bar serves to clarify or complicate the language.
> Does it profit authors to be required to remember where the vertical
> bar goes in the argument list?

We'll be using it the conventional way to separate the variate/outcome
(left of bar) from parameters (right of bar). So you'd only get
normal(y | mu, sigma), not any other bar placement. And in fact
for all of our current densities, there's only argument to the left. If
we had a bivariate distribtuion that we wrote p(x, y | theta), then we'd
have a case for putting two things on the left.

It's just the way (Bayesian) statisticians write these functions. And don't
worry, Andrew, we're not going to write normal(y; mu, sigma) as an alternative
for models intended for max likelihood estimation.

- Bob

>
> --
> Joshua N. Pritikin
> Department of Psychology
> University of Virginia
> 485 McCormick Rd, Gilmer Hall Room 102
> Charlottesville, VA 22904
> http://people.virginia.edu/~jnp3bc
>

[Bob Carpenter]

I'm not wedded to the vertical bar. Everyone seemed to think
it was a good idea in the Stan meeting.

- Bob

> On Jan 1, 2016, at 9:39 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> I agree about the vertical bar, seems mostly like another way to make syntax errors.
>

[Bob Carpenter]

> On Dec 31, 2015, at 11:20 AM, Sebastian Weber <sdw....@gmail.com> wrote:
>
> Hi Bob!
>
> This sounds like a huge thing - is it worth investing this effort?

The biggest effort's going to be in rewriting all the doc. And
the models, but we already want to do both of those things anyway,
so now's the time to change the language. We've been discussing this
for a while and we've been holding off on updating models for this
reason. It's terrible to be frozen waiting for a big change.

> Also, the deprecation period is hopefully very long,

The deprecation may be short --- as in introducing it in 2.10 and then dropping
it some time next year when we bring out Stan 3. But so far, we don't know.
My worry is that people will keep using the old notation.

The bigger worry is that any book, etc., with a Stan model will be
out of date as soon as we do this. My impression is that I'm the one
most concerned with breaking backward compatibility, but I may be wrong.

Stan 3's the time to change things if at all (though it doesn't necessarily
mean we have to break backward compatibility).

> i.e. Stan programs will look very different wrt to BUGS programs with this.

That's a feature, not a bug (see below).

> My comments:
>
> 1. Why introduce another logic with "<norm=true/false>"? Do you plan to introduce further options here, then it makes sense. If only for discriminating between normed/unnormed I would prefer a _norm naming convention thing.

Right now,

y ~ foo(theta);

drops normalizing constants but

increment_log_prob(foo_log(y, theta));

keeps them. We'd like to be able to control the behavior on
a case-by-case behavior.

The confusing bit's still going to be that our densities are
on the unconstrained scale with Jacobians, even if they get
the normalizing constants.

> 2. What is the rationale to introduce two new operators "<<" and ","?

The rationale was to make it clearer that everything in
Stan is just incrementing a log density. So instead of having
both sampling statements and increment_log_prob statements, we
just have the log density accumulator.

The problem we have now is that people interpret the "~" as random
number generation or as graphical model building and sometimes just
copying BUGS and hoping things will work. Instead, we want to emphasize
that a Stan program is just a regular imperative program defining
a log density.

- Bob

> Maybe I haven´t understood fully, but writing
>
> lp << density1 << density2;
>
> makes more sense to me.

The notation

lp << a, b, c, d;

is from Eigen. C++ output streams would be

lp << a << b << c;

which I don't think anyone's going to like.

> In any case, I suggest to try out this new scheme on people who are more statisticians than programmers,

Andrew's the one who's been lobbying most strongly for it!

> i.e. I think someone having a programming background will absorb this easily, but I am not so sure about those people without a programming background.

Stan's a programming language, so I'm not so worried about people
who don't think they're doing programming or haven't done it before.

- Bob

[Ben Goodrich]

On Friday, January 1, 2016 at 3:25:54 PM UTC-5, Bob Carpenter wrote:

I'm not wedded to the vertical bar.  Everyone seemed to think
it was a good idea in the Stan meeting.

I think we should do the vertical bar even if users don't like it. It is worth it to reinforce correct statistical thinking.

Ben

[Michael Betancourt]

Yes vertical bar.

Yes _ld suffix to language and math.

[Andrew Gelman]

+1

[Krzysztof Sakrejda]

On Friday, January 1, 2016 at 3:19:17 PM UTC-5, Bob Carpenter wrote:
> Do you think
>
> ld << normal_ld(y | mu, sigma);
>
> is a strong enough indicate that "normal_ld" is a log density?

Not by itself, but it's less likely to get confused for a link specification than _log. It will also give users something to search the manual for. Is there a Stan 3 longer-term-manual issue? The first mention of "_ld" in the manual should explain the convention, even it's just a text box rather than inline.

[Daniel Lee]

Vertical bar:
- will there be any ambiguity with the or operator?
- will expressions not be allowed on the left of the bar? In which case, the first question is moot.
- will expressions be allowed for all the arguments? This seems natural.

_ld:
+1

While we're at it, perhaps we should think about:
- a syntax for writing custom derivatives for functions
- a testing framework written in the language

[Krzysztof Sakrejda]

I tried searching the list for the background about the vertical bar/pipe discussion but it's not coming up, sorry if this covers well-trod ground.

On Saturday, January 2, 2016 at 9:25:28 AM UTC-5, Daniel Lee wrote:
> Vertical bar:
> - will there be any ambiguity with the or operator?

I've been trying to figure out how to express my objection and I think it's related to this. There already is a well-established separator for function arguments in programming and for better or worse it's a comma (almost?) everywhere. I don't think that programming syntax is the place to correct statistical thinking when it's going to cause syntax errors and mailing list traffic.

I'm trying to come up with another example where a programming language uses a syntactic choice to make a conceptual point that can't be inferred without consulting an outside source.... dplyr (well, magrittr) did this sort of thing in R by introducing a <- b %>% c as a way of stringing together arbitrary operations but a) that matches the way people (at least the "data science" people I know) think about data operations, it doesn't try to change it b) it builds on existing streaming concepts; and c) it's still screwy because the final assignment is on the ultimate LHS.

In any multi-level setting you can't even think about the LHS of the bar as an observable so is the goal here to get people to think/admit that they're specifying a model in terms of conditional distributions? Does forcing thinking in terms of conditional probabilities still make sense if specifying the model requires extra computations that increment the log probability through other functions? Is there a specific user mistake we're trying to avoid?

I get that most of the developers are on board with this notation, I just don't get it given the costs (mostly in terms of mailing list traffic where the answer is going to be some polite version of RTFM/RTFAQ).

Krzysztof

[Michael Betancourt]

ALL of the densities we implement are conditional densities. Densities are
particularly subtle objects because they are not just functions with two types
of inputs, the two inputs transform in different ways. This is of course related
to much of the confusion with Jacobians and the like.

[Aki Vehtari]

When discussing the new notation, please, include examples how it will affect
functions such as
normal_cdf
normal_cdf_log
normal_ccdf_log
Why not have a consistent approach, e.g.
normal_pdf (or normal_d?)
normal_lpdf (or normal_ld?)
normal_cdf (or ?)
normal_lcdf
normal_ccdf
normal_lccdf?

Also if optional arguments are allowed as normal_lpdf<norm=true>, why not then have also link functions as optional arguments, e.g.,
bernoulli_lpdf<link=probit> or bernoulli_lpdf<link=logit> or bernoulli_lpdf<link=linear> ?

Aki

[Bob Carpenter]

> On Jan 2, 2016, at 12:06 PM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> When discussing the new notation, please, include examples how it will affect
> functions such as
> normal_cdf
> normal_cdf_log
> normal_ccdf_log
> Why not have a consistent approach, e.g.
> normal_pdf (or normal_d?)
> normal_lpdf (or normal_ld?)
> normal_cdf (or ?)
> normal_lcdf
> normal_ccdf
> normal_lccdf?

That makes a lot of sense to me. And I prefer the suffix to
another <...> argument.

We also need things like differences of CDFs for full
truncation:

normal_diff_lcdf(y_low, y_high, mu, sigma)

=def= log(normal_cdf(y_high,mu,sigma) - normal_cdf(y_low,mu,sigma))

Right now, it's being implemented as above, but it could really
use better implementations. Then we could think about vectorizing
the truncated densities.

> Also if optional arguments are allowed as normal_lpdf<norm=true>, why not then have also link functions as optional arguments, e.g.,
> bernoulli_lpdf<link=probit> or bernoulli_lpdf<link=logit> or bernoulli_lpdf<link=linear> ?

We could. I don't want to get carried away and in particular don't
want expressions as values of those = things --- just template
parameters that resolve which function is called.

- Bob

[Bob Carpenter]

Oh, and the reason I don't like normal_pdf and that
I do like the vertical bar is that it is what you see
in stats papers for the function:

cauchy(y | mu, tau)

But on the log scale, we don't want to write

log cauchy(y | mu, tau) or log(cauchy(y | mu, tau))

because it's two functions (though I suppose I could fish it
out of the abstract syntax tree and reduce it during code
generation).

- Bob

> On Jan 2, 2016, at 12:06 PM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>

[Bob Carpenter]

That won't be part of the Stan 3 manual issue, but part of
changing the language issue.

The manual may go through a transition if I don't have
time to recode all the models and remove all mentions of
the deprecated operations.

- Bob

> On Jan 2, 2016, at 8:00 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>

[Bob Carpenter]

The vertical bar signals the normalization, so that

foo(a | b, c)

means that

INTEGRATE_A foo(a | b, c) d.a = 1

for any legal b, c, where A is the support (sometimes given b and c
as in uniform and pareto).

- Bob

[Bob Carpenter]

> On Jan 2, 2016, at 11:13 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> I tried searching the list for the background about the vertical bar/pipe discussion but it's not coming up, sorry if this covers well-trod ground.
>
> On Saturday, January 2, 2016 at 9:25:28 AM UTC-5, Daniel Lee wrote:
>> Vertical bar:
>> - will there be any ambiguity with the or operator?
>
> I've been trying to figure out how to express my objection and I think it's related to this. There already is a well-established separator for function arguments in programming and for better or worse it's a comma (almost?) everywhere. I don't think that programming syntax is the place to correct statistical thinking when it's going to cause syntax errors and mailing list traffic.
>
> I'm trying to come up with another example where a programming language uses a syntactic choice to make a conceptual point that can't be inferred without consulting an outside source....

I'm not sure what you mean by "inferred" here. Pretty much
everything in a programming language syntax needs to be looked
up.

We could also allow "," in this context, but I think that'd be
even more confusing.

R uses %*% for matrix multiplication, and you certainly need
to look that up, but I don't think it's trying to make a conceptual
point (other than that matrix multiplication is a different function
and they'd already decided that plain old * acted elementwise).

C/C++/Java use | vs. || and & vs. && to differentiate bitwise
vs. logical disjunction and conjuncation. But I think the reason
there is that they give different answers for integer inputs.

C++ uses a -> b to resolve a pointer, but a.b to resolve a reference.
All need to be looked up before using the language.

> In any multi-level setting you can't even think about the LHS of the bar as an observable

That's not the point; if we wanted to do that, we'd use the semicolon :-)

It's that the | indicates conditioning --- hence its effect on normalization.

> so is the goal here to get people to think/admit that they're specifying a model in terms of conditional distributions?

No, just to make it easier to read models!

> Does forcing thinking in terms of conditional probabilities still make sense if specifying the model requires extra computations that increment the log probability through other functions? Is there a specific user mistake we're trying to avoid?

We certainly don't force anyone to think in terms of conditional densities.
You can write everything out via increment-log-prob if you want.

> I get that most of the developers are on board with this notation, I just don't get it given the costs (mostly in terms of mailing list traffic where the answer is going to be some polite version of RTFM/RTFAQ).

I'm guessing this isn't going to be a big deal on the mailing list.
We'll only get our answer if we introduce the bar!

- Bob

[Bob Carpenter]

> On Jan 2, 2016, at 9:25 AM, Daniel Lee <bea...@alum.mit.edu> wrote:
>
> Vertical bar:
> - will there be any ambiguity with the or operator?
> - will expressions not be allowed on the left of the bar? In which case, the first question is moot.
> - will expressions be allowed for all the arguments? This seems natural.

Luckily, our "or" operator is "||", so it won't be a problem.
We don't have the bitwise or-operator "|". And all the parsing
is left-to-right, top-down, and greedy, so

foo(a || b | c)

will get parsed as

foo((a || b) | c)


> _ld:
> +1

Looks like everyone's in agreement (other than me), so we
can go with the wisdom of the crowd :-)

> While we're at it, perhaps we should think about:
> - a syntax for writing custom derivatives for functions
> - a testing framework written in the language

OK --- could you start a new thread, though, and not hijack this one?

- Bob

[Andrew Gelman]

Left-to-right, top down, and greedy. That pretty much describes the entire Stan attitude!

[Aki Vehtari]

On Saturday, January 2, 2016 at 7:35:50 PM UTC+2, Bob Carpenter wrote:
> We also need things like differences of CDFs for full
> truncation:
>
> normal_diff_lcdf(y_low, y_high, mu, sigma)

normal_ldcdf?
(it's log of differences and not diff of logs)

> > Also if optional arguments are allowed as normal_lpdf<norm=true>, why not then have also link functions as optional arguments, e.g.,
> > bernoulli_lpdf<link=probit> or bernoulli_lpdf<link=logit> or bernoulli_lpdf<link=linear> ?
>
> We could. I don't want to get carried away and in particular don't
> want expressions as values of those = things --- just template
> parameters that resolve which function is called.

One reason for this would be to avoid confusion of what _log in the end of function stands for (log of density or log link).

Aki

[Aki Vehtari]

On Sunday, January 3, 2016 at 12:45:44 AM UTC+2, Bob Carpenter wrote:
> > _ld:
> > +1
>
> Looks like everyone's in agreement (other than me), so we
> can go with the wisdom of the crowd :-)

I didn't agree as _ld is not consistent with _cdf and _ccdf, unless they are changed to _cd and _ccd?

Aki

[Aki Vehtari]

On Saturday, January 2, 2016 at 7:43:10 PM UTC+2, Bob Carpenter wrote:
> Oh, and the reason I don't like normal_pdf and that
> I do like the vertical bar is that it is what you see
> in stats papers for the function:
>
> cauchy(y | mu, tau)

The stas paper also use ~, which has been considered problematic in Stan.

If we have
<< normal(y | mu, sigma)
do we then have, e.g. (toy example)
<< normal(y | mu, sigma), normal_cdf(y_high,mu,sigma) - normal_cdf(y_low,mu,sigma)
or
<< normal(y | mu, sigma), log(normal_cdf(y_high,mu,sigma) - normal_cdf(y_low,mu,sigma))
another example
<< normal(y | mu, sigma), 1/sigma
or
<< normal(y | mu, sigma), -log(sigma)


I would love to use notation close to stats papers, but I don't see yet how we obtain consistency when combining built-in distributions and other expressions when using ~ or <<

Aki

[Bob Carpenter]

I think everyone is OK with:

* densities need some suffix

* it should not be "_log"

* log-scale should be explicitly indicated

I don't think there's a consensus on what the suffixes should
be, but Aki's right that we should look at the whole family.
I think his proposal makes a lot of sense:

pdf / cdf / ccdf / dcdf
lpdf / lcdf / lccdf / ldcdf

I somehow think the "f" is redundant --- I assume we can at
least assume our users know what a function is --- but then it's in
common usage.

Is everyone OK with that? And with using "pdf" for both density
and mass functions? (Yes, they're all densities with the right
math, however the textbooks say "pmf").

What I don't think there's consensus on is

* using a vertical bar, as in:

normal_lpdf(y | mu, sigma)

Also, I don't think there's consensus on the whole "<<" business,
as in:

ld << ...; // Michael
<< ...; // Aki
ld += ...; // Krzysztof

Krzysztof thinks the vertical bar and << will confuse non-programmers.

And should that "ld" match "pdf" or vice-versa?

I think some people are thinking we might add other accumulators like
lprior or llikelihood (there's that pesky "l" for "log" again), or
maybe "lbase" and "ldelta".

- Bob

[Michael Betancourt]

Aaaaaaah. Do we need official minutes? Almost all of these
issues were introduced in the last discussion of the proposed
changes.

> I don't think there's a consensus on what the suffixes should
> be, but Aki's right that we should look at the whole family.
> I think his proposal makes a lot of sense:
>
> pdf / cdf / ccdf / dcdf
> lpdf / lcdf / lccdf / ldcdf
>
> I somehow think the "f" is redundant --- I assume we can at
> least assume our users know what a function is --- but then it's in
> common usage.
>
> Is everyone OK with that? And with using "pdf" for both density
> and mass functions? (Yes, they're all densities with the right
> math, however the textbooks say "pmf”).

The problem with differentiating between pdfs and pmfs is that
people will get confused about whether the posterior is a pdf or
a pmf. It’s always a pdf for the models you can build in Stan,
but it’s not trivial to grok what goes where (pdfs over parameters,
hence prior or posterior, pdfs over data, pmfs over data). Then
again, being more explicit may provide an opportunity for more
clear error messages, such as when a user tries to use a pmf
over parameters.

Things become even more complicated when we use CDFs,
because they are fundamentally different from densities.

If we want to be more explicit and we don’t mind the extra typing
then I’d say go with

pdf / cdf / ccdf / dcdf
lpdf / lcdf / lccdf / ldcdf

with the understanding that not all may be implemented.

The f is technically redundant, but if we drop it then we’re
comparing pd and cd which are nonstandard and confusing
given the dual use of “d”.

> Krzysztof thinks the vertical bar and << will confuse non-programmers.

If anything “|” will confuse more programmers than users. And
“<<“ is no less intuitive for non-programmers than “+=“.

> And should that "ld" match "pdf" or vice-versa?

If we go with a more explicit choice then yes, it should become lpdf.

[Aki Vehtari]

On Sunday, January 3, 2016 at 2:18:39 AM UTC+2, Bob Carpenter wrote:
> Also, I don't think there's consensus on the whole "<<" business,
> as in:
>
> ld << ...; // Michael
> << ...; // Aki
> ld += ...; // Krzysztof

I'm fine with ld << ...;

If we want stats paper look, e.g.

ld << cauchy(tau | 0, 2.5),
normal(beta | 0, 2),
bernoulli_logit(y | x * (tau * beta));

then here << is better than +=, because += would imply summing distributions (not log densities). So the corresponding version with += would look like

ld += cauchy_lpdf(tau | 0, 2.5),
normal_lpdf(beta | 0, 2),
bernoulli_logit_lpdf(y | x * (tau * beta));

this latter form is more explicit and it migh be more natural if other than built-in distributions are used, e.g.

ld += normal_lpdf(y | mu, sigma) -log(sigma)

but maybe we could use << as

ld << normal(y | mu, sigma), 1/sigma

so that _lpdf's will be used if they are available and otherwise log of the expression is used. Yes, there is possibility of confusion by using both << and +=, but it would give possiblity for both stat paper look and more explicit log density construction.

Aki

[Michael Betancourt]

No, no, no. No automated processing, especially conditional processing. That is ripe for confusion. A stream operator is nice because it gives the notion of "adding" a new component to the target without specifying the representation. Yes, it will always be log densities for us, but it's nice to abstract the concept to reinforce pedagogy.

When moving beyond standard distributions the user is on their own in ensuring validity. Any processing on our part can give the user a false sense of security here which we definitely don't want.

[Krzysztof Sakrejda]

On Sunday, January 3, 2016 at 12:07:39 AM UTC-5, Michael Betancourt wrote:
> Aaaaaaah. Do we need official minutes? Almost all of these
> issues were introduced in the last discussion of the proposed
> changes.

I don't know about minutes but maybe wiki-level documentation of major design decisions would help. There's going to be some friction to having this many people involved.

>
> > I don't think there's a consensus on what the suffixes should
> > be, but Aki's right that we should look at the whole family.
> > I think his proposal makes a lot of sense:
> >
> > pdf / cdf / ccdf / dcdf
> > lpdf / lcdf / lccdf / ldcdf
> >
> > I somehow think the "f" is redundant --- I assume we can at
> > least assume our users know what a function is --- but then it's in
> > common usage.
> >
> > Is everyone OK with that? And with using "pdf" for both density
> > and mass functions? (Yes, they're all densities with the right
> > math, however the textbooks say "pmf”).
>
> The problem with differentiating between pdfs and pmfs is that
> people will get confused about whether the posterior is a pdf or
> a pmf. It’s always a pdf for the models you can build in Stan,
> but it’s not trivial to grok what goes where (pdfs over parameters,
> hence prior or posterior, pdfs over data, pmfs over data). Then
> again, being more explicit may provide an opportunity for more
> clear error messages, such as when a user tries to use a pmf
> over parameters.

I agree that being explicit about this material is a good opportunity to
make Stan usage clearer.

> Things become even more complicated when we use CDFs,
> because they are fundamentally different from densities.
>
> If we want to be more explicit and we don’t mind the extra typing
> then I’d say go with
>
> pdf / cdf / ccdf / dcdf
> lpdf / lcdf / lccdf / ldcdf

I'm on board with this. Also on board with differentiating pdf/pmf where relevant but I don't have a strong opinion on it.

[Krzysztof Sakrejda]

On Saturday, January 2, 2016 at 7:18:39 PM UTC-5, Bob Carpenter wrote:
> I think everyone is OK with:

Thanks for staying on top of doing the summaries Bob, at the moment it's the
only way I keep this discussion straight.

>
> * densities need some suffix
>
> * it should not be "_log"
>
> * log-scale should be explicitly indicated
>
> I don't think there's a consensus on what the suffixes should
> be, but Aki's right that we should look at the whole family.
> I think his proposal makes a lot of sense:
>
> pdf / cdf / ccdf / dcdf
> lpdf / lcdf / lccdf / ldcdf
>
> I somehow think the "f" is redundant --- I assume we can at
> least assume our users know what a function is --- but then it's in
> common usage.
>
> Is everyone OK with that? And with using "pdf" for both density
> and mass functions? (Yes, they're all densities with the right
> math, however the textbooks say "pmf").
>
> What I don't think there's consensus on is
>
> * using a vertical bar, as in:
>
> normal_lpdf(y | mu, sigma)
>
> Also, I don't think there's consensus on the whole "<<" business,
> as in:
>
> ld << ...; // Michael
> << ...; // Aki
> ld += ...; // Krzysztof

My though was more specific---I think we should stick to the more common operators:

ld += A + B + C;

where A/B/C can be any function that evaluates to something compatible with the ld accumulator. It lets us use a construct that appears in multiple other programming languages and avoids introducing a comma addition. I think the comma is really too clever.

>
> Krzysztof thinks the vertical bar and << will confuse non-programmers.

I think I figured out from Michael's comments that the bar is adding a piece of information about the function not available anywhere else (except the manual). I still don't understand if users are making a modeling mistake that this distinction helps to resolve but I'm ok with not knowing that.

> And should that "ld" match "pdf" or vice-versa?
>
> I think some people are thinking we might add other accumulators like
> lprior or llikelihood (there's that pesky "l" for "log" again), or
> maybe "lbase" and "ldelta".

> - Bob

[Aki Vehtari]

On Sunday, January 3, 2016 at 11:56:09 AM UTC+2, Aki Vehtari wrote:
> ld += cauchy_lpdf(tau | 0, 2.5),
> normal_lpdf(beta | 0, 2),
> bernoulli_logit_lpdf(y | x * (tau * beta));

Here I meant to have plus instead of comma if += is used, ie
ld += cauchy_lpdf(tau | 0, 2.5) +
normal_lpdf(beta | 0, 2) +

bernoulli_logit_lpdf(y | x * (tau * beta));

Aki

[Aki Vehtari]

On Sunday, January 3, 2016 at 12:29:46 PM UTC+2, Michael Betancourt wrote:
> No, no, no. No automated processing, especially conditional processing.

Constraints do automated processing, too :)

How about having
normal
to be alias for
normal_lpdf
?

Or do we want to strictly use _lpdf even if it makes the code look less like stat papers?

Aki

[Bob Carpenter]

> On Jan 3, 2016, at 8:27 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
> ...

> My though was more specific---I think we should stick to the more common operators:
>
> ld += A + B + C;
>
> where A/B/C can be any function that evaluates to something compatible with the ld accumulator. It lets us use a construct that appears in multiple other programming languages and avoids introducing a comma addition. I think the comma is really too clever.

I agree that we want the right-hand side to be any old expression.
So whatever is used for a density in place of A/B/C will just
be added to the log density.

+= is common, just not in R or BUGS or JAGS, so for
this particular audience, not so common.

I'm OK with not having a comma. You can think of it
as streaming in a vector, which is what we do now with

vector[N] x;
...
increment_log_prob(x)

which is equivalent to

for (n in 1:N)
increment_log_prob(x[n]);

So I guess another question is whether we allow vectors
on the right-hand-side of <<.

>> Krzysztof thinks the vertical bar and << will confuse non-programmers.
>
> I think I figured out from Michael's comments that the bar is adding a piece of information about the function not available anywhere else (except the manual). I still don't understand if users are making a modeling mistake that this distinction helps to resolve but I'm ok with not knowing that.

No modeling mistake, just the convention p(y | theta)
for conditional densities as opposed to p(y, theta) for
joint densities. And normal(y | mu, sigma) is the
conditional of y given mu and sigma in the sense that
for fixed sigma > 0 and any mu, it integrates to 1 over y.

- Bob

[Bob Carpenter]

> On Jan 3, 2016, at 4:56 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> On Sunday, January 3, 2016 at 2:18:39 AM UTC+2, Bob Carpenter wrote:
>> Also, I don't think there's consensus on the whole "<<" business,
>> as in:
>>
>> ld << ...; // Michael
>> << ...; // Aki
>> ld += ...; // Krzysztof
>
> I'm fine with ld << ...;
>
> If we want stats paper look, e.g.
>
> ld << cauchy(tau | 0, 2.5),
> normal(beta | 0, 2),
> bernoulli_logit(y | x * (tau * beta));
>
> then here << is better than +=, because += would imply summing distributions (not log densities). So the corresponding version with += would look like
>
> ld += cauchy_lpdf(tau | 0, 2.5),
> normal_lpdf(beta | 0, 2),
> bernoulli_logit_lpdf(y | x * (tau * beta));
>
> this latter form is more explicit and it migh be more natural if other than built-in distributions are used, e.g.
>
> ld += normal_lpdf(y | mu, sigma) -log(sigma)
>
> but maybe we could use << as
>
> ld << normal(y | mu, sigma), 1/sigma
>
> so that _lpdf's will be used if they are available and otherwise log of the expression is used. Yes, there is possibility of confusion by using both << and +=, but it would give possiblity for both stat paper look and more explicit log density construction.

I don't want that much context sensitivity. I really want
normal(y | mu, sigma) to be assignable, so that if
we allow

ld << normal(y | mu, sigma)

then we'd also allow expressions like

sum <- sum + normal(y[n] | mu, sigma);

and the normal would be on the log scale. But everyone's
convinced me that's too confusing and we'll need a suffix
to sort out which function associated with the normal distribution
that we mean.

Also, we want to control the way the log density is coded
and taking logs would leave too much up to the compiler.

- Bob

[Bob Carpenter]

> On Jan 3, 2016, at 9:09 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> On Sunday, January 3, 2016 at 12:29:46 PM UTC+2, Michael Betancourt wrote:
>> No, no, no. No automated processing, especially conditional processing.
>
> Constraints do automated processing, too :)
>
> How about having
> normal
> to be alias for
> normal_lpdf
> ?

That's really been my thinking all along.

But I think it's going to be confusing because
everyone's going to expect normal(y | mu, sigma)
to be on the linear scale because that's the convention.

> Or do we want to strictly use _lpdf even if it makes the code look less like stat papers?

I think we're going to have to.

- Bob

[Michael Betancourt]

Honestly? Most stat papers are terribly written so we shouldn't necessarily use them as a baseline. People (including me initially) might complain that lpdf is unnecessary, but it makes the model specification very clear. And we can tie it into a rich Stats intro in the manual to help with pedagogy, too.

[Bob Carpenter]

> On Jan 3, 2016, at 12:07 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Aaaaaaah. Do we need official minutes? Almost all of these
> issues were introduced in the last discussion of the proposed
> changes.

Right --- I was just trying to summarize where we were
at in terms of consensus vs. decisions still to be made.

I added a summary here:

https://github.com/stan-dev/stan/wiki/Stan-3-Density-Notation-and-Increments

Feel free to add or edit. I thought I'd try to recall my
own understanding to make sure everyone was on the same page.

We need to get this right this time --- we can't be mucking
with our fundamental notation more than once. Even with this,
there may be a user revolt.

>> I don't think there's a consensus on what the suffixes should
>> be, but Aki's right that we should look at the whole family.
>> I think his proposal makes a lot of sense:
>>
>> pdf / cdf / ccdf / dcdf
>> lpdf / lcdf / lccdf / ldcdf
>>
>> I somehow think the "f" is redundant --- I assume we can at
>> least assume our users know what a function is --- but then it's in
>> common usage.
>>
>> Is everyone OK with that? And with using "pdf" for both density
>> and mass functions? (Yes, they're all densities with the right
>> math, however the textbooks say "pmf”).
>
> The problem with differentiating between pdfs and pmfs is that
> people will get confused about whether the posterior is a pdf or
> a pmf. It’s always a pdf for the models you can build in Stan,

That's true of the unconstrained density used for sampling.
But the generated quantities can be discrete.

> but it’s not trivial to grok what goes where (pdfs over parameters,
> hence prior or posterior, pdfs over data, pmfs over data). Then
> again, being more explicit may provide an opportunity for more
> clear error messages, such as when a user tries to use a pmf
> over parameters.
>
> Things become even more complicated when we use CDFs,
> because they are fundamentally different from densities.

At least there's no PDF vs PMF distinction.

> If we want to be more explicit and we don’t mind the extra typing
> then I’d say go with
>
> pdf / cdf / ccdf / dcdf
> lpdf / lcdf / lccdf / ldcdf
>
> with the understanding that not all may be implemented.

> ...

>
>> Krzysztof thinks the vertical bar and << will confuse non-programmers.
>
> If anything “|” will confuse more programmers than users. And
> “<<“ is no less intuitive for non-programmers than “+=“.

Intuitive is all a matter of perspective.

I like "+=" if we're thinking of the "ld" or "lpdf" as a variable.
Too bad we don't use "=" for assignment.

I'm still worried about

lpdf << normal_lpdf(y | mu, sigma);

I can't tell if the duplication of "_lpdf" is an affordance or
a hindrance.

And it's so much longer than

<< normal(y | mu, sigma);

that I worry it'll be hard to read models.

>> And should that "ld" match "pdf" or vice-versa?
>
> If we go with a more explicit choice then yes, it should become lpdf.

OK --- I didn't quite get that in the Wiki page. I just used ld. There's
still the issue of maybe just leaving it blank.

- Bob

[Michael Betancourt]

I’m not a huge fan of keeping it blank.  I’m guessing the

appropriate “target_lpdf” would be too much to type, but

could we get away with “t_lpdf”?  I’d be fine with making

this “p_lpdf” or “post_lpdf” if we want to emphasize Bayes.

[Aki Vehtari]

On Monday, January 4, 2016 at 10:53:47 AM UTC+2, Bob Carpenter wrote:
> I added a summary here:
> https://github.com/stan-dev/stan/wiki/Stan-3-Density-Notation-and-Increments

Nice.

On Monday, January 4, 2016 at 9:28:17 AM UTC+2, Bob Carpenter wrote:
> everyone's going to expect normal(y | mu, sigma)
> to be on the linear scale because that's the convention.

The usual convention is also to return normalised density, but Stan returns unnormalised, which may add confusion (and has caused confusion as seen sometimes in Stan users list).

> cauchy_lpdf<norm=true>(y | mu, tau); // normalized
> cauchy_lpdf<norm=false>(y | mu, tau); // unnormalized
> cauchy_lpdf(y | mu, tau); // unnormalized

I just had an idea that maybe we could use BDA notation, ie, q for unnormalised density
cauchy_lpdf(y | mu, tau); // normalized
cauchy_lqdf(y | mu, tau); // unnormalized

And if we want less letters in all functions, we could use
_p, _lp, _q, _lq, _c, _lc, _cc, _lcc, _dc, _ldc, _r

and then we would have
cauchy_lp(y | mu, tau); // log, normalized
cauchy_lq(y | mu, tau); // log, unnormalized

or if we want to emphasize
cauchy_l_p(y | mu, tau); // log, normalized
cauchy_l_q(y | mu, tau); // log, unnormalized

I think these forms with df dropped, look nice with |

Aki

[Bob Carpenter]

Interesting idea. Still not sure what I make of it, though.
More inline.

> On Jan 4, 2016, at 7:11 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> On Monday, January 4, 2016 at 10:53:47 AM UTC+2, Bob Carpenter wrote:
>> I added a summary here:
>> https://github.com/stan-dev/stan/wiki/Stan-3-Density-Notation-and-Increments
>
> Nice.
>
> On Monday, January 4, 2016 at 9:28:17 AM UTC+2, Bob Carpenter wrote:
>> everyone's going to expect normal(y | mu, sigma)
>> to be on the linear scale because that's the convention.
>
> The usual convention is also to return normalised density, but Stan returns unnormalised, which may add confusion (and has caused confusion as seen sometimes in Stan users list).

Good point. I thought norm=true is a better default, but it's
going to be a killer for model efficiency and since we almost
always want norm=false for just sampling or optimziation, it's
going to be very verbose.

>> cauchy_lpdf<norm=true>(y | mu, tau); // normalized
>> cauchy_lpdf<norm=false>(y | mu, tau); // unnormalized
>> cauchy_lpdf(y | mu, tau); // unnormalized
>
> I just had an idea that maybe we could use BDA notation, ie, q for unnormalised density
> cauchy_lpdf(y | mu, tau); // normalized
> cauchy_lqdf(y | mu, tau); // unnormalized
>
> And if we want less letters in all functions, we could use
> _p, _lp, _q, _lq, _c, _lc, _cc, _lcc, _dc, _ldc, _r

I think this will cause no end of problems because of the
R / BUGS / JAGS convention:

dnorm : pdf
pnorm : cdf
qnorm : inverse cdf
rnorm : random draw

I really regret using lp__, which I was thinking of
as the p in pdf and pmf, not as a log (event) probability.

> and then we would have
> cauchy_lp(y | mu, tau); // log, normalized
> cauchy_lq(y | mu, tau); // log, unnormalized

I really like that this doesn't involve the verbose <norm=true>.

> or if we want to emphasize
> cauchy_l_p(y | mu, tau); // log, normalized
> cauchy_l_q(y | mu, tau); // log, unnormalized

The convention in the C/C++ library is mixed, between
things like log1p and lgamma, but I don't think we
need the underscore. Nor do we need to follow C++.

- Bob

[Andrew Gelman]

Yes, I agree with Bob that there’s a tension here. In BDA, normal(a|b,c^2) has a clearly defined meaning. In Stan, we’re changing things in two ways, and both are good: we’re moving from variance to sd parameterization (thus, normal(a|b,c)), and we’re working with log densities. In BDA notation: lp(theta|y) = const + lp(theta) + lp(y|theta). So as long as we make it clear that normal(a|b,c) is lp, we’re good. I agree with Bob that it would be good if normal(a|b,c) just represents this lp quantity. This is the cleanest way to go. And in this case I’d prefer not to have a _log suffix, as this would imply that not logging is the default. To be clean, we should never allow a density or probability function in Stan to compute a non-logged density. This just means we need to be able to handle zero density values using -Inf on the log scale. I assume we can do this?
A

> On Jan 4, 2016, at 2:27 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
> I don't want that much context sensitivity. I really want
> normal(y | mu, sigma) to be assignable, so that if
> we allow
>
> ld << normal(y | mu, sigma)
>
> then we'd also allow expressions like
>
> sum <- sum + normal(y[n] | mu, sigma);
>
> and the normal would be on the log scale. But everyone's
> convinced me that's too confusing and we'll need a suffix
> to sort out which function associated with the normal distribution
> that we mean.
>
> Also, we want to control the way the log density is coded
> and taking logs would leave too much up to the compiler.
>
> - Bob
>

> --

[Bob Carpenter]

Yes, we can use -inf on the log scale, but that
causes rejection and is always due to either a misspecified
model or arithmetic issues resulting in what Michael (and now
Stan) calls "divergent" transitions.

How do we make it clear that normal(y | mu, sigma) means
what everyone else writes as "log(normal(y | mu, sigma))"?
I was originally leaning toward having normal() be
a synonym (alias) for for normal_lpdf(), but now I'm
leaning toward it being too confusing if it's not packed
right into the name. They might not know what "lpdf" means
at first glance, but at least it'll indicate that there's
something going on.

Took me a while to unpack lp as log p() and the rest as Bayes's rule:

> lp(theta|y) = const + lp(theta) + lp(y|theta)

I don't see how that's relevant.

We don't have the non-log forms of the densities because they're
so unstable arithmetically. We do have the non-log forms of the CDFs,
but I think that's a mistake in retrospect for the same reason.
The only problem is that we don't have the kind of numerical algorithm
expertise (or perhaps interest) to write stable log CDFs and log CCDFs
with derivatives.

- Bob

[Andrew Gelman]

I think it would be best for everything including cdfs to be in log forms but we can discuss.
A

[Bob Carpenter]

Do you think we should get rid of the non-log forms altogether in
the language and make people apply exp() to the log-form if they
want it on the linear scale?

Everyone agrees CDFs and CCDFs should have arithmetically stable log forms.
(And that we should have inverse CDFs for fun with reparameterizations.)

Right now, our log forms are only as stable as the linear forms because
they just literally apply log(foo_cdf(...)) and log1m(foo_cdf(...)).

We need someone to code the robust versions.

- Bob

[Andrew Gelman]

Yes, my current inclination is to have everything on the log scale, as that makes all explanations and notation simpler! And, as you know, explanations and notation count for a lot!

[Michael Betancourt]

We’re not getting good log CDF implementations anytime soon unless
we hire an army of mathematical programmers.

[Aki Vehtari]

Could we use Boost distribution functions?
They give error rates for CDFs, see, e.g.
http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/dist_ref/dists/nc_t_dist.html

Boost doesn't include log versions, not even for pdf, but I guess most of the time they exponentiate as late as possible, so the code could be reused.

As the mainstream is to have linear scale functions, I think we need to be explicit about log. If we are explicit about log, the having linear version of the function should not be much additional work? They would be useful, e.g., in mixture models and in generated quantities.

Aki

[Bob Carpenter]

I think we just need one good one. We don't need to
get everything done.

- Bob

[Bob Carpenter]

We are leaning heavily on Boost. But they don't do derivatives.

- Bob

[Aki Vehtari]

Did you talk about notation in today's meeting? Missed the meeting again...

Aki

[Aki Vehtari]

On Tuesday, January 5, 2016 at 6:39:06 PM UTC+2, Bob Carpenter wrote:
> We are leaning heavily on Boost. But they don't do derivatives.

But you could copy the algorithms written in C++, and just add derivatives?
Reading the code it seems that have thought about precision of the algorithms used and explain the choices in comments.

Aki

[Andrew Gelman]

Hi, yes, we discussed and we agreed that it made sense to have all densities on the log scale by default, since all the expensive computations are done on the log scale, and it seems simpler for users to just have a single function call for each distribution.

Also, we agreed to move future Stan meetings to noon NY time.
A

> On Jan 5, 2016, at 12:59 PM, Aki Vehtari <Aki.V...@aalto.fi> wrote:
>
> Did you talk about notation in today's meeting? Missed the meeting again...
>
> Aki
>

[Bob Carpenter]

The way Boost is templated, we can autodiff the approximation
to function foo(x). But precision bounds on foo(x) don't tell
us anything about precision bounds on the derivative foo'(x)
resulting from autodiffing the approximation to foo().

The right thing to do is develop a new approximate algorithm for
foo'(x). But Boost hasn't done that. Ditto for inverse CDFs.

For example, differentiating beta_cdf(y, alpha, beta) w.r.t. y
is easy, but it's much trickier w.r.t. alpha and beta.

- Bob

[Aki Vehtari]

On Wednesday, January 6, 2016 at 3:46:01 AM UTC+2, Andrew Gelman wrote:
> Hi, yes, we discussed and we agreed that it made sense to have all densities on the log scale by default, since all the expensive computations are done on the log scale, and it seems simpler for users to just have a single function call for each distribution.

Please, tell more!

How about cdf, ccdf, dcdf and icdf? Would these functions work only on log scale, too? Currently normal_cdf is not on log scale, so the same suffix can't be used for log by default. Suffix _lcdf? Something else?

How about normalisation? Off by default? How about vector output (not summing) useful, e.g., for LOO (and waic)?

If
ld << normal(y|mu,sigma)
returns unnormalised summed log pdf, would you go for
log_lik <- normal<norm=true,vector=true>(y|mu,sigma)
to get normalised log pdfs in a vector?

Aki

[Bob Carpenter]

We didn't really have a quorum for this discussion.

Andrew's opinion during the meeting was:

a) always compute normalization terms for output of log-density;
never compute normalization where not needed in the model block;
continue to include normalization in generated quantities

b) everything is on the log scale and there simply aren't
non-log-scale functions, so they don't need to be marked

Both of these are contentious points.

Daniel suggested doing (a) without too much overhead by just recomputing
the model block once per iteration using double values.

If Andrew gets his way, that'd be:

normal_pdf(y | mu, sigma) : log density

normal_cdf(y ??? mu, sigma) : log cdf

normal_ccdf(y ??? mu, sigma) : log ccdf

normal_inv_cdf(y ??? mu, sigma) : log inverse cdf

This would change the behavior of the existing functions normal_cdf()
in a non-backward-compatible way.

Then I believe Andrew wanted to allow "normal" as an abbreviation
for "normal_pdf".

- Bob

[Andrew Gelman]

I thought that the most commonly-used form would be normal(a|b,c) which would be equivalent to normal_log_pdf.

But if we want to spell it out, I think we have to call it normal_log_pdf and normal_log_cdf.

No way we could use normal_pdf to stand for normal_log_pdf. That’s just too confusing.

A

[Ben Goodrich]

On Wednesday, January 6, 2016 at 4:09:33 PM UTC-5, Bob Carpenter wrote:

  normal_inv_cdf(y ??? mu, sigma) : log inverse cdf

Having an inverse CDF return something logged makes no sense (i.e. it could be negative unlogged for a normal distribution). Also, I'm not sure what the ??? is for.

normal_inv_cdf(p | mu, sigma);

would be fine.
 

[Andrew Gelman]

My main point is that if someone just writes normal(a|b,c) or binomial(a|b,c), etc., it should return the log density.  It’s fine with me to also allow explicit _pdf or _log_pdf or _cdf or _ccdf or _log_ccdf.

A

[Bob Carpenter]

The issue is whether we even have non-log versions of pdf, cdf, and ccdf.

I'd strongly prefer not to have multiple ways to write basic functions.
If both normal(a | b, c) and normal_log(a | b, c) are in play, then users
will see both and assume they are different functions, and naturally conclude
the former is not on the log scale. As Aki has pointed out, if we just
write normal(a | b, c), then everyone's going to assume it's *not* on the
log scale from looking at it. He's right, and this is a big big step.

Andrew --- please not the spacing. Compared to commas, we're losing a space.
And I don't like normal(a| b, c) or normal(a|b, c) or normal(a|b,c).

As to inverse CDFs, you're right that it doesn't make sense to have log
output. But what about logit input? Or do we have complementary inverse
cdfs to deal with the overflow to 1 problem on the 1 side of (0,1)?

- Bob

[Andrew Gelman]

I don’t think we should have a function called normal_log.

We can have normal_pdf, normal_log_pdf, normal_cdf, normal_log_cdf. Just plain normal is shorthand for normal_log_pdf.

Also, I don’t understand what you mean by “Compared to commas, we’re losing a space.” I thought these spaces had no effect on the Stan code?

[Ben Goodrich]

On Wednesday, January 6, 2016 at 5:09:29 PM UTC-5, Bob Carpenter wrote:

As to inverse CDFs, you're right that it doesn't make sense to have log
output.  But what about logit input?  Or do we have complementary inverse
cdfs to deal with the overflow to 1 problem on the 1 side of (0,1)?  

R has an option for log input, but I don't think it is that big a deal for Stan. If the probability is nearly 0 or 1, then the inverse CDF will probably not be very stable even if we give it more stable input.

Ben

[Bob Carpenter]

Won't that depend on how we write the code? I don't know anything
about this case, but I don't see any obstacle in principle to using
probabilities near 0 or 1 in inverse CDFs as long as the output doesn't
overflow or underflow.

But if we have something like gamma_inv_cdf(epsilon, a, b) and epsilon is
close to 0, then the output's going to be close to zero. Or if we write
beta_inv_cdf(epsilon, a, b) and epsilon is close to 1, then the output is
close to 1 and it cries out for a a complementary inverse CDF on the log
scale for input and output.

Anyway, we can put off CDFs for the time being because we don't have any
implemented or any path to implement them until someone has some time
and energy for it.

- Bob

[Bob Carpenter]

Sorry to confuse the issue --- I meant "_lpdf" or "_pdf", not "_log".

I don't want users to be able to write a "shorthand" and "longhand"
version of the same function. Let's just pick one and stick with it.

Stan's parser won't care about whitespaces and will let you write
any foolish syntactic convention you want (other than foregoing semicolons
[a la R] or dropping parens in favor of indenting [a la Python]).

But if you count the spaces in "normal(a, b, c)" and "normal(a | b, c)"
you'll see that the latter is one character wider. Size matters. As
does having enough space between things to read.

- Bob

[Ben Goodrich]

On Wednesday, January 6, 2016 at 5:25:46 PM UTC-5, Bob Carpenter wrote:

Won't that depend on how we write the code?  I don't know anything
about this case, but I don't see any obstacle in principle to using
probabilities near 0 or 1 in inverse CDFs as long as the output doesn't
overflow or underflow.  

I was just speculating that the stability in the tails of the output is probably worse than the stability in the input.
 

Anyway, we can put off CDFs for the time being because we don't have any
implemented or any path to implement them until someone has some time
and energy for it.

I agree we can decide later, although the existing ones (inv_Phi, logit) only input on the probability scale.

Ben

[Michael Betancourt]

+1

[Andrew Gelman]

Bob:
Now I’m confused. See below.
A

> On Jan 6, 2016, at 5:27 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> Sorry to confuse the issue --- I meant "_lpdf" or "_pdf", not "_log".
>
> I don't want users to be able to write a "shorthand" and "longhand"
> version of the same function. Let’s just pick one and stick with it.

On Jan 6, 2016, at 4:09 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Then I believe Andrew wanted to allow "normal" as an abbreviation
for “normal_pdf".

—

What I’d like is for “normal” to be the abbreviation for “normal_log_pdf”.

[Daniel Lee]

On Wed, Jan 6, 2016 at 6:30 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

What I’d like is for “normal” to be the abbreviation for “normal_log_pdf”.

What Bob, Michael, and I (and perhaps others) would like is there to be no abbreviation. Either we can have users write "normal_log_pdf" every time or "normal" -- we don't want to have both options kicking around.
 

[Andrew Gelman]

I’m confused because earlier Bob wrote this:  “Then I believe Andrew wanted to allow “normal" as an abbreviation for “normal_pdf”.”  And I didn’t think he objected.  All I did was clarify that normal should be an abbreviation for normal_log_pdf, not normal_pdf.

A

[Bob Carpenter]

I don't want "normal" to be an abbreviation for anything.

Either the function to compute the log pdf is "normal(y | mu, sigma)"
or it's "normal_lpdf(y | mu, sigma)" or it's something else.

Allowing both is going to be very confusing to users. They'll
immediately come to the conclusion there's a difference and try to
figure out which one they should use.

- Bob

[Bob Carpenter]

I'm objecting now that I've thought about it.

Maybe I shouldn't try to summarize other people's points
of view in emails.

- Bob

[Daniel Lee]

> On Jan 6, 2016, at 11:23 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> I'm objecting now that I've thought about it.

What are you objecting to? I'm lost.

>
> Maybe I shouldn't try to summarize other people's points
> of view in emails.

It's been super helpful for me.

[Bob Carpenter]

> On Jan 6, 2016, at 11:25 PM, Daniel Lee <bea...@alum.mit.edu> wrote:
>
>
>> On Jan 6, 2016, at 11:23 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>>
>> I'm objecting now that I've thought about it.
>
> What are you objecting to? I'm lost.

Having two names for the same function. We already do that
with max/fmax and with abs/fabs.

- Bob

[Andrew Gelman]

OK, we can discuss again. I see where you’re coming from. I guess it’s fine if normal(a | b, c) is the log density, then we’d also have normal_log_cdf(a | b, c) and normal_cdf(a | b, c). That works. And if anyone wants to compute the normal pdf (for example, if they’re defining a mixture model), they’d just have to bite the bullet and do exp(normal(a | b, c)).
Is that ok?

[Andrew Gelman]

Oh yeah, I hate that fabs thing!

[Ido Rosen]

On Wed, Jan 6, 2016 at 11:30 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

OK, we can discuss again.  I see where you’re coming from.  I guess it’s fine if normal(a | b, c) is the log density, then we’d also have normal_log_cdf(a | b, c) and normal_cdf(a | b, c).  That works.  And if anyone wants to compute the normal pdf (for example, if they’re defining a mixture model), they’d just have to bite the bullet and do exp(normal(a | b, c)).
Is that ok?

Why not just use function composition / cdf + log as higher-order functions?  It's almost the same amount of typing: normal_log_cdf -> log(cdf(normal(a|b,c)))) or logcdf(...) that gets resolved lazily(/symbolically?) to whatever normal_log_cdf(a|b,c) is.  It could be more flexible and require less aliases, maybe?

[Daniel Lee]

On Wed, Jan 6, 2016 at 11:30 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

OK, we can discuss again.  I see where you’re coming from.  I guess it’s fine if normal(a | b, c) is the log density, then we’d also have normal_log_cdf(a | b, c) and normal_cdf(a | b, c).  That works.  And if anyone wants to compute the normal pdf (for example, if they’re defining a mixture model), they’d just have to bite the bullet and do exp(normal(a | b, c)).
Is that ok?

Sorry, I lost a bit of what you were suggesting. I think:

normal(a | b, c) should be in R: dnorm(a, b, c, log = TRUE)

normal_cdf(a | b, c) should be in R: pnorm(a, b, c, log=TRUE)

These things will result in a compiler error because they won't exist:

normal_log(a | b, c)

normal_lpdf(a | b, c)

normal_log_cdf(a | b, c)

Are we on the same page?

And yes, if they really, really want things on the normal scale (what's this called? exponentiated scale? probability scale? sorry... it's late and I can't think of the appropriate terminology), they will have to exponentiate.

So far, we've been carrying around two notations and it's just made it harder to think about. Of course, we're not doing this with other functions, just with the probability distributions. It's not like the sin() operator will return the log of the sin(), but I think it's ok since we're treating conditional distribution calls as first class citizens in the language.

Right now, my opinion is that all our distribution functions should be on the log scale. This includes:

- probability distribution / mass function

- cumulative distribution function

- inverse cumulative distribution function

[Andrew Gelman]

Hey, that’s a cute idea but I don’t think it makes sense here.  Even setting aside the hassle for users of having a new kind of structure in the language, there’s the problem that I’d like normal(a | b, c) to be the log pdf already.  I think normal_cdf and normal_log_cdf are fine, but I think an expression like cdf(normal(…)) would just invite confusion!

[Daniel Lee]

Ok. Cool. I also object to having multiple names for the same thing.

[Andrew Gelman]

Actually I was proposing something slightly different:

normal(a | b, c) should be in R: dnorm(a, b, c, log = TRUE)

normal_cdf(a | b, c) should be in R: pnorm(a, b, c, log=FALSE)
normal_log_cdf(a | b, c) should be in R: pnorm(a, b, c, log=TRUE)

But your plan could work too.  I’d like to avoid the _pdf or _lpdf suffix in any case because that would clutter the model.  To me, the #1 thing is getting the commonly-used distributions in a clean format, then we worry about things like cdfs being logical.

As I was saying at meeting on Tues, the main Bayes formula is p(theta | y) propto p(theta) p(y | theta).

But now I’d like to write it as lp(theta | y) = const + lp(theta) + lp(y | theta).

In which case it’s clear that normal, poisson, binomial, etc., are just the lp’s.

As with our sd/variance notation, we’re making Stan inconsistent with BDA, and that’s too bad, but I think it’s worth it because of the transparency of the future Stan models.

[Ben Goodrich]

On Wednesday, January 6, 2016 at 11:51:19 PM UTC-5, Daniel Lee wrote:

Right now, my opinion is that all our distribution functions should be on the log scale. This includes:

- probability distribution / mass function

- cumulative distribution function

- inverse cumulative distribution function

Not for inverse CDF, which can be a negative number, and even if the distribution is defined as positive, its log is not what any user would expect.

Ben

[Daniel Lee]

Yes, you're right. I was thinking (1 - cdf), not inverse cdf. My fault.

 

Ben

[Daniel Lee]

On Wed, Jan 6, 2016 at 11:55 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Actually I was proposing something slightly different:

normal(a | b, c) should be in R: dnorm(a, b, c, log = TRUE)
normal_cdf(a | b, c) should be in R: pnorm(a, b, c, log=FALSE)
normal_log_cdf(a | b, c) should be in R: pnorm(a, b, c, log=TRUE)

But your plan could work too.  I’d like to avoid the _pdf or _lpdf suffix in any case because that would clutter the model.  To me, the #1 thing is getting the commonly-used distributions in a clean format, then we worry about things like cdfs being logical.

The decision regarding the cdfs must be made at the same time. We don't want to change the language twice.

For cdfs being logical, I was shooting for usage, example below.

 

As I was saying at meeting on Tues, the main Bayes formula is p(theta | y) propto p(theta) p(y | theta).

But now I’d like to write it as lp(theta | y) = const + lp(theta) + lp(y | theta).

I'm thinking about how to write a truncated model. I was really thinking about something like:

lp(mu, sigma | y, upper) = const + normal(y | mu, sigma) - normal_cdf(upper | mu, sigma)

I think we should have expressions for:

cdf

1 - cdf

cdf - cdf

So... something like:

  lp << normal(y | mu, sigma) - normal_cdf(upper |  mu, sigma);
  lp << normal(y | mu, sigma) - normal_1mcdf(lower |  mu, sigma);
  lp << normal(y | mu, sigma) - normal_cdf(lower, upper |  mu, sigma);

Or something similar? (please forgive me if I switched the truncation direction)

I think we also want to keep the T[] notation, but I could be convinced otherwise:

  normal(y | mu, sigma) T[, upper]

  normal(y | mu, sigma) T[lower, ]
  normal(y | mu, sigma) T[lower, upper]

[Andrew Gelman]

That works for me.  All on log scale makes these models clean.

[Bob Carpenter]

> On Jan 6, 2016, at 11:30 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> OK, we can discuss again. I see where you’re coming from. I guess it’s fine if normal(a | b, c) is the log density, then we’d also have normal_log_cdf(a | b, c) and normal_cdf(a | b, c). That works. And if anyone wants to compute the normal pdf (for example, if they’re defining a mixture model), they’d just have to bite the bullet and do exp(normal(a | b, c)).
> Is that ok?
>

We don't want to encourage anyone to ever use exp(normal(a | b, c)).
If they want a mixture model, they should use log_sum_exp as we
outline in the manual.

> As I was saying at meeting on Tues, the main Bayes formula is p(theta | y) propto p(theta) p(y | theta).
>
> But now I’d like to write it as lp(theta | y) = const + lp(theta) + lp(y | theta).

We really need a "propto" on the log scale. I've been trying to figure
out how to say it for years without writing "const +", which is awkward
because it uses a five-letter word where every other variable is a Greek or
roman letter.

> So... something like:
> lp << normal(y | mu, sigma) - normal_cdf(upper | mu, sigma);
> lp << normal(y | mu, sigma) - normal_1mcdf(lower | mu, sigma);
> lp << normal(y | mu, sigma) - normal_cdf(lower, upper | mu, sigma);

But we'll use "ccdf" instead of "1mcdf". But now the relation
is on the log scale, so with no logs in sight, that's

exp(foo_ccdf(y | theta)) = exp(1 - exp(foo_cdf(y | theta)))

Normally it's just

CCDF(...) = 1 - CDF(...)

> I think we also want to keep the T[] notation, but I could be convinced otherwise:
> normal(y | mu, sigma) T[, upper]
> normal(y | mu, sigma) T[lower, ]
> normal(y | mu, sigma) T[lower, upper]

Absolutely.

- Bob

[Bob Carpenter]

It's on the to-cleanup list. But right now, it's
deprecated opposite the way you like it. But I think
we've mastered enough C++ that we can take on the
namespace issues with plain-old max and min and abs.

- Bob

[Michael Betancourt]

I am strongly with Bob and Dan.  Shortcuts for "commonly-used distributions" is another way of saying "inconsistent behavior based on contexts that can change in time".  Especially when those contexts are based on Andrew's rapidly evolving perspectives!  We want the language to be explicit as possible so it will be robust to our changing understanding and best practices.