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