function pf=convpf(pf1, pf2) 
% convpf(pf1, pf2)  generates a probability function for the convolution of 
%                   the given probability functions 
% 
%    convpf is slow when one pf is discrete and the other is continuous 
% 
%    actually, it can be slow all the time 
% 
%    for a pair of discrete pmfs, the x values should be integers, or 
%    convpf will likely give you the wrong answer 
if ispf(pf1)==0 
   error('pf1 invalid') 
end 
if ispf(pf2)==0 
   error('pf2 invalid') 
end 
x1=xvals(pf1); 
x2=xvals(pf2); 
p1=pvals(pf1); 
p2=pvals(pf2); 
if isdiscrt(pf1)==0 & isdiscrt(pf2)==0 
   W1=x1(length(x1))-x1(1); 
   W2=x2(length(x2))-x2(1); 
   x2=(-x2(length(x2):-1:1)); 
   p2=p2(length(p2):-1:1); 
   pf2=[0 1; x2' p2']; 
   conp2=sample(pf2, -W2, W2, (W1+2*W2)/1000); 
   conx=[-W2:(W1+2*W2)/1000:W2]; 
   if ~(length(conp2)==length(conx)) 
      conx=conx(1:(length(conx)-1)); 
   end 
   slen=length(conp2); 
   conp1=sample(pf1, x1(1)-2*W2, x1(1)+W1+2*W2, (W1+2*W2)/1000); 
   newx=[x1(1)-W2:(W1+2*W2)/1000:x1(1)+W1+W2]; 
   if length(conp1)<length(conp2)+length(newx)-1 
      conp1=[conp1 0]; 
   end 
   newp=zeros(1,length(newx)); 
   for i=1:length(newx), 
      conp=(conp1(i:(i+length(conp2)-1))).*conp2; 
      P=trapz(conx, conp); 
      if P==[] 
         P=0; 
      end 
      newp(i)=P; 
   end 
   pf=[0 1; newx' newp']; 
elseif (isdiscrt(pf1)==0 & isdiscrt(pf2)==1) | (isdiscrt(pf1)==1 & isdiscrt(pf2)==0) 
   if isdiscrt(pf1)==1 
      conx1=x1; 
      conp1=p1; 
      conx2=x2; 
      conp2=p2; 
      W2=max(x1)-min(x1); 
      W1=max(x2)-x2(1); 
      cpf=pf2; 
      X=x2(1); 
   else 
      conx1=x2; 
      conp1=p2; 
      conx2=x1; 
      conp2=p1; 
      W2=max(x2)-min(x1); 
      W1=max(x1)-x1(1); 
      cpf=pf1; 
      X=x1(1); 
   end 
   newx=[X-W2:(W1+2*W2)/1000:X+W1+W2]; 
   newp=zeros(1,length(newx)); 
   for i=1:length(conx1), 
      cpf(2:(length(conx2)+1), 1)=(conx2+conx1(i))'; 
      P=sample(cpf, X-W2, X+W1+W2, (W1+2*W2)/1000); 
      if ~(length(P)==length(newp)) 
         P=[P 0]; 
      end 
      newp=newp+conp1(i).*P; 
   end 
   pf=[0 1; newx' newp']; 
else 
   TX=[x1 x2]; 
   X=[min(TX):max(TX)]; 
   P1=zeros(size(X)); 
   P2=zeros(size(X)); 
   index1=1; 
   index2=1; 
   for i=1:length(X), 
      if index1<=length(x1) 
         if x1(index1)==X(i) 
            P1(i)=p1(index1); 
            index1=index1+1; 
         end 
      end 
      if index2<=length(x2) 
         if x2(index2)==X(i) 
            P2(i)=p2(index2); 
            index2=index2+1; 
         end 
      end 
   end 
   newp=conv(P1, P2); 
   newx=[min(X):(min(X)+length(newp)-1)]; 
   if max(newp)>0 
      index1=1; 
      while newp(index1)==0 
         index1=index+1; 
      end 
      index2=length(newp); 
      while newp(index2)==0 
         index2=index2-1; 
      end 
      newp=newp(index1:index2); 
      newx=newx(index1:index2); 
   end 
   pf=[1 1; newx' newp']; 
end 
if inrange(ttlprob(pf), 0, 1e-4)==1 
   disp('Result is zero') 
elseif ispf(pf)==0 & isdiscrt(pf)==0 
   newp=newp./ttlprob(pf); 
   pf=[0 1; newx' newp']; 
end