Potential bug in hmmsinit.m

[HDhanis]

Hi,

I have been trying to implement HMM for fMRI data using STAN when last week a colleague recommended me this toolbox. For quick context I am interested in exploring how hidden states vary their transition probabilities with some population covariates (transition probabilities or any of its derivative properties), hence the implementation I have been trying to adapt is the pipeline proposed from the example "NeuroImage2021.m".

I am in one of my first attempts but I can't get past what seems to me a bug in hmmsinit.m.

To give you the order of events, this happens when running hmmmar for the first time. At some point it hits lines 192-194:

if isempty(options.Gamma) && isempty(options.hmm) % both unspecified
[hmm,info] = hmmsinit(data,T,options);
GammaInit = [];

These lines will call hmmsinit, because options.Gamma and options.hmm are not defined at this point. This is consistent with the examples of PNAS2017_fullpipeline, NeuroImage2021, and run_HMMMAR, where both options are undefined before the first call to hmmmar.

With hmmsinit.m being called a problem then arises in line 232:
subj_m_init(:,:,subj,assig(k)) = XG * Y(t,:);

"assig" is an optimal assignment matrix based on a cost matrix, hence assig will always be a matrix of ones and zeros. Since the line I mentioned is inside a for loop (line 230 for k = 1:K_i) that goes through assig, it will eventually always try to do something like subj_m_init(:,:,subj,0) which is not acceptable in matlab notation.

I had a look at the initialisation of subj_m_init (subj_m_init = zeros(npred,ndim,N,K);), and it seems to me that what should be happening in line 232 is that you are trying to do the matrix-assignment to the fourth-dimension index of the best cost-optimised state. If this interpretation is correct then line 232 should actually read:

subj_m_init(:,:,subj,find(assig(k, :))) = XG * Y(t,:);

Is my interpretation correct? I don't understand how noboby has come across this bug, because in theory it is not possible to pass this without fixing the assignment step. So I am thinking that outside the examples everybody must have been defining either options.Gamma or options.hmm? If that's a simple solution could you let me know what are their purpose so I know how to define them?

Thank you very much!
Herberto Dhanis

P.S: With regards to the error above, this problem is larger than what I described because all of this is inside a loop starting in line 126 for ii = 1:length(I), and the problem I described starts for ii >1. This is because for ii == 1, assig is actually just a vector rather than a diagonal matrix. I am testing a few different things to fix this script, but so far this seems t have solved the problem.

[vidaurre]

Hi Herberto,

First of all, many thanks for the very detailed diagnosis and explanation, this is appreciated.

I'm not been able to reproduce the problem actually. The following code for example runs fine. And to me, assig is always a vector with positions (not 0s or 1s).

I would also like to mention that we stopped supporting this version of the toolbox, in favour of the Python version, which code should also be cleaner:
https://github.com/vidaurre/glhmm
https://direct.mit.edu/imag/article/doi/10.1162/imag_a_00460/127499

And we'd very super happy to address any issue with that :)

X = randn(4000,3); T = 500*ones(8,1);
B = rand(3);
for i=1:length(T)
X((1:250)+sum(T(1:i-1)),:) = X((1:250)+sum(T(1:i-1)),:) * B;
for j=1:3
X((1:500)+sum(T(1:i-1)),j) = smooth(X((1:500)+sum(T(1:i-1)),j),10);
end
end
C = rand(3,10);
X = X * C;
options = struct();
options.K = 2;
options.Fs = 200;
options.pca = 0.95;
options.covtype = 'full'
options = struct();
options.K = 2;
options.Fs = 200;
options.covtype = 'full'
options.order = 0;
options.useParallel=0;
options.BIGNbatch = 4;
hmm = hmmmar(X,T,options);