uniforms <- runif(500)
tosses   <- as.numeric(I(uniforms > 0.5))
relfreq  <- cumsum(tosses)/(1:500)
plot(relfreq)
abline(0.5,0)
title(main="Illustration of the Weak Law of Large Numbers")

par(mfrow=c(2,2))
for(rep in 1:4)
{
    uniforms <- runif(500)
    tosses   <- as.numeric(uniforms > 0.5)
    relfreq  <- cumsum(tosses)/(1:500)
    plot(relfreq)
    abline(0.5,0)
}

par(mfrow=c(2,2))
for(rep in  1:4)
{
    cauchys <- rcauchy(500)
    xbar    <- cumsum(cauchys)/(1:500)
    plot(xbar)
    abline(0,0)
}

poisson.clt    <- function(k,n, parameter)
{
  samples.mean   <- rep(NA, k)
  for (i in 1:k)
  {
    samples.mean[i]   <- mean(rpois(n, lambda=parameter))
  }
  return(samples.mean)
}

test.pois <- poisson.clt(200,100,1)
summary(test.pois)
var(test.pois)
par(mfrow=c(1,3))
hist(test.pois)
boxplot(test.pois)
qqnorm(test.pois)
qqline(test.pois)

p<-c(0.01,0.1,0.3,0.5)
n<-c(10,100,1000,10000)

par (mfrow=c(4,4))
for (i in 1:4)
{
    for (j in 1:4)
    {
        mu <- n[j]*p[i]
        sd <- sqrt(mu * (1-p[i]))
        lo <- round(mu-3*sd)
        hi <- round(mu+3*sd)
        if (hi-lo<40)
        x <- seq(lo,hi,by=1) else
        x <- round(seq(lo,hi,len=40))
        pdf <- cbind(dbinom(x,n[j],p[i]),dpois(x,mu),dnorm(x,mu,sd))
        pdf[x<0,1:2]<- 0
        matplot(x,pdf,main=paste("p=", p[i],", n=",n[j],sep=""))
    }
}

estim.beta  <- function(k)
{
delta1      <- rep(NA, k)
delta2      <- rep(NA, k)
for (i in 1:k)
{
x   <- rbeta(500, shape1=2.5, shape2=5)
delta1[i] <- mean(as.numeric(I(0.2 <x & x<0.4)))
delta2[i] <- mean(sin(x)*exp(-x))
}
return(cbind(delta1,delta2))
}

delta.estim <- estim.beta(500)
apply(delta.estim,2,mean)
apply(delta.estim,2,var)
sqrt(apply(delta.estim,2,var))

buf<-function(n, d,l)
{

       R             <- rep(NA, n)
       x             <- runif(n, 0, d/2)
       theta         <- runif(n, 0, pi/2)
       y             <- (l/2)*cos(theta)
       R             <- ifelse(y > x,1,0)
       pi            <- cumsum(R)/(1:n)
       rho           <- l/d
       phi           <- pi/(2*rho)
       pi.hat        <- 1/phi
       pi.hat
       plot(1:n, pi.hat, type="l", xlab="Number of Simulations", ylab="Proportion of Hits")
}


theta.u <- .8
theta.i <- .4
prop.u <- .3
prop.i <- 1 - prop.u
theta <- prop.u * theta.u + prop.i * theta.i

sim.1 <- function() {
x <- rbinom(1,sampsize,theta)
return ( x / sampsize )
}

sim.2 <- function() {
n.u <- rbinom ( 1, sampsize, prop.u )
n.i <- sampsize - n.u
x.u <- rbinom ( 1, n.u, theta.u )
x.i <- rbinom ( 1, n.i, theta.i )
t.hat.u <- x.u / n.u
t.hat.i <- x.i / n.i
return ( prop.u * t.hat.u + (1-prop.u) * t.hat.i )
}

sim.3 <- function() {
n.u <- sampsize * prop.u
n.i <- sampsize * prop.i
x.u <- rbinom ( 1, n.u, theta.u )
x.i <- rbinom ( 1, n.i, theta.i )
t.hat.u <- x.u / n.u
t.hat.i <- x.i / n.i
return ( prop.u * t.hat.u + (1-prop.u) * t.hat.i )
}

sampsize <- 100
n.times <- 1000 
theta.hat <- matrix ( NA, n.times, 3 )
for ( i in 1:n.times ) {
theta.hat[i,1] <- sim.1()
theta.hat[i,2] <- sim.2()
theta.hat[i,3] <- sim.3()
}

print ( apply(theta.hat,2,mean) )
boxplot ( theta.hat ~ col(theta.hat) )