#
#  以下程式測試 Antithetic 估計　P(Cauchy > 2) 的機率
#  --> 第一種估計方法(Hit-or-miss)、第二種估計方法
#
# (1)
#a1=sum(x>2)/length(x) 
# a=(x>2)*1
#  b1=var(a)
#
# (2) 
#a2=sum(abs(x)>2)/length(x)
# aa=(x>2 | x< -2)*1
#  b2=var(aa)/4
#
t1=NULL
for (i in 1:1000) {
 a=runif(10000)
  b=1-a
   x1=qcauchy(a)
    x2=qcauchy(b)
   y1=sum(x1>2)/10000
  y2=sum(x2>2)/10000
 t1=cbind(t1,c(y1,(y1+y2)/2))
}
#
apply(t1,1,mean)
 a=apply(t1,1,var)
  a[1]/a[2]
#
#  前兩種方法的差異！
#
####################
####################
#
# 第三種估計方法
#  
#x1=runif(10000,0,2)
# y1=2/(pi*(1+x1^2))
#  z1=0.5-y1
#    a3=mean(z1)
#     b3=var(z1)
t1=NULL
 for (i in 1:1000) {
  x1=runif(10000,0,2)
   x2=2-x1
    y1=2/(pi*(1+x1^2))
     y2=2/(pi*(1+x2^2))
      z1=0.5-y1
      z2=0.5-(y1+y2)/2
     a1=mean(z1)
    a2=mean(z2)
   t1=cbind(t1,c(a1,a2))
  }
apply(t1,1,mean)
apply(t1,2,var)
#
# (3)-Control
x1=runif(10000,0,2)
 a=0.5-y1-0.15*(2*x1^2-8/3)+0.025*(2*x1^4-32/5)
  mean(a)
 var(a)
yy1=1/(pi*(1+x1^2))+1/(pi*(1+(2-x1)^2))
  zz1=0.5-yy1
    mean(zz1)
     var(zz1)

#
####################
####################
#
# 第四種估計方法
#  
# (4) & (4)-control, (4)-antithetic 
x2=runif(10000,0,0.5)
  y2=1/(2*pi*(1+x2^2))
    a4=mean(y2)
      b4=var(y2)
a=y2+0.312*(0.5*x2^2-1/24)-0.233*(0.5*x2^4-1/160)
 mean(a) 
  var(a) 
yy2=0.5/(2*pi*(1+x2^2))+0.5/(2*pi*(1+(0.5-x2)^2))
  mean(yy2)
   var(yy2)