function [ym,ysum,ns] = condexpvec(y,x,xrg)
	% [ym,ys] = condexp(y,x,xrg)
	% Computes ym(ii) = mean(y(:, x>=xrg(ii) & x < xrg(ii+1)) for every ii
	% additionally computes
	%          ys(ii) = sum(y(:, x>=xrg(ii) & x < xrg(ii+1)) for every ii
	% and      ns(ii) = sum(x>=xrg(ii) & x < xrg(ii+1))
	%
	% uses a fast algorithm based on the MATLAB builtin accumarray
	% supports conditionally averaging vector data
	% 
	
	x					= x(:);
	xrg				= xrg(:)';
	nx					= numel(x);
	sy					= size(y);
	assert(prod(sy(2:end)) == nx);
	
	n_edges			= length(xrg);
	n_bins			= n_edges-1;
	[~,bins]			= histc(x,xrg);

	ok					= bins > 0 & bins < n_edges;
	n_ok				= sum(ok);
	ns					= accumarray(bins(ok), ones(n_ok,1), [n_bins 1]);
	if n_ok == 0
		ysum			= zeros(sy(1), n_bins);
		ym				= nan(sy(1), n_bins);
		return;
	end

	is					= (1 : sy(1))'  * ones(1,n_ok);
	js					= ones(sy(1),1) * bins(ok)';
	ys					= y(:, ok);
	
	ysum				= accumarray({is(:) js(:)}, ys(:), [sy(1) n_bins]);
	ym					= bsxfun(@rdivide, ysum, shiftdim(ns,-1));
end