function demo_EM(p,N)

format compact;
n=zeros(3,2);
p1=1/4;
p2=1/4 + p/4; P2=p1+p2;
p3=1/2-p/4;
Ep=0;
Nloop=10;

z=rand(1,N);
I1=find(z<=p1);
I2=find(z>p1 & z<=P2);
I3=find(z>P2);
n(1,1)=length(I1);
n(2,1)=length(I2);
n(3,1)=length(I3);

%n(1,1)=round(p1*N);     
%n(2,1)=round(p2*N);    
%n(3,1)=N-n(1,1)-n(2,1); 

m1=n(1,1)+n(2,1);
m2=n(3,1);

p_true = 2*(m1-m2)/(m1+m2);

figure(1);
subplot(1,3,2);
plot(0,0,'o');
axis([0 10 0 1]);
hold on;
subplot(1,3,3);
L=m1*log(0.5+0.25*Ep)+m2*log(0.5-0.25*Ep);
plot(0,L,'ro');
hold on;

for i=1:Nloop
   
   n(1,2)=m1/(2+Ep);
   n(2,2)=m1*((1+Ep)/(2+Ep));
   n(3,2)=m2;
   
   Ep=(2*n(2,2)-n(3,2))/(n(2,2)+n(3,2));
   
   L=m1*log(0.5+0.25*Ep)+m2*log(0.5-0.25*Ep);
   
   figure(1);
   subplot(1,3,1);
   bar(n);
   subplot(1,3,2);
   plot(i,Ep,'o');hold on; 
   plot([0 1:Nloop],p*ones(1,Nloop+1),'g');
   plot([0 1:Nloop],p_true*ones(1,Nloop+1),'y');   
   subplot(1,3,3);
   plot(i,L,'ro');
   drawnow;
   %pause
  
   
end

p_true
Ep

figure(1);
subplot(1,3,2);
hold off;
subplot(1,3,3);
hold off;