function [f,mu,sigma] = pdf(typ, x, p1, p2)
   // Probability density function of several distributions
   //
   // Usage:
   //       [f, mu, sigma] = pdf(typ, x, p1, p2)
   // where
   //        typ:    Type of distribution:
   //                "nor":  normal
   //                "t":    t (Student)
   //                "chi2": chi square
   //                "F":    F (Fisher-Snedecor)
   //        x:      scalar, vector or array
   //        p1, p2: Parameters of the distribution
   //                (see table below)
   //
   // Parameter table (default between parenthesis)
   // Distribution    p1        p2      
   // ------------------------------------             
   // normal          mu (0)    sigma (1)  
   // t               nu (10) 
   // chi square      nu (20)
   // F               m  (3)    n     (10) 
   // ------------------------------------
   //
   // NOTES: 
   // 1) While in the normal distribution it is
   //    frequent to use a mean (mu) and standard 
   //    deviation (sigma) different from 0 and 1,
   //    for the other cases the data are always 
   //    transformed to the standard version.
   // 2) If no argument is provided the normal 
   //    distribution N(0,1) is computed on an 
   //    interval from -3 to 3
   //
   // -----------------------
   // Author: Federico Miyara
   // Date:   2025-05-04
   
   // -----------------------
   // Check arguments
   // -----------------------
   rhs = argn(2);
   if rhs==0
      typ = "nor";
      x = -3:0.1:3;
      p1 = 0;
      p2 = 1;
   else
      typs = ["nor" "t" "chi2" "F"];
      // Detect unsupported distributions
      if ~or(typ==typs)
         msg = "%s: Unsupported distribution"
         error(msprintf(msg, "pdf"))
      end       
   end
   if rhs==1
      if typ=="nor" | typ=="t"
         x = -3:0.1:3;
      elseif typ=="chi2"
         x = 0:0.1:20
      else
         x = 0:0.1:5
      end  
   end
   if rhs<3
      if typ=="nor"
         p1 = 0;
         p2 = 1;
      elseif typ=="t" | typ=="chi2" 
         p1 = 10;
      else
         p1 = 3
         p2 = 10
      end
   end
   if rhs<4
      if typ=="nor"
         p2 = 1;
      elseif typ=="F"
         p2 = 10
      else
      end
   end
   
   // -----------------------
   // Compute PDF
   // -----------------------
   select typ
   case "nor"
      // Mean and standard deviation
      mu = p1;
      sigma = p2;
      // Probability density (normal)
      f = 1/sqrt(2*%pi)/p2 * exp(-((x-p1)/p2).^2);
   case "t"
      // Degrees of freedom
      nu = p1;
      // Mean and standard deviation
      mu = 0;
      sigma = sqrt(nu/(nu-2));      
      // Probability density (t)
      f = gamma((nu+1)/2)/gamma(nu/2)/sqrt(%pi*nu) * ..
          (1 + x.^2/nu).^(-(nu+1)/2);
   case "chi2"
      // Degrees of freedom
      nu = p1;
      // Mean and standard deviation
      mu = nu;
      sigma = sqrt(2*nu);
      // Probability density (chi square)
      f =  1/2^(nu/2)/gamma(nu/2) * x.^(nu/2-1).*exp(-x/2);
      i0 = find(x<0)
      f(i0) = 0;
   case "F"
      // Degrees of freedom of numerator(m) and
      // denominator (n)
      m = p1;
      n = p2;
      // Mean and standard deviation
      mu = n/(n-2);
      sigma = sqrt(2*n^2*(m+n-2) / (m*(n-2)^2*(n-4)))
      // Probability density (F)
      f = 1./x/beta(m/2,n/2) .* ..
          sqrt((m*x).^m * n^n ./ (m*x + n).^(m+n));
      // For those indexes where x = 0 replace NaNs 
      // by the limit 
      i0 = find(x<=0)
      if x(1)==0
         if m==2 
            f(i0) = sqrt((m/n)^m)/beta(m/2,n/2);
         else
            f(i0) = 0;
         end
      end
   end
                      
endfunction
   
