p=0.187966*10^(-2);%modelul gilbert-elliot
q=0.080617;
pg=0.563331*10^(-3);
pb=0.417201;
m=12;
N=15;
a=logspace(0,2,13);
%a=1:13;
%crearea matricelor initiale
g=zeros(m+1,N);
b=zeros(m+1,N);
g(1,:)=(1-pg)*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
g(2,:)=pg*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
b(1,:)=(1-pb)*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
b(2,:)=pb*[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
%calculul elementelor prin recurenta
for k=2:N
   g(1,k)=g(1,k-1)*(1-p)*(1-pg)+b(1,k-1)*p*(1-pg);
   b(1,k)=b(1,k-1)*(1-q)*(1-pb)+g(1,k-1)*q*(1-pb);
end

for i=2:m,
   for j=2:N,
     g(i,j)=g(i,j-1)*(1-p)*(1-pg)+g(i-1,j-1)*(1-p)*pg+b(i,j-1)*p*(1-pg)+b(i-1,j-1)*p*pg;
     b(i,j)=b(i,j-1)*(1-q)*(1-pb)+b(i-1,j-1)*(1-q)*pb+g(i,j-1)*q*(1-pb)+g(i-1,j-1)*q*pb;
   end
end
%calculul ponderii de eroare
PgemN=(q/(p+q))*g(:,N).'+(p/(p+q))*b(:,N).';
loglog(a,PgemN,'*r')
t0=clock;
etime(t0,clock)