% Gamma PDF Distribution Analysis
% A.A. Trani (October 10, 1997)

% Functions call: gammads.m , gammadf.m

global A B

% Define the parameters of the distribution A and B

A = 1;
B = 2;

% Define upper and lower bounds to get a nice plot

npoints = 100;
parameter = 4;
means = A/B	;						% expected value of the distribution
sigma = sqrt(A/B^2)	;				% standard deviation of distribution
low = 0;
high = means + parameter* sigma;
interval = (high - low) / npoints;

% Define random variable x

x=low:interval:high;


% Define the function of the random variable x (PDF function)

fx = B^A * x.^(A-1) .* exp(-B .* x) / gamma(A);

% Plot the random variable x versus the PDF function

plot(x,fx)
xlabel('Random variable - x')
ylabel('f(x)')
grid

% Compute area under the curve f(x) from -infinity to infinity
% Change the values of the last two arguments if you 
% want another range of values

Pn=quad8('gammad',low,high);

  disp(' ')
  disp([blanks(5),'Area under the Gamma Distribution : '])
  disp(' ')
  disp([blanks(5),'     P(x>0)=',num2str(Pn)])
  disp(' ')
  
% Compute the first and second moments of the distribution
%
% The first moment is the mean of the distribution
% the second moment is the variance

mean_gamma = quad8('gammadf',0,100);
var_gamma = quad8('gammads',0,100);

  



