Моделирование случайной величины с гамма распределением

clear; clc;
 
% параметры
N = 10000; % кол-во полученных чисел с гамма распределением
k = 10; % параметр распределения (кол-во экспоненциальных чисел в сумме)
for i = 1:N
    k1 = randn(1,k); % нормальная величина (mu=0, sigma=1)
    k2 = randn(1,k); % тоже самое
    z = k1.^2 + k2.^2; % экспоненциальная величина (lambda=1+1=2)
    X(i) = sum(z); % сумма экспоненциальных
end
 
% проверяем шо мы тут натворили
% строим гистограмму
[n xn] = hist(X,50); % 50 диапазонов
y = n/(length(X)*(xn(2)-xn(1))); % считаем плотность по гистограмме
t = linspace(min(X),max(X),1000);
plot(t,pdf('gam',t,k,2), xn,y,'r')
legend('Теорет. плотность','Гистограмма полученных чисел')

clear; clc;
 
% параметры
N = 10000; % кол-во полученных чисел с гамма распределением
k = 10; % параметр распределения (кол-во экспоненциальных чисел в сумме)
L = 2; % параметр экспоненциального распределения
for i = 1:N
    r = rand(1,k); % равномерная величина (a=0, b=1)
    z = -L*log(r); % экспоненциальная величина (lambda=L)
    X(i) = sum(z); % сумма экспоненциальных
end
 
 
% или вообще в одну строку:
% X = sum( -L*log( rand(k,N) ) );
 
% проверяем шо мы тут натворили
% строим гистограмму
[n xn] = hist(X,50); % 50 диапазонов
y = n/(length(X)*(xn(2)-xn(1))); % считаем плотность по гистограмме
t = linspace(min(X),max(X),1000);
plot(t,pdf('gam',t,k,L), xn,y,'r')
legend('Теорет. плотность','Гистограмма полученных чисел')