
function y = normpdf_(x,mu,sigma)

% y = zeros( size(x) );
% d = x - mu;
% ind  = d < 3 * sigma;
% y(ind) = exp(-0.5 * (d(ind)./sigma).^2) ./ (sqrt(2*pi) .* sigma); 

y = exp(-0.5 * ((x - mu)./sigma).^2) ./ (sqrt(2*pi) .* sigma); 