Sum of lognormals?

[Geoff Davis]

Hi Stan folks -

I've got a measurement that I'm trying to model as a convex combination of k log-normals.  The weights in the convex combination are known and vary from measure to measure; what I'm trying to do is to estimate the parameters of the log-normals.  

I've used Stan happily for HLMs where everything is Gaussian, but I'm a bit puzzled about how to express this new kind of relationship since the log-normals don't conveniently collapse.

For concreteness, suppose

y_i = a_{i,1} x_{i,1} + a_{i, 2} x_{i, 2} + a_{i, 3} x_{i, 3}

where

x_{i, j} ~ lognormal(mu_j, sigma_j)

and I have priors on mu_j and sigma_j.

The values of y_i are all observed.

The values a_{i, j} are all known, and for all i, a_{i, 1} + a_{i, 2} + a_{i, 3} = 1.

A few possibilities:

1) Express x[i, 1] in terms of x[i, 2] and x[i, 3], e.g.

real diff;

for (i in 1:nsample) {

  x[i, 2] ~ lognormal(mu_2, sigma_2);

  x[i, 3] ~ lognormal(mu_3, sigma_3);

  diff <- (y[i] - a[i, 2] x[i, 2] - a[i, 3] x[i, 3]) / a[i, 1];

  diff ~ lognormal(mu_1, sigma_1);

}

2) Introduce a small source of additional noise with known variance.  Then I can declare something like

for (i in 1:nsample) {

  x[i, 1] ~ lognormal(mu_1, sigma_1);

  x[i, 2] ~ lognormal(mu_2, sigma_2);

  x[i, 3] ~ lognormal(mu_3, sigma_3);

  y[i] ~ lognormal(log(a[i, 1] * x[i, 1] + a[i, 2] * x[i, 2] + a[i, 3] * x[i, 3]), epsilon);

}

3) Something else?  I'm wondering if there's some sort of trick I need to do along the lines of the way Stan deals with simplexes.

Thanks!

Geoff

[Tamas Papp]

For (1), I think you would need a Jacobian correction. It may be the
most straightforward solution to code, but may not sample efficiently,
since diff will be < 0 for some parameters.

Whether (2) is appealing depends on what you are modeling. I would think
about making epsilon a parameter. Not sure about identification though,
do we know anything about x_{i,j} -- eg do they come from the same
distribution for a given j?

For (3) you might make

z_{i,j} = a_{i,j} x_{i,j} / y_i

and then z_i is in a 3-simplex. Then parametrize either in a_{i,j} or
x_{i,j}.

