Моделирование дискретной случайной величины

N=100;
p=(3^(i))/(4^(i+1));
k=[0:30];
for i = 1:N;
 z=rand(1);
 S=z;
 i=1;
 while S>0
 S=S-p(i);
 i=i+1;
 end
 x(j)=k(i-1)
end

clear, clc
 
k = 0:30;
p = (3.^k)./(4.^(k+1));
F = cumsum(p);
 
% теоретические характеристики
Mt = sum(p.*k);
Dt = sum(p.*(k-Mt).^2);
disp(['Теоретическое мат.ожидание: ',num2str(Mt)])
disp(['Теоретическая дисперсия: ',num2str(Dt)])
 
for N = [100 1000 10000]
    disp(' ')
    disp(['Длинна выборки, N: ',num2str(N)])
    x = zeros(1,N);
    for i = 1:N;
        S = rand;
        j = 1;
        while S>0 & j<=length(p)
            S = S-p(j);
            j = j+1;
        end
        x(i) = k(j-1); % находим 
    end
    Mx = mean(x); % мат.ожидание
    Dx = var(x); % дисперсия
    disp(['Экспериментальное мат.ожидание: ',num2str(Mx)])
    disp(['Экспериментальная дисперсия: ',num2str(Dx)])
    
    NN = zeros(size(k));
    for i = 1:length(k)
        NN(i) = sum( x==k(i) )/length(x); % считаем кол-во k-х элементов
    end
    
    figure % рисуем все
    stem(k,p,'sr'), hold on
    stem(k,NN,'b')
    title (['Плотность вероятности для N=',num2str(N)])
    legend('Теоретическая','Экспериментальная')
    xlabel('k')
    ylabel('p')
end