% Gamma distribution
clear all;

% Parameters
alpha = 2;
beta = 80;

% x and y values
x = 0:0.1:1000;
f = (1./(gamma(alpha)*beta^alpha)).*(x.^(alpha-1)).*exp(-x/beta);

% Random draws from a Gamma
y = random('Gamma',alpha,beta,100000,1);

% Graphs
hold on;
plot(x,f)
histogram(y,50,'Normalization','pdf');