home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / STBETA.C < prev    next >
Encoding:
Text File  |  1984-08-15  |  1.4 KB  |  68 lines
  1.    stbeta (n,m,x,h,s)
  2.  
  3.       /* this function computes the density function df = x**(m-1)*       */
  4.       /* (1-x)**(n-1)/b(m,n) and cumlative function of beta distribution. */
  5.  
  6.       int m,n;
  7.       float x,*h,*s;
  8.  
  9.    {
  10.  
  11.       int i,ii;
  12.       float fm,fn,fi,fs,h1,h2,hh,s1,s2,fs1,fs2;
  13.       extern double pow();
  14.  
  15.       fn = n;
  16.       fm = m;
  17.       *s = pow(x,fm)/fm;
  18.       fs = 1.;
  19.  
  20.       for(i = 1; i <= 1000; i++)
  21.       {
  22.        fi = i;
  23.        fs = fs * (fn - fi);
  24.        fs = fs / fi;
  25.        ii=i/2;
  26.        if((ii*2) == i) hh = 1.;
  27.         else hh  = -1.;
  28.        hh = hh * pow(x,(fm + i))/(fm + fi);
  29.        hh = hh * fs;
  30.        *s = *s + hh;
  31.        if((fn - fi) == 0.) break;
  32.       }
  33.       s1 = 1./fm;
  34.       fs1 = 1.;
  35.  
  36.       for(i = 1; i <= 1000; i++)
  37.       {
  38.        fi = i;
  39.        fs1 = fs1 * (fn - fi);
  40.        fs1 = fs1 / fi;
  41.        ii=i/2;
  42.        if((ii*2) == i) h1 = 1.;
  43.         else h1  = -1.;
  44.        h1 = h1/(fm + fi);
  45.        h1 = h1 * fs1;
  46.        s1 = s1 + h1;
  47.        if ((fn - fi) == 0.) break;
  48.       }
  49.       *s = *s / s1;
  50.       s2 = pow(x,(fm - 1.));
  51.       fs2 = 1;
  52.  
  53.       for(i = 1; i <= 1000; i++)
  54.       {
  55.        fi = i;
  56.        fs2 = fs2 * (fn - fi);
  57.        fs2 = fs2 / fi;
  58.        ii=i/2;
  59.        if((ii*2) == i) h2 = 1.;
  60.         else h2  = -1.;
  61.        h2 = h2 * pow(x,(fm + fi - 1.));
  62.        h2 = h2 * fs2;
  63.        s2 = s2 + h2;
  64.        if ((fn - fi) == 0.) break;
  65.       }
  66.       *h = s2 / s1;
  67.    }
  68.