function post = mixpost(net,y_exp,y_gate,t,sigma)
%mixpost Calculation of the posterior probabilities for a mixture of experts (or gated experts) model
%
%        POST = MIXPOST(NET,Y_EXP,Y_GATE,T) computes the posterior
%        probabilities. 
%        It takes the matrices of the output vectors of the experts 
%        Y_EXP, the matrix of the output vectors of the gating network 
%        Y_GATE, and the target vectors T. The choice of the conditional 
%        density of the experts corresponds to their output unit activation 
%        functions. POST has the same dimensions as Y_GATE (the values 
%        of which can be considered as the prior probabilities).
%
%        POST = MIXPOST(NET,Y_EXP,Y_GATE,T,SIGMA) if the conditional
%        density for the experts is an isotropic Gaussian (with one 
%        variance parameter for each expert), SIGMA contains these
%        variances.
%
%        See also
%        mix, mixem, mixfwd, gmmpost
%
%Copyright: Perry Moerland

ndata = size(t, 1);
nexp  = net.nexp;

%initialize posteriors
post = zeros(ndata,nexp);

switch net.actfn
  
  case 'linear'         % Linear outputs  
    
	for i=1:nexp
         % Calculate squared norm matrix, of dimension (ndata, 1) 
         n1 = t-y_exp(:,:,i);
         n2 = sum(n1.^2,2);
	 
	 switch net.tag
	   
	   case {'standard','kernel','mfa','mpca'}
	    % no variance parameters
            normal = (2*pi) .^ (net.nout/2);
            % compute the activations
            a = exp(-(n2/2))./ normal;
	    
  	   case {'variance','ker&var','mfa&var','mpca&var'}
	    % isotropic Gaussian 
	    normal = (2*pi*sigma(i)) .^ (net.nout/2);
            % compute the activations
            a = exp(-(n2/(2*sigma(i))))./ normal;
	    
	 end;
     
         post(:,i) = y_gate(:,i).*a;
        end;
 
  case 'logistic'       % Logistic outputs
    
    for i=1:nexp
     % Calculate cross-entropy norm matrix of dimension (ndata, 1) 
     n1 = (y_exp(:,:,i).^t).*((1-y_exp(:,:,i)).^(1-t));
     
     post(:,i) = y_gate(:,i).*n1;
    end;
   
  case 'softmax'        % Softmax outputs   
    
    for i=1:nexp
     % Calculate multinomial norm matrix of dimension (ndata, 1) 
     n1 = y_exp(:,:,i).^t;
     n2 = prod(n1,2);
     
     post(:,i) = y_gate(:,i).*n2;
    end;
   
end;

s = sum(post,2);
% Set any zeros to one before dividing
s = s + (s==0);
post = diag(1./s)*post;
