function [Ps, Pss, logPx] = hmmposteriors(x,hmm)
  n = length(x);
  ns = size(hmm.Pt,1);
  
  % loga(i,t) = log P(state_t=i AND prefix of output is x_1..x_t)
  loga = logalpha(x,hmm); % forward propagation
  % logb(i,t) = log P(sufix of output is x_t+1..x_n | state_t=i)
  logb = logbeta(x,hmm);  % backward propagation
  
  % logPx = log P(complete output sequence is x)
  logPx = logsumexp(loga(:,n));

  % logPsx(i,t) = log P(state_t = i AND complete output sequence is x)
  logPsx = loga+logb;
  % Ps(i,t) = P(state_t = i | complete output sequence is x)
  Ps = condprobexp(logPsx);
  % (*) Ps = exp(logPsx - logPx)

  if nargout>1,
    % logbx(i,t) = log P(sufix of output is x_t..x_n | state_t=i)
    logbx  = log(hmm.Po(:,x))+logb;
    % logPssx(i,j,t) = log P(state_t = i AND state_t+1=j AND output is x)
    %    = log P(state_t = i AND output prefix is x_1..x_t)
    %     +log P(state_t+1 = j | state_t = i)
    %     +log P( output sufix is x_t+1..x_n | state_t+1 = j)
    %    = loga(i,t) + log(Pt(i,j)) + logbx(j,t+1)
    logPssx = op(op(reshape(loga(:,1:n-1),ns,1,n-1),...
		    '+',reshape(logbx(:,2:n),1,ns,n-1)), ...
		 '+',log(hmm.Pt));
    % (Taking logbx(:,2:n) shifts it to the left, therby aligning loga(:,t) 
    % with logbx(:,t+1). The reshapes align the loga(i,t) with
    % logPssx(i,:,t) and logbx(j,t+1) with logPssx(:,j:t))
    % Pss(i,j,t) = P(state_t=i, state_t+1=j | complete output is x)
    Pss = reshape(condprobexp(reshape(logPssx,ns*ns,n-1)),ns,ns,n-1);
    % (The reshape is needed since we want to condition on the sum of the first
    % two dimensions)
    % (*) Pss = exp(logPssx-logPx)
  end;
    
% (*) These should be enough, without any need for condprobexp. By using
% condprobexp we, in effect, calculate logPx for each t seperately. All
% these different calculations of logPx /should/ be identical, but because
% of numerical issues, they might deviate a bit from one another. By
% calculating the conditional probabilities for time t by dividing by logPx
% calculated at time t, we ensure that sum_i Ps(i,t) = sum_i,j Pss(i,j,t) 
% = 1.
