Latent mixture with known parameters in one case?

[James Savage]

Hi all, 

Let's say I have a time series whose growth is either 0 or positive. Some values of the time series are

missing and so I want to use a state-space representation. 

The state model is: 

x(t) = x(t-1) with probability theta

x(t) = x(t-1) + r(t) with probability 1- theta

and the "growth bump" r(t) is always non-negative

r(t) ~ normal_{+}(0, sigma)

Another way of putting it is that r(t) is 0 with probability theta, and takes the above distribution with probability

1- theta. There's then a measurement model: 

observation(t) ~ normal(x(t), sd_measurement_error)

The question is how to implement the structure? At the moment, I'm hacking it, by simply stating that rather than

being precisely 0 in the first case, r(t) is small.  

model {

   // priors for theta and sigma

    target += log_sum_exp(log(theta) + normal_lpdf(r[t] | 0, 0.01), // r[t] very small rather than precisely 0

                                          log(1-theta) + normal_lpdf(r[t] | 0, sigma);

} 

Is there a better way of doing this? 

Cheers, 

Jim

[Daniel Lee]

Hi Jim,

I know in stats we overload the term "parameter" all the time. In Stan, "known parameters" are data, not parameters. Just wanted to clarify that terminology. Anything known is data, anything unknown is parameters.

I think you're complicating things a little bit. Just based on what you're describing, it should be:

// missing loops, other important stuff

  target += log_sum_exp(log(theta) + log(1),

                                          log(1-theta) + normal_lpdf(r[t] | 0, sigma);

Note: I wrote out log(1) explicitly to hopefully clarify that you're saying:

x(t) = x(t-1) with probability theta

When you write it, you should probably drop that term.

Daniel

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[James Savage]

Hey Daniel - 

If r(t) in the first state is equal to 0, then isn't that a Dirac delta function (ie. has a density of +infty, not 1?). 

Or perhaps I'm not understanding something. 

Jim

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Daniel Lee]

On Mon, Sep 26, 2016 at 9:40 AM, James Savage <jav...@gmail.com> wrote:

Hey Daniel - 

If r(t) in the first state is equal to 0, then isn't that a Dirac delta function (ie. has a density of +infty, not 1?). 

True statement.

 

Or perhaps I'm not understanding something. 

I think there are two ways of looking at this:

1) you've written x(t) = x(t-1) with probability theta. Conditional on theta (if we treated this as a latent discrete 0/1 variable), that's saying Pr(x(t) = x(t-1)) = 1.

2) for dirac delta (in a super hand-wavy way): density is +infy, but the probability is 1. It's discontinuous. It's the same as saying: Pr(x(t) = x(t-1)) = 1 and Pr(x(t) != x(t-1)) = 0. (which is the same as #1... we can't ever say x(t) != x(t-1))

 

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Michael Betancourt]

Think more generatively — r[t] is always there, it’s just doesn’t always

affect x.  For example,

for(t in 2:T) {

  real mu_1;

  real mu_2;

  r[t] ~ normal(0, sigma); // r[t] are parameters declared to be positive

  mu_1 = x[t - 1];

  mu_2 = x[t - 1] + r[t];

  target += log_sum_mix(theta, normal_lpdf(x[t] | mu_1, sigma_meas, normal_lpdf(x[t] | mu_2, sigma_meas);

}

[Daniel Lee]

Ah... so Michael's model is different than the one I described. And that's cause I'm not sure what Jim has described.

Mike's missed a couple of close paren to make that parse... this should be clearer:

 target += log_sum_mix(theta, normal_lpdf(x[t] | mu_1, sigma_meas), normal_lpdf(x[t] | mu_2, sigma_meas));

So Mike is thinking of there being sigma_meas contributing to both terms. This might be much more natural than what I described. This will probably describe your model better.

Daniel

}

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
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+unsubscribe@googlegroups.com.

[James Savage]

Mike - you little ripper. Works a treat (after fixing parentheses). Also, isn't it log_mix()?  

[Avraham Adler]

On Monday, September 26, 2016 at 11:16:20 AM UTC-4, Daniel Lee wrote:

Ah... so Michael's model is different than the one I described. And that's cause I'm not sure what Jim has described.

Mike's missed a couple of close paren to make that parse... this should be clearer:

 target += log_sum_mix(theta, normal_lpdf(x[t] | mu_1, sigma_meas), normal_lpdf(x[t] | mu_2, sigma_meas));

So Mike is thinking of there being sigma_meas contributing to both terms. This might be much more natural than what I described. This will probably describe your model better.

Daniel

On Mon, Sep 26, 2016 at 10:40 AM, Michael Betancourt <betan...@gmail.com> wrote:

Think more generatively — r[t] is always there, it’s just doesn’t always

affect x.  For example,

for(t in 2:T) {

  real mu_1;

  real mu_2;

  r[t] ~ normal(0, sigma); // r[t] are parameters declared to be positive

  mu_1 = x[t - 1];

  mu_2 = x[t - 1] + r[t];

  target += log_sum_mix(theta, normal_lpdf(x[t] | mu_1, sigma_meas, normal_lpdf(x[t] | mu_2, sigma_meas);

}

Perhaps I'm misunderstanding, but according to Michael's parameterization we have two simultaneous processes as it were, one always affected by r[t] and one never affected. Doesn't that means when that we "jump" from the non-affected to the affected regime, we get the benefit of all the compounded r[t] for all t? But if I understood the original question, there is really only one process which is constant theta of the time and has a single jump (1-theta) of the time.

More concretely, let mu_0 = 1, r[0] = 0.1, r[1] = 0.2, and r[2] = 0.3.

Let's now assume that we trigger theta at time 2.

According to my original understanding, the mu progression is 1, 1, 1.3. According to my understanding of Michael's model it is 1. 1. 1.6,

Avi

[Michael Betancourt]

Mike - you little ripper.

Australian slang is the weirdest.

Works a treat (after fixing parentheses). Also, isn't it log_mix()?  

If I used the right syntax then it would be too easy.

[Bob Carpenter]

Agree with Michael that you want to think generatively. What
you really have is a discrete parameter:

z(t) ~ bernoulli(theta);

which will be 1 if x(t) != x(t-1) [equivalently r(t) != 0]
and 0 otherwise.

Marginalizing out the z(t) is exponential in the number
of time steps because there are 2^T paths for t in 1:T.

I don't see how you can fit the model in less than exponential
time in the number of time steps in Stan.

- Bob

[Michael Betancourt]

Won’t it be linear because the z(t) are independent?

Also, the model we eventually used was constant
as one theta was used to turn on or off all of the z(t).
To get the linear version there would be one theta
for each t.

[Bob Carpenter]

I don't see how to do it linearly because of the mass
arising from not moving when z(t) = 0.

- Bob