/* program to generate plots of density
   functions for noncentral f-distributions
   and overlay them on plots of the denstiy
   for the corresponding central f distribution */

proc iml;

    start fden;
    z = j(1000,3);

    do  i= 1 to 1000;
    x = i/200;

    a=5; b=20;              /*  Set the degrees of freedom  */

    d=3;                     /*  Enter the noncentrality parameter
                                 t(B)*t(X)*X*b/(2*(error variance))  */

                             /*  evaluate the central F density  */
    fab = (b**(b/2))*(a**(a/2))*(x**((a-2)/2));
    fab = fab/(((b+a*x)**((b+a)/2))*gamma(b/2)*gamma(a/2)/gamma((b+a)/2));

                             /*  evaluate noncentral f density  */

    sum1 = 1;
    term=1;
    p1 = ((d*a*x)/(b+a*x));
    p2 = 1;
    nt = 100;
    do j = 1 to nt;
      term = term*p1*(a+b+2*(j-1))/((a+2*(j-1))*j);
      sum1 = sum1 + term;
      end;
    fabnc = fab*exp(-d)*sum1;
    z[i,] = x||fab||fabnc;
    end;

    create set1 from z;
       append from z;

    finish;


  run fden;


filename graffile pipe 'lpr -Dpostscript';
goptions gsfmode=replace gsfname=graffile;
goptions cback=white ctext=black targetdevice=ps300
         rotate=landscape;


  axis1 label = (h=3 r=0 a=90 f=swiss 'density')
        length = 5 in
        color = black width=8.0 style=1
        order = 0 to 0.8 by .2;

  axis2 label = (h=2.5 f=swiss 'f')
        length = 5.5 in
        color=black width=8.0 style=1
        order = 0 to 5 by 1;

  symbol1 v=none i=spline l=1 h=2 c=black width=8;
  symbol2 v=none i=spline l=3 h=2 c=black width=8;
  symbol3 v=none i=spline l=9 h=2 c=black width=8;
  symbol4 v=none i=spline l=16 h=2 c=black width=8;

proc gplot data=set1;
  plot (col2 col3)*col1 / overlay vaxis=axis1 haxis=axis2;
  title h=3 f=swiss  'Densities for Central and Noncentral';
  title2 h=3 f=swiss  'F Distributions with (5,20) df';
  title3 h=3 f=swiss '(noncentrality parameter = 3.0)';
  run;