Daniel@0: function s = logsumexp(a, dim)
Daniel@0: % Returns log(sum(exp(a),dim)) while avoiding numerical underflow.
Daniel@0: % Default is dim = 1 (rows) or dim=2 for a row vector
Daniel@0: % logsumexp(a, 2) will sum across columns instead of rows
Daniel@0: 
Daniel@0: % Written by Tom Minka, modified by Kevin Murphy
Daniel@0: 
Daniel@0: if nargin < 2
Daniel@0:   dim = 1;
Daniel@0:   if ndims(a) <= 2 & size(a,1)==1
Daniel@0:     dim = 2;
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: % subtract the largest in each column
Daniel@0: [y, i] = max(a,[],dim);
Daniel@0: dims = ones(1,ndims(a));
Daniel@0: dims(dim) = size(a,dim);
Daniel@0: a = a - repmat(y, dims);
Daniel@0: s = y + log(sum(exp(a),dim));
Daniel@0: %i = find(~finite(y));
Daniel@0: %if ~isempty(i)
Daniel@0: %  s(i) = y(i);
Daniel@0: %end