DLM for Proportions/Multinomial Model

319 views
Skip to first unread message

Devin Caughey

unread,
Jan 16, 2014, 10:59:35 AM1/16/14
to stan-...@googlegroups.com
Hello all,

I have a question that is more statistical that technical, so I apologize if this is not appropriate for the Stan list. I am wondering about how best to parameterize a dynamic linear model with a multinomial observation equation, where the parameters of interest are the time-specific population proportions for the different categories. The approach I've tried so far is to apply logit transformation to the proportions and specify the Gaussian DLM on that scale. Here is the code:

data {
    int<lower=1> T; // number of time periods
    int<lower=1> K; // number of categories
    int<lower=0> counts[T, K];
}
parameters {
    simplex[K] props[T]; // period-specific category proportions
    real<lower=0> tau; // innovation standard deviation
}
model {
    for (t in 1:T) {
        if (t > 1) {
            for (k in 1:K) {
                logit(props[t, k]) ~ normal(logit(props[t - 1, k]), tau); // Gaussian on logit scale
            }
        }
        counts[t] ~ multinomial(props[t]); // Multinomial sampling model
    }
}

This model seems to work well in the simulations I've run, but I am curious whether there is a better way to approach this problem. If not, I'll just stick with this. Thanks for your help!

Devin


Bob Carpenter

unread,
Jan 16, 2014, 11:26:41 AM1/16/14
to stan-...@googlegroups.com
I'll leave the statistical modeling advice to others, but
I'd like to point out that your model isn't quite correct as formulated.
The problem is that you need to add the Jacobian of the inverse
logit transform to the total log probability.

That'd be

for (k in 1:K) {
logit(props[t, k]) ~ normal(logit(props[t - 1, k]), tau); // Gaussian on logit scale
increment_log_prob(log(inv_logit(props[t,k])) + log1m(inv_logit(props[t,k])));
}

Here the Jacobian adjustment is the log of the derivative of the
inverse transform,

log | d/dp inv_logit(p) |

= log( inv_logit(p) * (1 - inv_logit(p) )

= log(inv_logit(p)) + log1m(inv_logit(p))

where p = props[t,k]. We discuss the Jacobian for the logit in the
section of the manual on constraints and transforms (near the end).

A more typical formulation I've seen is to parameterize props itself on the
logit scale (i.e., unconstrained vector[K] here), and then use softmax
to transform that to the parameter to plug into the multinomial. Then you
don't need to fiddle with Jacobians.

BUT, and this is a big but for now, our current softmax implementation has a memory
leak. If your system's small, it probably won't matter.
We fixed it on the 2.1.1 hotfix branch on GitHub, but we haven't released the
code with the fix yet, as we had the same kind of memory leak in a few other
functions that we're going to fix. In fact, I'm doing that in the other window
I have open now :-)

You could also speed the model up a bit by vectorizing the call to normal, but I'd
wait until the model's all working to try that.

- 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.
> For more options, visit https://groups.google.com/groups/opt_out.

Devin Caughey

unread,
Jan 16, 2014, 11:36:15 AM1/16/14
to stan-...@googlegroups.com
Ah, thank you (again) Bob. I must admit that I haven't quite wrapped my mind around the role of Jacobians and the total log probability, which is something I'd better rectify. Thanks again.

Devin


You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

Bob Carpenter

unread,
Jan 16, 2014, 1:42:21 PM1/16/14
to stan-...@googlegroups.com
The chapter of the manual starts from scratch, so hopefully
it'll be reasonably easy to follow. It's
pretty easy to see if you draw the densities out and
think about probabilities in intervals.

- Bob
>>> stan-users+...@googlegroups.com <mailto:stan-users+...@googlegroups.com>.
>>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>> --
>> You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com
>> <mailto:stan-users+...@googlegroups.com>.

Devin Caughey

unread,
Jan 16, 2014, 2:59:38 PM1/16/14
to stan-...@googlegroups.com
Thanks Bob. Section 17 of the Stan manual is particularly helpful.

Just to follow up, I rewrote my earlier model in two different ways. The first is similar to the first in that it entails a change of variables and thus requires a Jacobian adjustment:

