function pdf=dfwei(a, b, samples) 
% dfwei(a, b, samples)  Generates a Weibull pdf with parameters a, b. 
%    samples sets the size of the numerical representation, and is optional. 
%    The default sample size is 1000. 
is nargin==2 
   samples=1000; 
end 
if a<=0 
   error('invalid value for a') 
end 
if b<=0 
   error('invalid value for b') 
end 
V=sqrt(((1/a)^(2/b))*(gamma(1+2/b)-(gamma(1+1/b))^2)); 
X=5*V; 
while a*b*(X^(b-1))*exp(-a*(X^b)) > 1e-6 
   X=X+V; 
end 
x=[0:X/samples:X]; 
p=(a*b).*(x.^(b-1)).*exp(-a.*(x.^b)); 
pdf=[0 1; x' p']; 
if ispf(pdf)==0 
   p=p./ttlprob(pdf); 
   pdf=[0 1; x' p']; 
end