## Monte Carlo estimate of area of a quarter circle of unit radius
## Find sample size in K1.
## Find number of repetitions in K2.
## Use c1-c6 and K96-K99 as scratch.
## Put sample proportions in c3.
## Put rescaled sample proportions in c4.


gmacro
monte

name k1 'sample size' k2 'number of repetitions'
print 'sample size' 'number of repetitions' 

erase c1-c6  # kill current contents of columns

do k99 = 1: k2
Random k1 c1 c2;
	uniform 0 1.
let c3(k99)= mean(c1*c1 + c2*c2 <= 1)
enddo


## the values in c3 should be observations on
## Binomial(sample size = k1, prob = pi/4)/k1

let k96 = 3.14159/4 # true of area, pi/4,  of quarter circle
let k97 = sqrt(k96*(1-k96)/k1)
let k98 = mean(c3)
let k99 = stdev(c3)
name k96 'true p' k97 'true std.dev'&
	k98 'sample proportion' k99 'sample std. dev.'

print 'sample size' 'number of repetitions' &
	'true p'  'true std.dev'&
	'sample proportion''sample std. dev.'

## subtract theoretically true proportion,
## divide by true std. dev
let c4 = (c3- k96)/k97
name c4 'rescaled sample proportions'

## calculate points for standard normal density
Set c5
	 1( -3 : 3 / .1 )1
End.
PDF c5 c6;
	Normal 0 1.

histogram c4;
	density;
	Line c5 c6;  # draw standard normal density
	overlay;
	nintervals 20.

endmacro