function normdenhist(data, mu, sigma, numbins)

if nargin<4; numbins = 10; end
% plot histogram and a normpdf
  n = length(data);
  binwidth = range(data)/numbins;
  edg= min(data):binwidth:max(data);
  [count, bin] = histc(data, edg);
  h = bar(edg, count./(n*binwidth), 'histc');
  %set(h, 'facecolor', [0.8 0.8 1]); % change the color of the bins
  %set(h, 'edgecolor', [0.8 0.8 1]);
  
  dx= (max(data)-min(data))/(n*100);
  
  tmp = 1;
  if (min(data)<0), a = - tmp*min(data);else a = tmp*min(data); end
  if (max(data)<0), b = - tmp*max(data);else b = tmp*max(data); end

  x= (min(data) - a):dx:(max(data) + b);
  density = normpdf(x,mu,sigma);
  
  
   coeff_norm = (max(count./(n*binwidth)))/max(density);
%    hold on,plot(x,coeff_norm*density,'r');
     hold on,plot(x, density,'r','linewidth',2)
     xlabel('x'); ylabel('f(x)')