New Meeting Format

[Michael Betancourt]

With the increased number of participants at each dev meeting
we need to get a little more organized and unless there are
any objections I’m going to change things up on Tuesday.

In order to get important business out of the way before the
meeting becomes too tangential the meeting will now follow a
rigid progression:

1. Roundtable

Everyone very briefly reviews that they’ve been up to in the past
week (less than a minute, maybe more if they have been particularly productive).
This could be removed if there is enough support, but I think it’s
nice to touch base with everyone.

2. Unresolved Topics

These are topics carried over from the past meeting — they won’t
be applicable this week unless there are important issues that
were raised last week. Let me know.

3. New Topics

These are topics of discussion that need to be addressed immediately
or before the next meeting, including questions about modeling or
development, well-posed requests for Stan features, or decisions that
have to be made. New topics not immediately addressed at the meeting
will become “Unresolved topics” to be reviewed the following meeting.

4. Open Discussion Topics

This is what our meetings typically devolve into right now, only hopefully
we can make it a bit more organized. After all proposed discussion
topics have been covered we can open the meeting to other topics.

WHAT THIS MEANS FOR YOU:

If you hate this then let me know and we can try to find a compromise.

If there is anything that you want to have addressed next meeting
add it under “New Topics” at https://github.com/stan-dev/stan/wiki/Stan-Development-Meeting-Agenda.
If there are any any topics that you would like to discuss add then under
“Open Discussion Topics” at the same place.

This page is not meant for recording minutes or providing an archive
of meeting discussions — it’s just for organizing the meetings week
to week.

[Bob Carpenter]

Thanks. I'm in favor of more structure and this
looks good.

- Bob

> --
> 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.
>
>

[Aki Vehtari]

Great! I had been thinking to ask something like this!

Aki

[Krzysztof Sakrejda]

Sounds good to me, I appreciate focused meetings. Krzysztof

[Andrew Gelman]

All sounds good to me.
A

> On Mar 5, 2016, at 1:34 PM, Michael Betancourt <betan...@gmail.com> wrote:
>

[Michael Betancourt]

Reminder that if you want to discuss something at the meeting tomorrow make sure to add an item
at http://github.com/stan-dev/stan/wiki/Stan-Development-Meeting-Agenda _before_ the meeting.

[Aki Vehtari]

Could there be brief summary of some of the topics discussed for those who were not able to join the meeting or had to leave early? As the topics have been cleared from the agenda it seems that Stan 3 and Stan book questions received satisfactory answers, which I would be interested to know, too :)

Aki

[Michael Betancourt]

I’m keeping the meeting on track, not keeping minutes. If someone
wants to do that go ahead, but they should also figure out a way to
archive meeting minutes.

[Bob Carpenter]

Sure, if we can get a volunteer to take notes.

The only major decisions I'm aware of were:

1. use = interchangeably with <- and deprecate <- (and already
merged on develop since the meeting)

2. Daniel's going to take over the command infrastructure.

There was inconclusive discussion on when and what Stan 3 are
going to entail, and also about putting someone in
charge of making decisions at some level (package, repo, etc.).

I think the book's moot until we decide what's going into
Stan 3. I don't want to have to write it twice.

