# Here's an easy way to sample a Bernoulli trial if the
# probability of success is 1/2:

> sample(0:1, 1)
[1] 1
> sample(0:1, 1)
[1] 0
> sample(0:1, 1)
[1] 1
> sample(0:1, 1)
[1] 1
> sample(0:1, 1)
[1] 0

# Of course, if we want more values, we'll have to sample
# with replacement:

> sample(0:1, 10)
Error in sample.int(length(x), size, replace, prob) : 
  cannot take a sample larger than the population when 'replace = FALSE'
> sample(0:1, 10, replace=TRUE)
 [1] 0 0 0 1 0 1 0 0 0 1

# Here, we take a large sample of Bernoulli trials to assess
# via simulation the things we had said about the probability
# distribution based on just thinking about it:

> sample(0:1, 100000, replace=TRUE) -> CoinTosses
> head(CoinTosses)
[1] 1 0 0 0 0 1
> table(CoinTosses)
CoinTosses
    0     1 
50190 49810 
> barplot(table(CoinTosses))
> abline(h=0)

# We argued that the mean of the probability distribution
# should be 1/2:

> mean(CoinTosses)
[1] 0.4981

# Same for the standard deviation:

> sd(CoinTosses)
[1] 0.4999989

# Our Pearson skew index won't work because order statistics
# don't work well with discrete variables:

> pskew(CoinTosses)
[1] 2.988607

# But the third moment skew statistic shows that the
# distribution is symmetric:

> mskew(CoinTosses)
[1] 0.007600169

# The iqe won't work because, again, it is based on 
# order statistics:

> iqr(CoinTosses)
[1] 1

# Here, we go through a similar exercise for the uniform
# distribution over the interval (0,1):
 
> Spinner <- runif(100000)
> head(Spinner)
[1] 0.37189407 0.60158209 0.25346424 0.90511831 0.03558984 0.45839622
> hist(Spinner)
> hist(Spinner, freq=FALSE)

# The mean, as predicted, is 1/2:

> mean(Spinner)
[1] 0.4988157

# The standard deviation is the square root of 1/12:

> sd(Spinner)
[1] 0.2887953
> sqrt(1/12)
[1] 0.2886751

# The median and iqr are both 1/2:

> median(Spinner)
[1] 0.4971595
> iqr(Spinner)
[1] 0.5009139

# Skew indices indicate symmetry:

> pskew(Spinner)
[1] 0.01720422
> mskew(Spinner)
[1] 0.009131002
>