Unless there is a strong reason not to (again, depends on your model, of
which I don't have domain-specific knowledge), I would go for (2). I
have spent about a week recently coding up clever transformations for a
model with constraints similar to (1)/(3) (in spirit; they were actually
nonlinear), only to find that (2) makes more sense and is easier to
extend, and I could assign meaning to the error.

Best,

Tamas

[Bob Carpenter]

I would not recommend 2, because I don't think you can easily fake
a hard constraint with a soft constraint.

The problem for (1) isn't so much that you need a Jacobian with this transform

> x[i,1] <- (y[i] - a[i, 2] x[i, 2] - a[i, 3] x[i, 3]) / a[i, 1];

but rather that it's that you're trying to provide three distributions for
x[i,1:3] on what are effectively two variables, x[i,2:3] and a funciton of them.
So you need to work out what the implications are. It's like you're doing
something like this:

y[1] ~ foo(...);
y[2] ~ foo(...);
y[1] + y[2] ~ foo(...);

It's much easier to work with a model conceptually if you can imagine
the generative process leading to the data.

You can see an instance of that with the way we handle simplexes under the hood.
If you declare

parameters {
simplex[5] theta;
...

then there are really 4 underlying variables. Under the hood Stan does
the transform and takes care of the Jacobian. Then if you do this

model {
theta ~ dirichlet(alpha);

everything's OK because the Dirichlet has support over the simplexes.
But you're constraining things more and using mismatched distributions that aren't
distributions over your constrained variables. It'd be like for the simplex doing this

theta ~ normal(0, 1);

and then I'm not sure what it'd mean. Someone else will hopefully know and
set us both straight here.

- Bob

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

[Geoff Davis]

Thanks to both of you for the feedback.  Yes, I'm having a hard time thinking about how to handle the constraint.  

Here's more about the generative process (good suggestion):

* I've got a set of blood samples.  Each sample is a mixture of K cell types (monocytes, lymphocytes, etc), and I've measured the fraction of each cell type in each sample.  I'm assuming for now that these measurements are noiseless.  The quantity a_{i, j} is the fraction of cell type j in sample i.

* I measure a substance of interest in each sample (my y_i terms).  I know that each cell type produces a different amount of the substance (x_{i, j} is the amount of substance produced by cell type j in sample i), and y_i = \sum_j a_{i, j} x_{i, j}.  There is variation in the amount of substance produced per cell type from sample to sample, and I know a priori that the distribution is log-normal.  Think of it as biological variation from person to person, and assume that the samples all come from different people.  More concretely, if I were able to get samples from a set of people that were purified so they contained only a single cell type, I know a priori that the amount of substance I'd measure in each sample would be distributed roughly log-normally.  (Each cell type j has its own log-normal parameters)

* There's measurement error in the mix as well, and I know a priori that the errors are multiplicative.

The model in which I've introduced the extra error term I think of as being a model in which I've got both biological and technical variation represented explicitly.  If I let the measurement error term be a parameter, I get identifiability problems, since increasing the technical variation is roughly equivalent to increasing all the sources of biological variation (i.e. the variances in the log-normals for individual cell types).

The constrained model is an attempt to assume the measurement error can be subsumed into the biological variation terms.

Adding in the extra multiplicative error term is an attempt to model the measurement error explicitly, and my making it a known quantity is a way to get around identifiability issues.

Does that make more sense?

Geoff

[Geoff Davis]

Bob, your comment spurred an interesting thought:

One way to think about the constrained joint distribution of the x[i, j] for sample i is that it's the unconstrained distribution restricted to the hyperplane on which the constraint is satisfied, i.e. where x[i,1] = (y[i] - a[i, 2] x[i, 2] - a[i, 3] x[i, 3]) / a[i, 1].  The headache is that what I want is a conditional probability where the condition is the hyperplane constraint, i.e. I need to normalize the constrained joint distribution, and the normalization factor depends on the parameters of interest.  If I'm understanding correctly, I need to compute the normalization each step of the sampling process along the lines of what's described in section 13.10 of BDA3 (pp 345-6).  (The soft constraint is similar but involves integrating a different function.)

Does this make sense?

[Bob Carpenter]

Yes, I think that's exactly what's going on. And I think that's
what happens in the Dirichlet normalization to let it operate over
simplexes (K dimensional simplex involves only K-1 dimensions
underlyingly).

I think if you declare the appropriate constraints on your x
parameters (or transform them there from raw), then you can just
put the regular generative density on x and everything will work out.
You should try that with simulated data to confirm

BUT (and this is a biggy), my understanding of the intricacies of
conditional densities and the geometry of lower-dimensional submanifolds
is sketchy at best.

I still never understood why constraining a vector to be ordered,
then giving each component a normal gives you the same distribution
as giving a bunch of unconstrained vector components normal priors
and then sorting,

parameteters
ordered_vector[4] x;
model
x ~ normal(0, 1);

VS

parameters
vector[4] x_raw;
transformed parameters
ordered_vector[4] x;
x <- sort(x_raw);
model
x_raw ~ normal(0, 1);

- Bob

[Bob Carpenter]

(Sorry if this is a dupe --- it was sitting unsent in my outbox.)

That makes it much clearer what's going on in terms of parameters
and the generative process.

So I think your (1) is on the right track:

> 1) Express x[i, 1] in terms of x[i, 2] and x[i, 3], e.g.
>
> real diff;
> for (i in 1:nsample) {
> x[i, 2] ~ lognormal(mu_2, sigma_2);
> x[i, 3] ~ lognormal(mu_3, sigma_3);
> diff <- (y[i] - a[i, 2] x[i, 2] - a[i, 3] x[i, 3]) / a[i, 1];
> diff ~ lognormal(mu_1, sigma_1);
> }

The main issue you're going to have with this given that
y and a are known but x is not, is that you have an implicit
constraint

y[i] > a[i,2] * x[i,2] + a[i,3] * x[i,3]

that follows from (diff > 0) and (a > 0).

Because Stan doesn't let you custom constrain each member of a
matrix like x, the only recourse I can see is to write your own
custom transforms and apply the appropriate Jacobians.

That is, constrain

0 < x[i,2] < y[i] / a[i,2]

0 < x[i,3] < (y[i] - a[i,2] * x[i,2]) / a[i,3]

That means declaring x_raw unconstrained, then inverse logit
and scale transforming it and applying the Jacobians.

Then you're in a position to use the diff formulation. But it's
a lot of work. Maybe someone else will see a clever way to encode
this.

- Bob

[Daniel Lakeland]

On Tuesday, September 15, 2015 at 2:39:29 PM UTC-7, Bob Carpenter wrote:

I would not recommend 2, because I don't think you can easily fake
a hard constraint with a soft constraint.

The problem for (1) isn't so much that you need a Jacobian with this transform

>   x[i,1] <- (y[i] - a[i, 2] x[i, 2] - a[i, 3] x[i, 3]) / a[i, 1];

but rather that it's that you're trying to provide three distributions for
x[i,1:3] on what are effectively two variables, x[i,2:3] and a funciton of them.
So you need to work out what the implications are.  It's like you're doing
something like this:

  y[1] ~ foo(...);
  y[2] ~ foo(...);
  y[1] + y[2] ~ foo(...);

It's much easier to work with a model conceptually if you can imagine
the generative process leading to the data.

Perhaps, but I believe that the above type of model is a perfectly valid declarative model. After some recent discussions here I posted to my blog a variety of posts on declarative vs generative bayesian models.

http://models.street-artists.org/?s=declarative

I'd be interested to hear your thoughts Bob, as you may have some experience with declarative vs imperative programming, and I think the analogy is quite good.

I'd be very interested to hear what a run of Geoff's model declared as in his (1) does. The big issue as I see it is only that he may have difficulty with sampling because of the possibility for diff to be negative and so samples to be rejected as falling off the edge of the world... Depending on his data, this might only be a theoretical problem and not a practical one.

As for the need for a jacobian with the diff statement, I've managed to convince myself that no jacobian is needed in this statement, as it's a statement of the probability associated to a transformed data point conditional on the values of the parameters. It's not even sensible to talk about the jacobian because the transform takes one data value and 5 parameters to a new "transformed data" value. Presumably if the parameters are correct then the transformed data value is in the high probability region of the lognormal, and if they're incorrect then the transformed value is NOT in the high probability region, and so the probability of those parameters is highly downweighted, thereby giving you the posterior inference you want.

After drinking the full blown Cox Theorem version of Bayesianism, it seems one is pretty free to use declarative probabilistic statements in modeling.

[Bob Carpenter]

Andrew and Michael and Ben are the three you really want on this case.

As to the Jacobian calculation, the function maps (x[i,2], x[i,3], y[i])
to (x[i,2], x[i,3], x[i,1]), so you can just work it out (I put it in
that order to make sure the Jacobian's triangular and has an easy determinant).
I think you're right here that the Jacobian term drops out in Stan (but
you should check that --- my Jacobian calcs are error prone without extensive
unit testing).

As to the blog, I don't believe there are "declarative" programming
languages in practice. In order to program Prolog well (as outlined,
for instance, in O'Keefe's awesome book, The Craft of Prolog), you have
to know how it's evaluated. If you don't, you'll get infinite loops at
worst, and inefficient programs at best. Order matters in practice, if not in theory,
because evaluation is depth first. Same for SQL --- you need to show the plan
in order to evaluate how good a program you get from your code.

The proof for Prolog is really simple. You code up an array with
binary search access --- it's where the log term comes from. In practice,
though, things are much worse in that some algorithms are nearly impossible
to convert to even moderately efficient Prolog (all that theory ignores
important little facts like memory locality and garbage collection overhead).

This used to be my bread and butter --- computation theory and logic programming :-)
I even wrote a book on the topic.

I was expecting the analogy to be that the generative model was declarative.
The truly declarative approach here, in my opinion, is laying out the joint
density. The generative approach is one way to do that, the imperative approach
is another, but in both cases, at least in Stan, all you're doing is defining
the joint density.

I put an edited form of these comments on the blog, or I should say, put them
into the moderation queue.

- Bob

[Daniel Lakeland]

On Wednesday, September 23, 2015 at 3:46:37 PM UTC-7, Bob Carpenter wrote:

Andrew and Michael and Ben are the three you really want on this case.

As to the Jacobian calculation, the function maps (x[i,2], x[i,3], y[i])
to (x[i,2], x[i,3], x[i,1]), so you can just work it out (I put it in
that order to make sure the Jacobian's triangular and has an easy determinant).  
I think you're right here that the Jacobian term drops out in Stan (but
you should check that --- my Jacobian calcs are error prone without extensive
unit testing).

I see that I was confused in thinking the x values were data.

still, the statement:

diff ~ lognormal(foo,bar)

is a statement about x[i,3] because x[i,3] *is* diff, so there is no transformation, only a deterministic connection between the value of x[i,3] and some other things that are known conditionally. If you were going to do

log(diff) ~ normal(...) for example, you'd need a jacobian correction.

at least, that's how I read it. I admit to being a little distracted at the moment, so hopefully Geoff will delve deeper to confirm as needed.

[Daniel Lakeland]

On Wednesday, September 23, 2015 at 5:26:37 PM UTC-7, Daniel Lakeland wrote:

On Wednesday, September 23, 2015 at 3:46:37 PM UTC-7, Bob Carpenter wrote:

Andrew and Michael and Ben are the three you really want on this case.

As to the Jacobian calculation, the function maps (x[i,2], x[i,3], y[i])
to (x[i,2], x[i,3], x[i,1]), so you can just work it out (I put it in
that order to make sure the Jacobian's triangular and has an easy determinant).  
I think you're right here that the Jacobian term drops out in Stan (but
you should check that --- my Jacobian calcs are error prone without extensive
unit testing).

I see that I was confused in thinking the x values were data.

still, the statement:

diff ~ lognormal(foo,bar)

is a statement about x[i,3]....

SORRY: x[i,1], typo.

[Daniel Lakeland]

On Wednesday, September 23, 2015 at 5:26:37 PM UTC-7, Daniel Lakeland wrote:

at least, that's how I read it. I admit to being a little distracted at the moment, so hopefully Geoff will delve deeper to confirm as needed.

My small child distraction is itself distracted, so I thought I might help clarify my thinking here:

If you know the value of x[i,2], x[i,3], and the a[i,j] values etc, then you can calculate x[i,1] and it's exactly equal to "diff".

Taking the full Cox/Bayes view of bayesian inference, where the probability values are relative credibility measures, you can then ask:

How credible is it that the x[i,2], x[i,3], and a[i,j] values really are what the sampler is telling me they currently are, given the observed y value? Well, it's more or less credible based on not what the value of y is, but what the value of "diff" is, which is a function of y and the other parameters. And the relative credibility of the diff values is given by the lognormal density function, hence

diff ~ lognormal(...)

in the end, as Bob says, you wind up defining a density over the parameters, and you can work out what it is, but it's somewhat more complicated than usual.

So, although it's not *generative* because it doesn't specify that "y" is a sample from some surrogate mathematical probability density (ie. it doesn't have y ~ foo_density(...)), it is a meaningful bayesian model in the same way that ABC models are meaningful bayesian models, even though the accept/reject criterion is not a generative one either.

[Michael Betancourt]

Whaaaaaaaaa?

I can’t follow this discussion at all, but let me try to clarify a few points that have come up.

First — the Cox/Jaynes perspective has absolutely no impact on modeling. It’s just an (incomplete)

motivation for why we can always place inductive learning into a probabilistic framework.  For more

details on the incompleteness of Cox/Jaynes in general state spaces see http://arxiv.org/abs/1507.06597.

Second — there’s nothing “generative” or “imperative” about our posterior distributions.  If

we had an imperative description of the posterior then we wouldn’t need HMC — we’d just

use ancestral sampling to get exact samples!  It is far more accurate to say that posterior

distributions in Stan are “declarative” in that the user declares a posterior (via a density)

and the code draws samples without the user having to specify how samples are drawn.

This is _always_ how Stan works.

The confusion arises because we advocate thinking generatively when building up a

_likelihood_.  Users typically have strong information about how the data were

generated, so this generative perspective makes it easier to build up robust likelihoods

representative of the true underlying process.  But, again, in Stan you don’t have to

specify a target distribution as prior + generative likelihood.  You can specify whatever

you want provided it is a self-consistent probability density.  We just find that prior +

generative likelihood is the easiest way to guarantee such a valid density.

Finally — another obfuscating aspect here is the difference between probability distributions

(measures) verses probability densities.  Densities are a parameterization-dependent

representation of a distribution, and the Jacobian comes around when comparing densities

from different parameterizations.  Jacobians arise in this example because you’re trying

to specify densities for parameterizations that are inconsistent with the parameterizations

assumed by the Stan probability densities.  Ultimately, if x ~ some_distribution then we

don’t know the distribution for f(x) so you have to define that yourself.  If you get it wrong

then you’ve specified a density for a different probability distribution!

[Daniel Lakeland]

On Thursday, September 24, 2015 at 1:46:25 AM UTC-7, Michael Betancourt wrote:

Whaaaaaaaaa?

I can’t follow this discussion at all, but let me try to clarify a few points that have come up.

First — the Cox/Jaynes perspective has absolutely no impact on modeling.

Well, I see we are in complete disagreement here. In particular the Cox perspective gives motivation for the MEANING of a distribution, and so it has EVERYTHING to do with modeling, since modeling is about creating a meaningful connection between the world and some mathematical statements. The Cox perspective tells you what you can put a distribution on (basically anything quantifiable), and what that distribution means (a degree of plausibility, not a relative frequency of occurrence, in particular, we need not match the empirical frequency of data points to the distribution function).

On the other hand, I completely agree with you when it comes to implementing a sampler. After the modeling, the mathematics is the mathematics. And Stan's whole purpose is to create a machine that explicitly connects the Bayesian plausibility measures to frequency measures in the output of Stan. The "problem" here is that Stan makes it so easy to connect the model to the code that the difference is sometimes hard to see  (that's a compliment to you all!)

>Second — there’s nothing “generative” or “imperative” about our posterior distributions.

There are however "generative" vs "declarative" methods of modeling. This discussion is about the step before you write down a Stan model, or rather, because Stan makes it so easy to translate ideas into code, it winds up looking more like interpreting the meaning (in terms of what it implies for the world) of Stan statements.

The confusion arises because we advocate thinking generatively when building up a

_likelihood_.  Users typically have strong information about how the data were

generated, so this generative perspective makes it easier to build up robust likelihoods

representative of the true underlying process.  But, again, in Stan you don’t have to

specify a target distribution as prior + generative likelihood.  You can specify whatever

you want provided it is a self-consistent probability density.

Yes, this is exactly the case. Here I'm arguing thinking about Stan statements as declarative statements which are valid Bayesian modeling statements (they are meaningful for describing the real-world process) even though they **don't define a likelihood** (at least, they don't define a likelihood based on data coming from an iid sampling process). They do, however, define a weighting function for the prior which acts like a likelihood (it downweights the parameter values which are inconsistent with the probabilistic statement of fact). From the Cox perspective, this is a perfectly natural thing to do, whereas from other perspectives of probability, it's meaningless because probability = frequency of occurrence, so you'd have to write statements of the form data ~ my_distribution because that's got a meaningful connection to frequencies.

In this case, Geoff's model (1) may well run just fine even though the statement

diff ~ lognormal(...)

doesn't define the usual type of iid likelihood.

It's worth coding the model and trying it out before going down other roads, because the only thing "wrong" with it is that it might fall off the edge of the world, when diff becomes negative. but that might not happen given his data and priors.

I hope that clarifies my position here. And please note: I'm just trying to help Geoff and others get something both meaningful and useful.

[Tamas Papp]

On Thu, Sep 24 2015, Daniel Lakeland <dlak...@street-artists.org> wrote:

> The confusion arises because we advocate thinking generatively when
> building up a _likelihood_. Users typically have strong information
> about how the data were generated, so this generative perspective
> makes it easier to build up robust likelihoods representative of the
> true underlying process. But, again, in Stan you don’t have to

The only time I have strong information about how my data was generated
is when I simulate said data myself for testing some code.

I don't know about other fields, but in social sciences such notions of
"strong information" about the DGP are rapidly shipwrecked by the
laughter of the gods.

Best,

Tamas

[Daniel Lakeland]

Even when you do have strong information, as Geoff may have, it is sometimes strong information on an expression which involves data and parameters together (as diff does here) and that is still information we can use in a legitimate way (that is, in a way which has a meaningful interpretation and which helps define a bayesian posterior even in the absence of an iid sampling story). In this particular case, in the presence of y the x values are not independent, but we can still say that plausible values for "diff" are given by the high probability region of a lognormal.


 

[Michael Betancourt]

> Well, I see we are in complete disagreement here. In particular the Cox perspective gives motivation for the MEANING of a distribution, and so it has EVERYTHING to do with modeling, since modeling is about creating a meaningful connection between the world and some mathematical statements. The Cox perspective tells you what you can put a distribution on (basically anything quantifiable), and what that distribution means (a degree of plausibility, not a relative frequency of occurrence, in particular, we need not match the empirical frequency of data points to the distribution function).

I didn’t say anything otherwise. You can assign whatever distribution you want, but
how you _manipulate_ that distribution is highly constrained by the rules of probability.

> Yes, this is exactly the case. Here I'm arguing thinking about Stan statements as declarative statements which are valid Bayesian modeling statements (they are meaningful for describing the real-world process) even though they **don't define a likelihood** (at least, they don't define a likelihood based on data coming from an iid sampling process). They do, however, define a weighting function for the prior which acts like a likelihood (it downweights the parameter values which are inconsistent with the probabilistic statement of fact). From the Cox perspective, this is a perfectly natural thing to do,

No. No it’s not.

Cox-Jaynes tells you that uncertainty of logical statements is equivalent to Komologorov
measure theory (up to the aforementioned technical details that are still somewhat
unresolved). At that point any learning has to be consistent with those rules of measure
theory and there is a huge wealth of literature showing that Bayes theorem (prior + likelihood)
is the only self-consistent means of updating an initial distribution with data. The road
to alternative parameter updates is paved with fuzzy sets.

You can specify whatever distribution you want in Stan and provided it admits a smooth
density with respect to the Lebesgue measure corresponding to the given parameterization
then Stan will draw samples (hopefully successfully). But unless that posterior arises from
a proper (not necessarily unique) prior + likelihood breakdown then it will not be a well-defined
inference procedure. Again, we don’t enforce this but it also means that you’re on your own
when it comes to interpreting the results.

[Michael Betancourt]

> The only time I have strong information about how my data was generated
> is when I simulate said data myself for testing some code.
>
> I don't know about other fields, but in social sciences such notions of
> "strong information" about the DGP are rapidly shipwrecked by the
> laughter of the gods.

Strong is relative here. In physics it can be extremely detailed, for
example the transition from incident particle to EM shower to photons to
photodetector current to op amp voltage to readout. Social sciences
may not be anywhere near as detailed, but considering a generative
narrative for, say, poll non-response or population exchangeability is
still the best way to identify a robust likelihood.

I’ll side eye any analysis claiming a generative model from quantum
mechanics to neuronal activity to human decision to poll response
too, but I’ll also suspect any analysis that uses a default likelihood
with no further consideration.

[Bob Carpenter]

I really like this part:

> ... posterior

> distributions in Stan are “declarative” in that the user declares a posterior (via a density)
> and the code draws samples without the user having to specify how samples are drawn.
> This is _always_ how Stan works.
>

> ... we advocate thinking generatively when building up a

> _likelihood_. Users typically have strong information about how the data were
> generated, so this generative perspective makes it easier to build up robust likelihoods
> representative of the true underlying process.

For me, it's this second step of building up reasonable likelihoods
(and priors) and then being able to do inference with them that makes Bayesian
modeling so attractive. Being able to account for how the data came
about both in terms of the underlying process (science) and the
measurement is critical. Basically what Andrew outlines in the first
bit of the first chapter of Bayesian Data Analysis.

> But, again, in Stan you don’t have to
> specify a target distribution as prior + generative likelihood. You can specify whatever
> you want provided it is a self-consistent probability density. We just find that prior +
> generative likelihood is the easiest way to guarantee such a valid density.

We need to somehow make that clearer. I've been trying
in talks to emphasize how the sampling notation is just
a syntactic shorthand.

- Bob

[Bob Carpenter]

Maybe we mean something different by "strong" here. From Andrew
et al.'s work on voting, there are modeling decisions such as which
demographic information to include (age, sex, education, income,
state of residence, etc.), where to include multilevel factors (effect
of income might vary by state of residence and have a correlated slope
and intercept), when to include time-series data (maybe voting varies
over time based on what politics was like when you came of age in
addition to the age trend itself), when to include interactions (interaction
of income and education), and so on.

Another example is Aki's take on the birthday problem (what calendar day
kids are born on). It's largely a social issue for humans, so Aki took
into account seasonal effects, day of week effects, special holiday effects,
a general historical trend, and so on).

In both cases, this isn't so much generative as it is regression based
(what a machine learning person would call "discriminative" in contrast
to "generative", to which Andrew objects for reasons explained in BDA [short
answer is that you can and should think of the whole model as generative
with the predictor model marginalized out due to an independence assumption]).

The main motivation for building Stan was that there were such intuitions
that Andrew et al. had that couldn't be fit in other software.
(And they're still not super easy to fit in Stan, but it's possible.)

- Bob




>
> Best,
>
> Tamas

[Daniel Lakeland]

On Thursday, September 24, 2015 at 8:50:33 AM UTC-7, Michael Betancourt wrote:

No.  No it’s not.

Cox-Jaynes tells you that uncertainty of logical statements is equivalent to Komologorov
measure theory (up to the aforementioned technical details that are still somewhat
unresolved).  At that point any learning has to be consistent with those rules of measure
theory and there is a huge wealth of literature showing that Bayes theorem (prior + likelihood)
is the only self-consistent means of updating an initial distribution with data.

I disagree, ABC is a hugely effective procedure in some circumstances, but it doesn't involve a likelihood in which data is the result of IID draws from a distribution, yet the typical accept-reject method can be thought of as updating a prior with the additional probabilistic information assigned as:

distance(summary(fake_data(parameters)),summary(real_data)) ~ uniform(0,epsilon)

for some complicated function "summary" and some complicated function "fake_data" (often a black box physical simulator) and some choice of distance function.

After the model is specified, measure theory takes over. So if you're interested in using the above ABC statement, but you want to compute with a slightly different distance function that's easier to calculate, you'll need to jacobian transform the uniform density for example...

But before the model is specified we're talking a modeling process, there is rarely one true model. So what's important is how we choose and interpret models. One GREAT way to choose models is to use a generative formalism to describe them, but it's not the only way.

I think to avoid talking past each other it's probably important here to define what is meant by a likelihood, and the usual definition is p(data | parameters) and we have to assign that, it's not an automatic thing you can prove from first principles. And, there's no reason to think that you need to say express that in terms of IID sampling, or restrict yourself to the data space. why not work in some kind of complicated distance between the summaries of the generated simulations and the data?

p(distance(summary(fake_data(parameters)),summary(real_data)) | parameters) ??? if you know the answer to what the plausible values are for that expression then you can define a density over that distance which expresses the plausibility, and get a proper (proper in a Cox/Jaynes interpretation) posterior for parameters.

Perhaps you'd just include this in "likelihoods" so we may be talking past each other based on definition of what is meant by "likelihood".

Or perhaps you disagree that assigning a distribution to distance(summary(fake_data(parameters)),summary(real_data)) is a valid thing to do. I would be very interested to clarify.

[Daniel Lakeland]

On Thursday, September 24, 2015 at 1:08:01 PM UTC-7, Daniel Lakeland wrote:

On Thursday, September 24, 2015 at 8:50:33 AM UTC-7, Michael Betancourt wrote:

No.  No it’s not.

Cox-Jaynes tells you that uncertainty of logical statements is equivalent to Komologorov
measure theory (up to the aforementioned technical details that are still somewhat
unresolved).  At that point any learning has to be consistent with those rules of measure
theory and there is a huge wealth of literature showing that Bayes theorem (prior + likelihood)
is the only self-consistent means of updating an initial distribution with data.

Another way to put my point is that measure theory doesn't get to tell you what you are going to model. It only gets to tell you *after* you specify the model, what you can do with it while keeping it invariant and logically consistent. If you have some data that you don't know much about, but you know a lot about some function of that data conditional on some values of some parameters, you can model that function.

That's exactly what Geoff's model (1) does. it assigns plausibilities to an expression "diff" which according to him under his modeling assumptions and conditional on the proper values of the parameters, looks like a number pulled from a lognormal distribution. This doesn't violate any rules of logic or measure theory, it is just a choice to model something other than the observed data point. Since the observed data point is always a function of some other stuff anyway (position of a needle on a gauge is a function of a current and a spring constant and some geometry and whatnot) this should be pretty obvious.

[Bob Carpenter]

This mailing list and especially this thread are probably
not the best venue for discussing the theory of probability
modeling. Let me just sum up what Stan does so that doesn't
get lost in all of this discussion.

Stan only requires you to define a the log of p(theta|y)
where p(theta|y) is your posterior up to a proportion.
The requirements are that it be differentiable in order
for HMC to work, have support on its declared constraints
for Stan to work, be defined in our language, so that we
can compute derivatives, and to be proper (finite integral).
Sometimes even with all this, the geometry will be too difficult
to adapt to or sample from, which will show up in the
convergence diagnostics.

You can define p(theta|y) any way you want. A Bayesian
would often define the posterior in terms of a likelihood
and prior via Bayes's rule as

p(theta | y) propto p(y | theta) p(theta)

For a Bayesian, nothing requires p(y | theta) to involve
a vector of data y that is i.i.d. It can be anything.

If you prefer penalized likelihoods and point estimation, you can think
of a model as defining a log penalized likelihood, as in

LL_penalized(y | theta) = log p(y | theta) + pen_fun(theta)

where pen_fun(theta) is some penalty function. Then the optimizers
will maximize this to find the penalized MLE.

- Bob

[Daniel Lakeland]

On Thursday, September 24, 2015 at 2:05:39 PM UTC-7, Bob Carpenter wrote:

This mailing list and especially this thread are probably
not the best venue for discussing the theory of probability
modeling.  Let me just sum up what Stan does so that doesn't
get lost in all of this discussion.

Thanks Bob, I think that was a good summary.

I hope my enthusiasm isn't considered noise by the list readers. In particular, I really did just want to help Geoff, because he had formulated a model in a way very similar to something I had tried to do recently and I had enough confusion over it that I went out and wrote a bunch of blog posts to discuss the issue, got a little feedback, and managed to convince myself that a model like Geoff's (1) is perfectly valid and interpretable under some interpretations of applied probability theory.

So, all trying to be in the spirit of helping people with modeling.

[Geoff Davis]

While I really like the declarative model building in Stan, it had me thinking about the problem in the wrong way.  I've enjoyed the philosophical discussion, and some of the comments spurred me to look at the problem a bit differently.  I think I've figured out what I need to do to move forward (but alas it involves writing my own MCMC code rather than using Stan).  Thanks, all!

[Bob Carpenter]

I still don't think Stan is "declarative" in any useful
sense. It uses imperative programming constructs (loops,
conditionals, local variables, etc.), and uses them to define
a log density, which it can then sample from. That's really
all Stan's doing.

If you want to sample from that model using custom MCMC code,
you'll still need to be careful with the constraints on the
x parameters imposed by your model:

y[i] ~ lognormal(log(a[i, 1] * x[i, 1] + a[i, 2] * x[i, 2] + a[i, 3] * x[i, 3]), epsilon);

And you still need to respect any required Jacobians, because
that's just probability theory. If the a > 0 and the x > 0
by construction, there shouldn't be a problem fitting this in Stan.

Good luck!

- Bob

[Daniel Lakeland]

On Wednesday, September 30, 2015 at 1:41:07 PM UTC-7, Geoff Davis wrote:

While I really like the declarative model building in Stan, it had me thinking about the problem in the wrong way.  I've enjoyed the philosophical discussion, and some of the comments spurred me to look at the problem a bit differently.  I think I've figured out what I need to do to move forward (but alas it involves writing my own MCMC code rather than using Stan).  Thanks, all!

Geoff, I did some thinking about this problem in comments over on my blog. But I'll give a summary over here:

let's ignore the weights for the moment and say that y is the sum of three variables

y = a + b + c

you declare

a ~ lognormal(...)
b ~ lognormal(...)

and then you want a likelihood p(y | a,b)

but c ~ lognormal(...) as well, and p(y | a,b,c) = delta_function(y-a-b-c).

so, to get p(y|a,b) you can marginalize away c...

p(y| a,b) = integrate(p_c(c) * p(y | a,b,c),dc) = integrate(p_c(c) * delta_function(y-a-b-c), dc) = p_c(y-a-b)

so:

y-a-b ~ lognormal(pc_parameters)

this is exactly the logic behind your model (1) and may well sample in Stan just fine if it stays away from the "edge of the world" since a,b,y-a-b all need to be positive.

I'd recommend trying to run your (1) and seeing what happens.

[Daniel Lakeland]

On Wednesday, September 30, 2015 at 2:24:55 PM UTC-7, Bob Carpenter wrote:

I still don't think Stan is "declarative" in any useful
sense.  It uses imperative programming constructs (loops,
conditionals, local variables, etc.), and uses them to define
a log density, which it can then sample from.  That's really
all Stan's doing.

It's not declarative in a programming sense, but it is declarative in a modeling sense. You "declare" the density you want to sample from, and it figures out how to do the sampling. It's a little like declaring the facts you know to prolog and then it figures out how (depth-first-search) to find the valid solutions to your problem.

So, in Stan, you don't describe the technique for getting the samples, you just describe the thing you want to sample from.

However, in Stan, the way you declare the density is to use iterative programming techniques to describe the density computation, so that's different from BUGS/JAGS for example where you declare the structure of a directed graph or whatever.