data {
    int<lower=1> T; // number of time periods
    int<lower=1> K; // number of categories
    int<lower=0> counts[T, K];
}
parameters {
    simplex[K] props[T]; // period-specific category proportions
    real<lower=0> tau; // innovation standard deviation
}
transformed parameters {
    vector[K] logit_props[T]; // logit of proportions
    for (t in 1:T) {
        for (k in 1:K) {
            logit_props[t, k] <- logit(props[t, k]);
            // Adjust for logit tranform by incrementing the log probability
            // by the log of the derivative of the inverse-logit
            increment_log_prob(log(inv_logit(props[t, k])) +
                               log1m(inv_logit(props[t, k]))); 
        }
    }
}
model {
    for (t in 1:T) {
        if (t > 1) {
            logit_props[t] ~ normal(logit_props[t - 1], tau); //vectorized
        }
        counts[t] ~ multinomial(props[t]);
    }
}

I also coded up a version with a variable transformation using softmax, which doesn't require an adjustment:

data {
    int<lower=1> T; // number of time periods
    int<lower=1> K; // number of categories
    int<lower=0> counts[T, K];
}
parameters {
    vector[K] logit_props[T]; // logit of proportions
    real<lower=0> tau; // innovation standard deviation
}
transformed parameters {
    simplex[K] props[T]; // period-specific category proportions
    for (t in 1:T) {
        props[t] <- softmax(logit_props[t]);
    }
}
model {
    for (t in 1:T) {
        if (t > 1) {
            logit_props[t] ~ normal(logit_props[t - 1], tau);
        }
        counts[t] ~ multinomial(props[t]);
    }
}

~Devin
>>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>> --
>> You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com

Bob Carpenter

unread,
Jan 17, 2014, 8:55:56 AM1/17/14
to stan-...@googlegroups.com
Sorry --- I just realized I gave you the Jacobian the wrong way around.
I'm used to going from unconstrained parameters to parameters,
which involves the inverse logit. You're going the other way,
so you need to increment with the log absolute derivative of
the logits, not the inverse logits, that is:

log | d/dp logit(p) |

= log | d/dp log(p/(1-p)) |

= log | d/dp log(p) - d/dp log(1-p) |

= log | 1/p + 1/(1-p) |

= log ( 1/(p * (1-p)) )

= - log (p * (1 - p))

= -[log(p) + log(1-p)]

= -[log(p) + log1m(p)]

I'd check my work on that derivative, too.

- Bob
> > On Jan 16, 2014, at 11:26 AM, Bob Carpenter <ca...@alias-i.com <javascript:> <mailto:ca...@alias-i.com
> >>> stan-users+...@googlegroups.com <javascript:> <mailto:stan-users+...@googlegroups.com <javascript:>>.
> >>> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.
> >>
> >> --
> >> You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list"
> group.
> >> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe
> <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe>.
> >> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com <javascript:>
> >> <mailto:stan-users+...@googlegroups.com <javascript:>>.
> >> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.
> >
> > --
> > 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 <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

Bob Carpenter

unread,
Jan 17, 2014, 9:02:46 AM1/17/14
to stan-...@googlegroups.com
I should've also added that the way to check all of this is
to simulate data according to the model and make sure the
Stan model can recover the parameters.

- Bob

Devin Caughey

unread,
Jan 17, 2014, 1:25:22 PM1/17/14
to stan-...@googlegroups.com

Hi Bob,

Thanks for the correction. Your math is correct. The updated code is:

    int<lower=1> T; // number of time periods
    int<lower=1> K; // number of categories
    int<lower=0> counts[T, K];
}
parameters {
    simplex[K] props[T]; // period-specific category proportions
    real<lower=0> tau; // innovation standard deviation
}
transformed parameters {
    vector[K] logit_props[T]; // logit of proportions
    for (t in 1:T) {
        for (k in 1:K) {
            logit_props[t, k] <- logit(props[t, k]);
            // Adjust for logit tranform by incrementing the log probability
            // by the log of the derivative of the logit function
            increment_log_prob(-(log(props[t, k]) + log1m(props[t, k])));
        }
    }
}
model {
    for (t in 1:T) {
        if (t > 1) {
            logit_props[t] ~ normal(logit_props[t - 1], tau);
        }
        counts[t] ~ multinomial(props[t]);
    }
}


Devin

>>      >>> stan-users+...@googlegroups.com <javascript:> <mailto:stan-users+unsub...@googlegroups.com <javascript:>>.
>>      >>> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.
>>      >>
>>      >> --
>>      >> You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list"
>>     group.
>>      >> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe
>>     <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe>.
>>      >> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com <javascript:>
>>      >> <mailto:stan-users+unsub...@googlegroups.com <javascript:>>.

Devin Caughey

unread,
Jan 17, 2014, 1:27:44 PM1/17/14
to stan-...@googlegroups.com
Code without the first line cut off:

