Homec4science
Diffusion Sparse coupled logistic regression (master) Browse

Sparse_logistic_regression/LogisticRegressionFramework/Framework/Miscellaneous/RSCHMM_Compute_HMM.m
c5f5e4464253master

Last Change
Raw File

RSCHMM_Compute_HMM.m
View Options

function [Transitions,Starts,Gammas,MU,SIGMA] = RSCHMM_Compute_HMM(Data,epsilon,MU,SIGMA)
%%% OUTPUTS TO FILL AND RETURN
% Number of time points per subject (reasonable number)
n_DP = size(Data{1},2);
% Number of networks per subject
N = size(Data{1},1);
n_subjects = length(Data);
% Will contain all transition probabilities across networks
Transitions = cell(1,N);
% Will contain all start probabilities across networks
Starts = cell(1,N);
% Gamma Probabilities across networks
Gammas = cell(1,N);
%%% RUNNING THE PHMM ITSELF
% We want to run this for all different systems
for i = 1:N
% We initialize the log-likelihood at the minimal value, and same for the
% difference of LL across iterations (maximal at first)
LL_before = - realmax;
Delta_LL = realmax;
% Transition probabilities (between two states)
TRANS = 1/2*ones(2,2);
% Starting probability (no reason to favour one)
START = [1/2;1/2];
disp(['Estimating network ',num2str(i),'/',num2str(N)]);
% Runs the HMM as long as we are not converging
while Delta_LL > epsilon
% Expectation step
% We will compute the parameters for all sequences separately
for m = 1:n_subjects
% Dynamic programming computation of alpha
ALPHA{m}(:,1) = START.*Get_Pdf_2STATE(Data{m}(i,1),MU(:,i),SIGMA(:,i));
Scaling(1) = 1/(sum(ALPHA{m}(:,1))+realmin);
ALPHA{m}(:,1) = ALPHA{m}(:,1)*Scaling(1);
for t = 2:n_DP
ALPHA{m}(:,t) = (ALPHA{m}(:,t-1)' * TRANS)' .* Get_Pdf_2STATE(Data{m}(i,t),MU(:,i),SIGMA(:,i));
Scaling(t) = 1/(sum(ALPHA{m}(:,t))+realmin);
ALPHA{m}(:,t) = ALPHA{m}(:,t)*sum(Scaling(1:t));
end
% Dynamic programming computation of beta
BETA{m}(:,n_DP) = ones(2,1)*Scaling(n_DP);
for t = n_DP-1:-1:1
BETA{m}(:,t) = TRANS * (Get_Pdf_2STATE(Data{m}(i,t+1),MU(:,i),SIGMA(:,i)) .* BETA{m}(:,t+1));
BETA{m}(:,t) = min(BETA{m}(:,t)*sum(Scaling(t:n_DP)),realmax);
end
% ALPHA and BETA are arrays, with each cell a sequence (i.e. a
% subject in our case), and dimension n_hidden_states x
% n_timepoints
% Computation of the likelihood of the full data
LL_tmp(m) = sum(ALPHA{m}(:,n_DP).*BETA{m}(:,n_DP));
% Computation of gamma
GAMMA{m} = (ALPHA{m}.*BETA{m})./repmat(sum((ALPHA{m}.*BETA{m}),1),2,1);
% Computation of zeta
for t = 1:n_DP-1
% Dimension i x j
tmp_zeta = ( repmat(ALPHA{m}(:,t),1,2) .* TRANS .* repmat(Get_Pdf_2STATE(Data{m}(i,t+1),MU(:,i),SIGMA(:,i))',2,1) .* repmat(BETA{m}(:,t+1)',2,1) );
ZETA{m}(:,:,t) = tmp_zeta/sum(sum(tmp_zeta,2));
end
clear Scaling
end
% Full likelihood of the dataset
LL = sum(log(LL_tmp));
% Maximization step
% Probability to start somewhere
START = zeros(2,1);
for m = 1:n_subjects
START = START + GAMMA{m}(:,1);
end
START = START/n_subjects;
% Transitions
TRANS = zeros(2,2);
for m = 1:n_subjects
TRANS = sum(ZETA{m},3)./repmat(sum(GAMMA{m},2),1,2);
end
% Mean values
MU(:,i) = zeros(2,1);
denom_gamma = zeros(2,1);
for m = 1:n_subjects
for t = 1:n_DP
MU(:,i) = MU(:,i) + GAMMA{m}(:,t) .* Data{m}(i,t);
denom_gamma = denom_gamma + GAMMA{m}(:,t);
end
end
MU(:,i) = MU(:,i)./(denom_gamma+realmin);
% Variance
SIGMA(:,i) = zeros(2,1);
for m = 1:n_subjects
for t = 1:n_DP
SIGMA(:,i) = SIGMA(:,i) + GAMMA{m}(:,t) .* ((Data{m}(i,t)-MU(:,i)).^2);
end
end
SIGMA(:,i) = SIGMA(:,i)./(denom_gamma+realmin);
Delta_LL = abs(LL_before - LL);
LL_before = LL;
disp(['Log-likelihood: ',num2str(LL_before),' and Delta log-likelihood: ',num2str(Delta_LL)]);
end
Transitions{i} = TRANS;
Starts{i} = START;
Gammas{i} = GAMMA;
clear START
clear TRANS
clear Delta_LL
clear LL_before
clear LL
clear LL_tmp
clear ALPHA
clear BETA
clear GAMMA
clear ZETA
clear Scaling
end
end
c4science · Help