/* 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;