Моделирование случайной величены с бета-распределением методом Неймана

clear, clc
 
% параметры распределения:
alpha = 2;
betta = 3;
% интервал по х:
a = 0;
b = 1;
% формула плотности вероятности
w = @(x) pdf('beta',x,alpha,betta); % используем встроенную формулу
 
x = linspace(a,b,500); % разбиваем диапазон на 500 точек
W = w(x); % считаем плотность
mW = max(W); % находим макс. значение
 
N = 1000; % длинна выборки
X = zeros(1,N); % пустой массив под числа
n = 1; % счетчик
while n<=N % пока не наберем нужно кол-во чисел
    % генерируем случайные координаты:
   ay = mW*rand;
   ax = a + (b-a)*rand;
   if ay <= w(ax) % если точка попала под график
       X(n) = ax; % сохраняем координату
       n = n+1; % увеличиваем счетчик
   end
end
% проверяем, шо мы там на генерировали:
[yn2 xn2] = hist(X,50); % получаем гистограмму
yn2 = yn2/(length(X)*(xn2(2)-xn2(1))); % считаем плотность по гистограмме
 
% рисуем:
plot (x,W,'r',... %теор
    xn2,yn2,'b') % эксперимент
grid on
title( ['Бетта-распределение с параметрами \alpha=',num2str(alpha),' , \beta=',num2str(betta)] )
legend('Теоретическая','Экспериментальная')