data {
    int<lower=1> T; // number of time periods
    int<lower=1> K; // number of categories
    int<lower=0> counts[T, K];
}
parameters {
    simplex[K] props[T]; // period-specific category proportions
    real<lower=0> tau; // innovation standard deviation
}
transformed parameters {
    vector[K] logit_props[T]; // logit of proportions
    for (t in 1:T) {
        for (k in 1:K) {
            logit_props[t, k] <- logit(props[t, k]);
            // Adjust for logit tranform by incrementing the log probability

Devin Caughey

unread,
Jan 29, 2014, 10:30:52 PM1/29/14
to stan-...@googlegroups.com
One more follow-up on this thread: Upon further reflection, I came to the conclusion that it would be more direct and possibly more efficient to use a Dirichlet distribution for the dynamic component of this model. I therefore re-wrote the model as follows:

data {
    int<lower=1> T; // number of time periods
    int<lower=1> K; // number of categories
    int<lower=0> counts[T, K];
}
parameters {
    simplex[K] props[T]; // period-specific category proportions
    real<lower=1> priorN; // implied "sample size" of the previous period
}
model {
    for (t in 1:T) {
        if (t > 1) {
            props[t] ~ dirichlet(priorN * props[t - 1]); // transition equation
        }
        counts[t] ~ multinomial(props[t]);
    }
}

Note that the parameter "priorN" plays roughly the same role as the innovation standard deviation parameter "tau" in the earlier version of the model.

~Devin

Bob Carpenter

unread,
Jan 30, 2014, 10:17:34 AM1/30/14
to stan-...@googlegroups.com
Neat. Did it work as expected?

Both these models remind me of Blei and Lafferty's dynamic topic models, which are
also time-series multinomials, although modeled quite differently:

http://www.cs.cmu.edu/~lafferty/pub/dtm.pdf

Instead of modeling the evolution of the parameters for a Dirichlet-Multinomial,
they use a logistic normal that evolves the actual multinomial
parameter on the logit scale. That's more like your first model.

- Bob
> >> >>> stan-users+...@googlegroups.com <javascript:> <mailto:stan-users+...@googlegroups.com
> <javascript:>>.
> >> >>> For more options, visit https://groups.google.com/groups/opt_out
> <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out
> <https://groups.google.com/groups/opt_out>>.
> >> >>
> >> >> --
> >> >> You received this message because you are subscribed to a topic in the Google Groups "stan users
> mailing list"
> >> group.
> >> >> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe
> <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe>
> >> <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe
> <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe>>.
> >> >> To unsubscribe from this group and all its topics, send an email to
> stan-users+...@googlegroups.com <javascript:>
> >> >> <mailto:stan-users+...@googlegroups.com <javascript:>>.
> >> >> For more options, visit https://groups.google.com/groups/opt_out
> <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out
> <https://groups.google.com/groups/opt_out>>.
> >> >
> >> > --
> >> > 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 <javascript:>.
> >> > For more options, visit https://groups.google.com/groups/opt_out
> <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out

Devin Caughey

unread,
Jan 31, 2014, 7:26:21 PM1/31/14
to stan-...@googlegroups.com
Yes, the Dirichlet parameterization actually worked extremely well, at least for the initial application. When I tried applying the model to a problem where the proportions were very small, however, I started running into "vanishing density" problems, which I gather from previous posts is the consequence of getting too close to parameter constraints. I'm working on debugging them now, but I may re-post with more questions later.

Devin
>              >>      >>> stan-users+...@googlegroups.com <javascript:> <mailto:stan-users+unsub...@googlegroups.com
>             <javascript:>>.
>              >>      >>> For more options, visit https://groups.google.com/groups/opt_out
>             <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out
>             <https://groups.google.com/groups/opt_out>>.
>              >>      >>
>              >>      >> --
>              >>      >> You received this message because you are subscribed to a topic in the Google Groups "stan users
>             mailing list"
>              >>     group.
>              >>      >> To unsubscribe from this topic, visit
>             https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe
>             <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe>
>              >>     <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe
>             <https://groups.google.com/d/topic/stan-users/IpFRiQFFYps/unsubscribe>>.
>              >>      >> To unsubscribe from this group and all its topics, send an email to
>             stan-users+...@googlegroups.com <javascript:>
>              >>      >> <mailto:stan-users+unsub...@googlegroups.com <javascript:>>.
Reply all
Reply to author
Forward
0 new messages