'plots for goodness-of-fit tests
'for version 4
'12/18/2000 h.


'--------------------------------------------------------------------------------------
'subroutines
'--------------------------------------------------------------------------------------

'evaluate weibull density at xt and store in fx
subroutine eval_weibull(series xt, series fx, scalar m, scalar s, scalar a)
	series zt = (xt - m)/s
	fx = (a/s) * zt^(a-1) * exp( -(zt^a) )
endsub

'evaluate Gaussian density at xt and store in fx
subroutine eval_norm(series xt, series fx, scalar m, scalar s)
	series zt = (xt - m)/s
	!pi = @acos(-1)
	fx = (1/@sqrt(2*!pi)/s) * exp( -(0.5*zt^2) )
endsub

'--------------------------------------------------------------------------------------
'main program
'--------------------------------------------------------------------------------------

'overlay theoretical and empirical distribution
workfile edfplot u 1 500

'set random number generator
rndseed(type=mt) 1234567
'simulate data
series x = @rchisq(5)

'kernel density estimate
!ngrid = 150			'number of points to evaluate
do x.kdensity(!ngrid, b=2.3, o=kdxy)

'make sure workfile is large enough
if (@obsrange < !ngrid) then
	expand 1 !ngrid
endif

'convert stored estimates into series
smpl 1 !ngrid
series xt
series ft_kdensity
group group1 xt ft_kdensity
mtos(kdxy, group1)

'evaluate theoretical distribution at xt
series ft_dist
call eval_norm(xt, ft_dist, @mean(x), @stdev(x))

graph g1.xyline xt ft_dist ft_kdensity
g1.option linepat	' need to set linepat
g1.setelem(1) legend()
g1.setelem(2) legend(normal)
g1.setelem(3) legend(kernel density)
g1.setelem(1) lcolor(blue) lpat(dash1)
g1.setelem(2) lcolor(red) lpat(solid) 
g1.legend columns(1) -inbox position(2.3,0)

'qq-plot
freeze(g2) x.qqplot(n)
g2.setelem(1) legend()
g2.addtext(0.1,0.1) QQ-plot 

'combine two graphs
graph graph1.merge g2 g1
graph1.align(2,0,0)
show graph1

'display goodness-of-fit tests
show x.edftest(type=normal,showopts)

