I’m working on inference on a HMM with cmdstan
and, following Bales’s example, I need to define the likelihood at all states to be used in the hmm_marginal
function to integrate out the hidden states.
The data is defined as row_vector[N] y;
and the parameters as
vector[Q] Mu;
real<lower=0.0> Sigma;
matrix[Q, N] log_omega;
then the likelihoods
for(n in 1:N) {
for(i in 1:Q){
log_omega[i, n] = normal_lpdf(y[n] | Mu[i], Sigma);
}
}
the number of states Q
and of observed points N
could be quite large so I was wondering if there’s a more efficient way of coding this step rather than a double loop.
thank you for your time!
In principle, yes, but you would need to write a custom likelihood function equivalent to the normal_lpdf, but which returns a vector. The existing function is vectorized, but it sums all the log likelihods and returns a scalar. This would work:
functions {
row_vector normal_lpdf_vec(row_vector y, real mu, real sigma) {
return -((y - mu)/sigma)^2 / 2 - log(sigma) - log2()/2 - log(pi())/2;
}
}
// other parts
transformed parameters {
// other parts
for(i in 1:Q){
log_omega[i, ] = normal_lpdf_vec(y | Mu[i], Sigma);
}
}
But I doubt that would have a significant effect, as loops in C++ are quite efficient. For the example in the link you posted, I updated the code with such a vectorized version and the times are identical. But that’s with small N and Q, so maybe with more it would work. Feel free to try it
PS: some of the operations in normal_lpdf_vec are the same across loops (because sigma is teh same and the constants), so you could take these out and only apply them to the final log_omega matrix. Maybe that speeds it up a bit more
2 Likes
Thanks for working on this, @Ven_Popov
I tried the vectorised function with Q~40 and N~5000 but, as you anticipated, there is no significant improvement.I’ll rethink the approach and see if I can come up with a better way of proceeding.
thanks again!
Davide
1 Like