mixture syntax

[Bob Carpenter]

Rob and I are going to tackle feature

add mixture support functions
https://github.com/stan-dev/stan/issues/970

as an example of how to code functions with autodiff.

*** QUESTION ***: what should the function be called?

The feature request proposes two functions:

1. binary mixture

mix(real lambda, real lp1, real lp2)
= log_sum_exp(log(lambda) + lp1, log1m(lambda) + lp2)

where lambda must be in (0,1).

and

2. general mixture

mix(vector theta, vector lp)
= log_sum_exp(log(theta) + lp)

where theta is a K-simplex and lp a K-vector.

This would let us write a normal 2-mixture as:

increment_log_prob(mix(lambda,
normal_log(y[n], mu[1], sigma[1]),
normal_log(y[n], mu[2], sigma[2])));

instead of what we have now:

increment_log_prob(
log_sum_exp(log(lambda) + normal_log(y[n], mu[1], sigma[1]),
log1m(lambda) + normal_log(y[n], mu[2], sigma[2])));

and it'll reduce memory and speed up gradients.

If we call the function mix_log, we could even write:

lambda ~ mix(normal_log(y[n], mu[1], sigma[1]),
normal_log(y[n], mu[2], sigma[2])));

but that seems a bit weird, because it'd be super-confusing to
write the whole thing out as

for (n in 1:N)
lambda ~ mix(normal_log(y[n], mu[1], sigma[1]),
normal_log(y[n], mu[2], sigma[2])));

even though it'd do the right thing with the log prob accumulator
and is a less verbose syntax than

for (n in 1:N)
increment_log_prob(mix_log(lambda,
normal_log(y[n], mu[1], sigma[1]),
normal_log(y[n], mu[2], sigma[2])));

What I'd like to have eventually is something vectorized like this:

y ~ finite_mixture<normal>(lambda, mu, sigma);

or some other syntax with a function-like argument "normal" that'd
expand to the above. This would be a bit tricky on the syntactic side
to get everything to line up, but certainly possible. It'd yield a new
kind of statement (or looked at from a user perspective, a new kind of
way for writing distributions down).

Or if we're going away from the sampling statement syntax, something
like:

increment_log_prob(finite_mixture<normal>(y, lambda, mu, sigma));

Then finite_mixture... is a special kind of expression.

Maybe we should just drop the "increment_"? And maybe capitalize?

LOG_PROB(finite_mixture<normal>(y, lambda, mu, sigma));

It'd then just be read as marking erms in the overall log probabilty function.

- Bob

[Aki Vehtari]

> What I'd like to have eventually is something vectorized like this:
>
> y ~ finite_mixture<normal>(lambda, mu, sigma);

I just wanted to say that this would be really cool!

Aki