program='F_illustration'
ps=.1
psdirect,program,ps,xsize=13,ysize=5,/vfont
nwin

;make pdf of F-distribution
;F_PDF(V, Dfn, Dfd) returns cumulative 

dfn=20 & dfd=20
x=findgen(1000)*.005   & apu=x*0. & pdf=apu
for i=0,n_elements(x)-1 do begin
    apu(i)=f_pdf(x(i),dfn,dfd)
    IF(I NE 0) then pdf(i)=(apu(i)-apu(i-1))/(x(i)-x(i-1))
endfor

plot,x,pdf,xtitle='test statistic F',ytitle='probability',$
  title='Sampling distribution dfn=20,dfd=20,alpha=0.05',$
  charsize=0.8

;F_CVF(P, Dfn, Dfd) returns value Fcut,   prob(F>Fcut)=P
alpha=0.05
xupper=F_cvf(alpha/2,dfn,dfd)
xlower=F_cvf(1-alpha/2,dfn,dfd)             ;these two forms
xlower=1./F_cvf(alpha/2,dfd,dfn)            ;give the same result

indl=where(x lt xlower)
xapu=[0,x(indl),reverse(x(indl))]
yapu=[0,pdf(indl),pdf(indl)*0]
polyfill,xapu,yapu
xyouts,.4,.1,'!7a!3/2'

indl=where(x gt xupper)
xapu=[x(indl(0)),x(indl),reverse(x(indl))]
yapu=[0,pdf(indl),pdf(indl)*0]
polyfill,xapu,yapu

xyouts,3.5,.1,'!7a!3/2'
psdirect,program,ps,/stop


;-------------------------------------------------------------
psdirect,program+'_a',ps,xsize=13,ysize=6
!p.multi=[0,2,1]

X = [257, 208, 296, 324, 240, 246, 267, 311, 324, 323, 263, $
305, 270, 260, 251, 275, 288, 242, 304, 267]
Y = [201,    56, 185, 221, 165, 161, 182, 239, 278, 243, 197, $
271, 214, 216, 175, 192, 208, 150, 281, 196]
res=fv_test(x,y)

histo_f,x,0,400,25,xt1,yt1,/nos
histo_f,y,0,400,25,xt2,yt2,/nos

nwin
plot,xt1,yt1,psym=10,yr=[0,10]
oplot,xt2+2,yt2,psym=10,lines=2,thick=3
xyouts,10,9,'F='+string(res(0))
xyouts,10,8.5,'prob(1-sided)='+string(res(1))

;remove one 'bad' Y point
Y=y(where(Y ne 56))
res=fv_test(x,y)

histo_f,x,0,400,25,xt1,yt1,/nos
histo_f,y,0,400,25,xt2,yt2,/nos

nwin
plot,xt1,yt1,psym=10,yr=[0,10]
oplot,xt2+2,yt2,psym=10,lines=2,thick=3
xyouts,10,9,'F='+string(res(0))
xyouts,10,8.5,'prob(1-sided)='+string(res(1))

!p.multi=0

psdirect,program+'_a',ps,/stop
end