I asked Andrew how GMO was different from MML and didn't
understand the response (we've been going in circles in email
on the same topic --- statisticians and computer scientists
don't speak the same pseudocode).

I asked Andrew why he wanted GMO when we had ADVI and
I think it came down to his not understanding and/or trusting
ADVI (Andrew's comments are always subtly hedged, so I'm
sure I'm doing him and ADVI a disservice trying to summarize).

- Bob

[Michael Betancourt]

I asked Andrew how GMO was different from MML and didn't
understand the response (we've been going in circles in email
on the same topic --- statisticians and computer scientists
don't speak the same pseudocode).  

I asked Andrew why he wanted GMO when we had ADVI and
I think it came down to his not understanding and/or trusting
ADVI (Andrew's comments are always subtly hedged, so I'm
sure I'm doing him and ADVI a disservice trying to summarize).

1) The two algorithms answer two different questions.

ADVI tries to estimate _all_ marginal _means_.

MML tries to estimate _any_ marginal expectation with

respect to _a few chosen parameters_ (typically hyperparameters

in a hierarchical model).

2) Consequently the two algorithms lives on different points

of the speed/accuracy frontier.  Which is better for a given

application depends on the details of that application, which

means we’re entering into “expert-level knowledge needed”

territory.

[Bob Carpenter]

I thought MML wasn't about expectations, but was a max-max
algorithm, as in first doing the marginal max of high-level
params

phi* = argmax_phi p(phi | y),

then computing the max of the low-level params

alpha* = argmax_alpha p(alpha, phi* | y)

Where's the expectation? Anyway, no rush on answers if you've
got better things to do---I can eventually sort this out with
enough back-and-forth with Andrew and Dustin. The problem I have
with the way statisticians write algorithms is that it's never
clear to me what they're trying to compute.

And then I thought once they got (phi*, alpha*) that there
was a subsequent step of importance-adjusted draws from the
multivariate approximation at (phi*, alpha*) with covariance
equal to inverse Hessian, so that the mean returned by GMO
was not necessarily going to be (phi*, alpha*) after the importance
weighting.

So I'm still baffled by the algorithm.

- Bob

[Dustin Tran]

Michael has it right. MML only concerns the marginal max of “high-level params” in that equation you write.

To do this, we do gradient descent. Gradient descent requires gradients of the marginal posterior density. This is hard because it’s a gradient of an expectation. Thus there’s an “E-step” in dealing with that expectation.

I think your confusion is about the relationship between EM and MM. You can view EM-type procedures as doing an “MM”, where computing expectations corresponds to maximizing something (for example, variational EM is an approximate “MM”).

Dustin

[Michael Betancourt]

Remember that everything in statistics (really everything in measure theory)
is expectations. Modes are only useful as a means to approximate
expectations, for example by approximating a target distribution with
a delta function at the mode,

\int dx p(x) f(x) \approx \int dx \delta_{x*} (x) f(x) = f(x*)

where

x* = argmax p(x).

Of course this isn’t parameterization invariant and is an atrocious estimator
in high dimensions unless there is an _insane_ amount of symmetry which
statistically translates to little or no corrections.

MML methods try to take a high dimensional problem and reduce it to a
lower dimensional one. In particular, they partition the parameters x
into two sets, a low dimensional set y and a high dimensional set z,
x = (y, z). If we could marginalize out the y analytically then we could
find the mode of p(y) which might be a useful approximation to some
expectations. In practice we can’t do that, but we can approximate it,
which is what MML methods do.

The key computational step in MML methods is estimating the
conditional p(z | y) and then estimating

p(y) = \int dz p(z | y) p(y).

The conditional and the expectation are typically done with Gaussian
approximations. The importance sampling step in GMO just tries to
correct these approximations to get more accurate estimates of the
marginal density.

[Andrew Gelman]

Hi I don’t know if this will help, but . . .

1.  MML vs GMO

MML is a goal.  GMO is an algorithm.

MML is marginal maximum likelihood (more generally, marginal maximum penalized likelihood).  GMO is gradient-based marginal optimization, an algorithm (or class of algorithms) for reaching the MML solution using a particular approach to computing a simulation-based estimate of the gradient of the marginal likelihood function.  The GMO algorithm is similar to the ADVI algorithm in that it approximates the gradient directly, in contrast to the usual MML algorithms with approximate the marginal lilkelihood and then go from there.  It turns out it can make sense to go for the gradient directly, as we can then make use of the path sampling identity.

MML is to GMO as VB is to ADVI.

2.  Averaging and optimization in GMO

Conditional on the hyperparameters, GMO computes a simulation-based estimate of the gradient of the marginal likelihood.  This is the gradient of an integral, thus an integral of a gradient, thus it’s an integral.  But the way we simulation-estimate this integral is to start with a conditional mode-based approximation.  That’s why I say that GMO has an external optiimzation routine and an internal optimization routine. Really, it has an internal numerical integration routine, but it operates by first performing a conditional optimization.  It would be possible to run GMO without doing an internal optimzation; for example, the internal integral could be estimated using Nuts.  I’ve just been assuming that the mode-based approximation will be faster.

So, the internal integration is a necessary part of GMO.  The internal optimzation is not necessary, but it’s how we’ve been doing things.  The internal optimization is only there as a way to set up the internal integration.

[Bob Carpenter]

> On Mar 8, 2016, at 9:49 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> Hi I don’t know if this will help, but . . .
>
>
> 1. MML vs GMO
>
> MML is a goal. GMO is an algorithm.
>
> MML is marginal maximum likelihood (more generally, marginal maximum penalized likelihood). GMO is gradient-based marginal optimization, an algorithm (or class of algorithms) for reaching the MML solution using a particular approach to computing a simulation-based estimate of the gradient of the marginal likelihood function. The GMO algorithm is similar to the ADVI algorithm in that it approximates the gradient directly, in contrast to the usual MML algorithms with approximate the marginal lilkelihood and then go from there. It turns out it can make sense to go for the gradient directly, as we can then make use of the path sampling identity.
>
> MML is to GMO as VB is to ADVI.

If I understood all the correctiosn to my language earlier,
MML is just this bit

phi* = argmax_phi p(phi | y)

and the following bit is something I shouldn't call "MML":

alpha* = argmax_alpha p(alpha, phi* | y)

but which has no name.

> 2. Averaging and optimization in GMO
>
> Conditional on the hyperparameters, GMO computes a simulation-based estimate of the gradient of the marginal likelihood. This is the gradient of an integral, thus an integral of a gradient, thus it’s an integral. But the way we simulation-estimate this integral is to start with a conditional mode-based approximation. That’s why I say that GMO has an external optiimzation routine and an internal optimization routine. Really, it has an internal numerical integration routine, but it operates by first performing a conditional optimization. It would be possible to run GMO without doing an internal optimzation; for example, the internal integral could be estimated using Nuts. I’ve just been assuming that the mode-based approximation will be faster.
>
> So, the internal integration is a necessary part of GMO. The internal optimzation is not necessary, but it’s how we’ve been doing things. The internal optimization is only there as a way to set up the internal integration.

I also got that part of the algorithm, at least at a high
level, by analogy to ADVI, which Alp wrote in language I
could follow.

But I thought GMO went even further in taking the MML + alpha* thing
point estimate (phi*, alpha*) and created a multivariate normal
approximation of the unconstrained posterior with mean (phi*, alpha*)
and covariance given by the inverse Hessian at (phi*, alpha*).

Then you take a bunch of draws from that multivariate normal and
importance weight them, then transform back to the constrained space.
Now if I take the importance-weighted average of the constrained draws,
I'm not guaranteed to get (phi*, alpha*) as the posterior mean.

So that'd be a rather different algorithm than MML. But
ADVI really is just an algorithm for computing a particular
variational approximation and otherwise, the solution is the
usual variational approximation.

- Bob

>
>
>> On Mar 8, 2016, at 5:25 PM, Michael Betancourt <betan...@gmail.com> wrote:
>>
>>
>>> I asked Andrew how GMO was different from MML and didn't
>>> understand the response (we've been going in circles in email
>>> on the same topic --- statisticians and computer scientists
>>> don't speak the same pseudocode).
>>>
>>> I asked Andrew why he wanted GMO when we had ADVI and
>>> I think it came down to his not understanding and/or trusting
>>> ADVI (Andrew's comments are always subtly hedged, so I'm
>>> sure I'm doing him and ADVI a disservice trying to summarize).
>>
>> 1) The two algorithms answer two different questions.
>>
>> ADVI tries to estimate _all_ marginal _means_.
>>
>> MML tries to estimate _any_ marginal expectation with
>> respect to _a few chosen parameters_ (typically hyperparameters
>> in a hierarchical model).
>>
>> 2) Consequently the two algorithms lives on different points
>> of the speed/accuracy frontier. Which is better for a given
>> application depends on the details of that application, which
>> means we’re entering into “expert-level knowledge needed”
>> territory.
>
>

[Dustin Tran]

> But I thought GMO went even further in taking the MML + alpha* thing
> point estimate (phi*, alpha*) and created a multivariate normal
> approximation of the unconstrained posterior with mean (phi*, alpha*)
> and covariance given by the inverse Hessian at (phi*, alpha*).

GMO maximizes the marginal likelihood p(y | phi) to get phi^*. As you said, this is the problem of MML; everything else is an extension off that point estimate.

re:covariance. It’s an extension. As Bayesians, we view phi^* as a point mass approximation to the marginal posterior density p(phi | y). A cheap way to get uncertainty is to estimate the covariance of p(phi | y) using the inverse Hessian evaluated at phi^*. We can use this with phi^* to assume a normal approximation to p(phi | y), with mean phi^* and covariance H^{-1}(phi^*). This is like how in lm, it outputs both regression coefficients and their standard errors. Note: this still marginalizes over alpha.

Dustin

[Andrew Gelman]

> On Mar 8, 2016, at 10:07 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>

>> On Mar 8, 2016, at 9:49 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>>
>> 2. Averaging and optimization in GMO
>>
>> Conditional on the hyperparameters, GMO computes a simulation-based estimate of the gradient of the marginal likelihood. This is the gradient of an integral, thus an integral of a gradient, thus it’s an integral. But the way we simulation-estimate this integral is to start with a conditional mode-based approximation. That’s why I say that GMO has an external optiimzation routine and an internal optimization routine. Really, it has an internal numerical integration routine, but it operates by first performing a conditional optimization. It would be possible to run GMO without doing an internal optimzation; for example, the internal integral could be estimated using Nuts. I’ve just been assuming that the mode-based approximation will be faster.
>>
>> So, the internal integration is a necessary part of GMO. The internal optimzation is not necessary, but it’s how we’ve been doing things. The internal optimization is only there as a way to set up the internal integration.
>
> I also got that part of the algorithm, at least at a high
> level, by analogy to ADVI, which Alp wrote in language I
> could follow.
>
> But I thought GMO went even further in taking the MML + alpha* thing
> point estimate (phi*, alpha*) and created a multivariate normal
> approximation of the unconstrained posterior with mean (phi*, alpha*)
> and covariance given by the inverse Hessian at (phi*, alpha*).
>
> Then you take a bunch of draws from that multivariate normal and
> importance weight them, then transform back to the constrained space.
> Now if I take the importance-weighted average of the constrained draws,
> I'm not guaranteed to get (phi*, alpha*) as the posterior mean.
>
> So that'd be a rather different algorithm than MML. But
> ADVI really is just an algorithm for computing a particular
> variational approximation and otherwise, the solution is the
> usual variational approximation.

At each “external” step (that is, guess of phi), GMO approximates the gradient of the integral p(phi,alpha|y) d alpha. My particular idea of a GMO algorithm (so, I guess we could call this “GMO,” or at least the first recorded flavor of GMO) indeed first computed alpha*, then the normal approximation, then simulations, importance weighting, etc. Don’t worry about constrained space for now; I’m assuming we’re entirely working in the unconstrained space. The algoithm as written, with importnace sampling, gives a simulation-consistent estimate of the gradient of the integral. So it really is a simulation-based approach to getting MML. In practice, of course we don’t have infinitely many simulation draws so it’s only an approximation to MML, but that’s the same as ADVI and any other simulation-based integration.

[Dustin Tran]

There’s a technical snippet that’s worth noting (to Andrew):

If we apply the path sampling identity to take the gradient of the marginal density (not the log marginal density), then there is no longer a ratio of expectations but one expectation. Importance sampling the one expectation gives us a simulation-consistent and also unbiased estimate of that gradient; therefore we can say it performs exact MML. Importance sampling a ratio of expectations as we currently do leads to a biased but simulation-consistent estimate (Appendix E in the write-up has the math). In practice, there’s no reason to do this over the ratio of expectations because optimizing on the original scale is more numerically unstable than optimizing on the log scale. But it’s nice to know that this variant exists.

Dustin

[Bob Carpenter]

OK, I see now. I didn't realize the importance
sampling was for the integral --- I thought it
was for approximate draws from the posterior.
And now I get GMO:MML::ADVI:VI --- we're just
talking about p(phi | y) here.

What you do with the phi* you get out of GMO
is your business!

Which brings up the question of what you (or lme4)
do for alpha* if MML only gives you a
Laplace approximation to p(phi | y).

And is what people call "empirical Bayes" just
full Bayes for p(alpha | phi*, y)?

- Bob

P.S. Dustin --- I'm putting this in Bayesian terms
and avoiding max likelihood and std errors on purpose,
because we're going to have a prior and use the usual
Bayesian manipulations to get the MML estimate as the
approximate posterior mean if the prior is uniform.

[Andrew Gelman]

We’ll get our alpha* based on our phi*, or else based on our 1000 importance-weighted draws of phi. If this phi* or these importance-weighted draws are just approximations (because calculated based on an approximation of the internal integral) then, sure, alpha* won’t be quite right.

As to what people call “empirical Bayes,” that’s a flexible term that means different things to different people. But, yes, it often corresponds to plugging in a point estimate of the hyperparameter.

I’ve been talking about GMO rather than MML because we have typically talked about algorithms.

Thus, GMO, ADVI, NUTS, RMC

Or MML, VB, full Bayes, full Bayes.

[Dustin Tran]

> And is what people call "empirical Bayes" just
> full Bayes for p(alpha | phi*, y)?

Following Andrew’s point, I don’t know a good name for this. Another phrase thrown around is "Bayesian inference with data-dependent priors”. It is shorthand for "Bayesian inference with priors whose hyperparameters are point-estimated from the data”, which is what Bob wrote above in the case phi are prior hyperparameters. The general notion of “data-dependent priors” is also technically true for full Bayes, since the full posterior p(alpha, phi | y) has “data-dependent priors”. So I don’t know what to call it (if we should even call it something).

> P.S.  Dustin --- I'm putting this in Bayesian terms

> and avoiding max likelihood and std errors on purpose

Cool, we’re on the same page.

Dustin

[Allen B. Riddell]

Creating a wiki page for each set of minutes would be one easy way.

I don't mind taking minutes on the days which I can attend (e.g., days
when the agenda includes PyStan-related stuff)

[Bob Carpenter]

I still don't quite understand what GMO computes (specifically
how you get alpha*).

At least with HMC and NUTS I knew the goal was to draw a sample
from the posterior!

It's traditional in computer science to define the function
being computed before introducing an algorithm to compute it.
Working backwards to the function from the algorithm can be painful,
especially for those of us who don't know what MML or lme4 was doing.

- Bob

[Andrew Gelman]

With HMC and NUTS I knew the goal was to draw a sample from the posterior.

With GMO the goal is to compute the phi that maximizes p(phi | y) = int p(phi, alpha | y) d alpha.

There are also variants of GMO that compute some approximation. But basic GMO aims for that argmax.

Knowing about lme4 is important if we’re interested in efficiency improvements for GMO. But you don’t need to know about lme4 to understand GMO.

A