## box model expected value and standard error simulation
## ADM 10/3/2007

## specifics
the.box <- c(0,2,3,4,6)
draws <- 25

## perform simulation
G <- 10000
holder <- rep(NA, G)

for(g in 1:G) {
   holder[g] <- sum(sample(the.box, draws, replace=TRUE))
}

## plot the results
hist(holder, freq=FALSE, breaks=50)
plot(density(holder))

## estimate the expected value from the simulation
mean(holder)

## estimate the standard error from the simulation
sd(holder)

## estimate the probability of the sum between 50 and 100
sum(holder >= 50 & holder <= 100) / G

## estimate the probability of the sum between 80 and 90
sum(holder >= 80 & holder <= 90) / G

## same computations based on analytic results using the Normal distribution
pnorm(100, mean=75, sd=10) - pnorm(50, mean=75, sd=10)
pnorm(90, mean=75, sd=10) - pnorm(80, mean=75, sd=10)