Mathcad to Matlab - случайные процессы и NPS

Question

Я пытаюсь преобразовать простой файл Mathcad в скрипт MATLAB, но результаты, которые я получаю, просто странные. Есть ли какая-то разница в функциях rand? Я забыл что-то важное, что отличается от этих двух сред?

Вот как выглядит файл Mathcad:

И это код, который я сейчас использую:

clear all
clc

N=100;
data=poissrnd(1000,N,N);
data2=zeros(N);
for l=1:N
    for k=1:N
       drozp=zeros(1,3);
       eta_fot=data(l,k);
       for m=1:eta_fot
          ran=rand(1);
          if ran<0.25
              drozp(1)=drozp(1)+1;
          elseif ran>=0.75
              drozp(3)=drozp(3)+1;
          else
              drozp(2)=drozp(2)+1;
          end
       end
       if (k>1) && (k<N)
          data2(l,k)=data2(l,k)+drozp(2);
          data2(l,k-1)=data2(l,k-1)+drozp(1);
          data2(l,k+1)=data2(l,k+1)+drozp(3);
       elseif k==1
          data2(l,k)=data2(l,k)+drozp(2);
          data2(l,N)=data2(l,N)+drozp(1);
          data2(l,k+1)=data2(l,k+1)+drozp(3);
       else
          data2(l,k)=data2(l,k)+drozp(2);
          data2(l,1)=data2(l,1)+drozp(3);
          data2(l,k-1)=data2(l,k-1)+drozp(1);
       end

    end
end

data2 = data2 - mean(data2(:));
data_w = (1/sqrt(N))*fftshift(fft(data2, [], 2));

NPS1=(abs(data_w)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;

plot(NPS)

Теперь я знаю, что он может работать быстрее, изначально этот блок:

for m=1:eta_fot
          ran=rand(1);
          if ran<0.25
              drozp(1)=drozp(1)+1;
          elseif ran>=0.75
              drozp(3)=drozp(3)+1;
          else
              drozp(2)=drozp(2)+1;
          end
       end

выглядело так:

[~,R] = histc(rand(1,eta_fot),cumsum([0;w(:)./sum(w)]));
R=a(R);
drozp=histc(R,a)

где:

a=1:3;
w=[0.25,0.5,0.25];

Однако я попытался приблизиться к файлу Mathcad. С другой стороны, результаты, которые я получаю (независимо от того, каким образом я это делаю) выглядят просто плохо. Смотри картинку ниже.

Единственный способ получить результаты, аналогичные результатам из файла Mathcad, - это изменить весь цикл rand на "жесткий" drozp(1)=0.25*eta_fot, drozp(2)=0.5*eta_fot а также drozp(3)=0.25*eta_fot но тогда это больше не имеет ничего общего с вероятностью. Затем снова цикл rand выглядит (по крайней мере, функционально) как тот, который используется в Mathcad, поэтому я понятия не имею, что не так.

Я что-то упускаю? Что-то я не смог понять? Я в недоумении.