% FILE  hlp051.m  begins.
clear; close all;
N=1000;
u=rand(1,N); v=rand(1,N); 
z1=sqrt(-2*log(u)).*cos(2*pi*v); 
z2=sqrt(-2*log(u)).*sin(2*pi*v);
data=z1+3;
a=histogram(data,40); % 40 is number of bins
hlp=a.BinEdges;
y=a.Values;    % sum(y)=1000; length(hlp) =1+length(y)
y=y/sum(y);    %sum(y) = 1
x=hlp(2:end);  % length(x) =length(y)
fmodel=@(t,lam) (lam(1)*exp(-lam(2)* (t-lam(3)).^2));
fobj=@(lam) norm(fmodel(x,lam)-y);
lam0=[-1,1,2];
v1=fobj(lam0);  % Value before fitting
lam=myparf(x, y,fmodel,lam0)
v2=fobj(lam);   % Value after fitting
figure
fprintf('fobj  before/after fitting: %f   %f\n',v1,v2);
yfit=fmodel(x,lam);
axes('FontSize',[15],'FontWeight','bold'); hold on;
bar(x,y)
hold on
title(['fobj: before =' num2str(v1,3) ', after = ' ...
       num2str(v2,3)])
plot(x,yfit,'k-',x,y,'k.','MarkerSize',20,'LineWidth',2);  
disp(['fmodel >data: ' num2str(sum(yfit>y)) '; fmodel <= data: ' num2str(sum(yfit<=y))])
grid on;
ax=axis; x1=ax(1)+0.1*(ax(2)-ax(1)); y1=ax(3)+0.9*(ax(4)-ax(3));
text(x1,y1, ['lam= ' mat2str(lam,3)],'FontWeight','bold','FontSize',[15]);
xlabel(func2str(fmodel))
%lam(2) =1/(2*sig^2)
sighat=1/sqrt(2*lam(2));
fprintf('muhat, sighat = %f   %f\n', lam(3),sighat)



function w=myparf(xdata1,ydata1,mymodel, laminit)
% E.g. fmodel=@(x,lam)exp(-lam(1)* x)+lam(2)*exp(-lam(3)*x);
fobj=@(lam) norm(mymodel(xdata1,lam)-ydata1);
[lam fval eflag] = fminsearch(fobj,laminit);
w=[lam ];
end
% FILE  parfhisto.m  ends.
