#This list of commands generates 100 samples of size 100 each
# and then creates 4 alternative estimates of 1-pcauchy(2).  Theoretical
# and empirical estimates of standard errors are computed
x=matrix(rcauchy(10000),ncol=100)
u2=matrix(runif(10000,0,2),ncol=100)
up5=matrix(runif(10000,0,.5),ncol=100)
#Construct estimates from each trial; we could also estimate standard
#errors for each trial using MC integration
theta1=apply(x,1,function(x){mean(x>2)})
theta2=apply(x,1,function(x){.5*mean(abs(x)>2)})
theta3=apply(u2,1,function(x){.5-mean(2/(pi*(1+x^2)))})
theta4=apply(up5,1,function(x){mean(1/(2*pi*(1+x^2)))})
vref=1-pcauchy(2)
par(mfcol=c(2,2))
hist(theta1,main="Histogram of Theta 1", xlim=c(.1,.2))
abline(v=vref,col="red")
hist(theta2,main="Histogram of Theta 2", xlim=c(.1,.2))
abline(v=vref,col="red")
hist(theta3,main="Histogram of Theta 3", xlim=c(.1,.2))
abline(v=vref,col="red")
hist(theta4,main="Histogram of Theta 4", xlim=c(.1,.2))
abline(v=vref,col="red")
par(mfcol=c(1,1))
paste("Empirical variance for Theta 1 is", var(theta1))
paste("Theoretical variance for Theta 1 is",.126/100)
paste("Empirical variance for Theta 2 is", var(theta2))
paste("Theoretical variance for Theta 2 is",.054/100)
paste("Empirical variance for Theta 3 is", var(theta3))
paste("Theoretical variance for Theta 3 is",.0285/100)
paste("Empirical variance for Theta 4 is", var(theta4))
paste("Theoretical variance for Theta 4 is",.0000955236/100)