function f = GEV_pdf ( x, param )
% f = GEV_pdf ( x, [ m , s, k ] )
% compute the PDF of the generalized extreme value distribution 
% with parameters m, s, and k
%  f = (1/s)*(1+k*z)^(-1-1/k) * exp(-(1+k*z).^(-1/k) ); where  z = ( x - m ) / s;
% usually -0.5 < k < +0.5

 m = param(1);	% shift
 s = param(2);	% scale
 k = param(3);	% shape

 z = ( x - m ) / s;

 f = real ( (1/s)*exp(-(1+k*z).^(-1/k)) * (1+k*z).^(-1-1/k) );

 f(find(1 + k*z < 0)) = eps;

% GEV_pdf ---------------------------------------------------- 
% HP Gavin 31 Feb 2012 2019-11-20



