N=200;
Nb=15; %ilosc binow
x=sum(rand(200,4),2);
range=[min(x) max(x)+.001];
dx=(range(2)-range(1))/Nb;
bins=range(1):dx:range(2);
centers=range(1)+dx/2:dx:range(2);
subplot 211
h=histc(x,bins);
hb=bar(bins,h,'histc');
set(gca,'xlim', range, 'xtick', [range(1):dx:range(2)], ...
    'xticklabel',[],'xgrid','on','ytick',[])
set(hb(1),'facecolor','y');
hold on
mu=mean(x);
sigma=sqrt(var(x));
for i=1:length(centers)
    Pi(i)=normcdf(bins(i+1),mu,sigma)-normcdf(bins(i),mu,sigma );
end
hc=bar(centers,N*Pi,.2)
set(hc(1),'facecolor',[.5 .5 .5]    );

subplot 212
for i=1:Nb
    chi2(i)=(h(i)-N.*Pi(i)).^2 / (N*Pi(i));
end

s_chi2=sum(chi2);

t=0:0.1:2*Nb;
plot(t, chi2pdf(t,Nb-3));
set(gca, 'xtick', s_chi2, 'xgrid', 'on');
p=chi2cdf(s_chi2, Nb-3);
if p>.5 p=1-p; end
text(1, 0.094, sprintf('statystyka \\chi^2 dla %d stopni swobody', Nb-3))
text(s_chi2+.5, 0.01, sprintf('p=%.2f', p))
