Формулы Бернулли, биномиальные коэффициенты, суммы рядов

Pe=10^(-2);
m=0:1:30;
Py=2*Pe*(1-Pe);
syms m r;
Pme=symsum((factorial(m)./(factorial(r).*factorial(m-r))).*(Py).^r.*(1-Py).^(m-r),r,(m+1)/2,m);
S=[2 2 2 1 1 1 1 1 1 1];
for k=1:10
    P(k)=0;
    for r=1:S(k)
        if mod(r,2)~=0
            P(k)=(factorial(S(k))/(factorial(r)*factorial(S(k)-r)))*(Pme)^r*(1-Pme)^(S(k)-r)+P(k);
        end;
    end;
end;
Pk=P;
t=3;
Pmult=1;
for k=1:10
    Pmult=Pmult*(1-2*Pk(k)*(1-Pk(k)));
end;
P=Pmult^(t-1);
plot(m,P)

clear, clc
m = 0:30;
Pe = [1e-2, 1e-3, 1e-4, 1e-5]; 
 
Pme = zeros(length(m),length(Pe)); % инициализируем
for j = 1:length(Pe)
  Py = 2*Pe(j)*(1-Pe(j));
  for i = 1:length(m)
     r = fix((m(i)+1)/2):m(i);
     Pem(i,j) = sum( (factorial(m(i))./(factorial(r).*factorial(m(i)-r))) .* Py.^r .* (1-Py).^(m(i)-r) );
   end
   str{j}=['Pe=10^',num2str(log10(Pe(j)))]; % дедаем подпись
end
plot(m, Pe)
grid on
legend(str)

clear, clc
m = 1:2:29;
Pe = [1e-2, 1e-3, 1e-4, 1e-5]; 
S=[2 2 2 1 1 1 1 1 1 1];
t=3;
 
Pme = zeros(length(m),length(Pe)); % инициализируем, создаем матрицу из нулей размером m  x Pe
Pk = zeros(length(S)); % инициализируем, создаем строку из нулей размером 1 x 10
P = zeros(length(m),length(Pe)); % инициализируем, создаем матрицу из нулей размером m  x Pe
Pmult = zeros(length(m),length(Pe)); % инициализируем, создаем матрицу из нулей размером m  x Pe
for j = 1:length(Pe)
  Py = 2*Pe(j)*(1-Pe(j)); % БПУ
  for i = 1:length(m)
     r = (m(i)+1)/2:m(i);
     Pme(i,j) = sum( (nchoosek(m(i),r)) .* Py.^r .* (1-Py).^(m(i)-r) );
     %Вероятность ошибочной оценки символа на выходе МЭ
     
     %Вероятность ошибки определения координат:
     for k=1:length(S)
         Pk(k)=0;
         for r=1:S(k)
             if mod(r,2)~=0
                 Pk(k)=(nchoosek(S(k),r))*((Pme(i,j))^r)*((1-Pme(i,j))^(S(k)-r))+Pk(k);
             end;
         end;
     end;
     
     %Вероятность правильной синхронизации:
     Pmult(i,j)=1;
     for k=1:length(S)
         Pmult(i,j)=Pmult(i,j)*(1-2*Pk(k)*(1-Pk(k)));
     end;
     P(i,j)=(Pmult(i,j))^(t-1); 
     
  end;
   str{j}=['Pe=10^',num2str(log10(Pe(j)))]; % делаем подпись
end
subplot(1,2,1); plot(m, Pme)
subplot(1,2,2); plot(m, P)
grid on
legend(str)

Error using nchoosek (line 24)
The second input has to be a non-negative integer.
 
Error in Ps2 (line 17)
     Pme(i,j) = sum( (nchoosek(m(i),r)) .* Py.^r .* (1-Py).^(m(i)-r) );