Home » Software Development » R: Bootstrap confidence intervals

### About Mark Needham # R: Bootstrap confidence intervals

I recently came across an interesting post on Julia Evans’ blog showing how to generate a bigger set of data points by sampling the small set of data points that we actually have using bootstrapping. Julia’s examples are all in Python so I thought it’d be a fun exercise to translate them into R.

We’re doing the bootstrapping to simulate the number of no-shows for a flight so we can work out how many seats we can overbook the plane by.

We start out with a small sample of no-shows and work off the assumption that it’s ok to kick someone off a flight 5% of the time. Let’s work out how many people that’d be for our initial sample:

```> data = c(0, 1, 3, 2, 8, 2, 3, 4)
> quantile(data, 0.05)
5%
0.35```

0.35 people! That’s not a particularly useful result so we’re going to resample the initial data set 10,000 times, taking the 5%ile each time and see if we come up with something better:

We’re going to use the sample function with replacement to generate our resamples:

```> sample(data, replace = TRUE)
 0 3 2 8 8 0 8 0
> sample(data, replace = TRUE)
 2 2 4 3 4 4 2 2```

Now let’s write a function to do that multiple times:

```library(ggplot)

bootstrap_5th_percentile = function(data, n_bootstraps) {
return(sapply(1:n_bootstraps,
function(iteration) quantile(sample(data, replace = TRUE), 0.05)))
}

values = bootstrap_5th_percentile(data, 10000)

ggplot(aes(x = value), data = data.frame(value = values)) + geom_histogram(binwidth=0.25)``` So this visualisation is telling us that we can oversell by 0-2 people but we don’t know an exact number.

Let’s try the same exercise but with a bigger initial data set of 1,000 values rather than just 8. First we’ll generate a distribution (with a mean of 5 and standard deviation of 2) and visualise it:

```library(dplyr)

df = data.frame(value = rnorm(1000,5, 2))
df = df %>% filter(value >= 0) %>% mutate(value = as.integer(round(value)))
ggplot(aes(x = value), data = df) + geom_histogram(binwidth=1)``` Our distribution seems to have a lot more values around 4 & 5 whereas the Python version has a flatter distribution – I’m not sure why that is so if you have any ideas let me know. In any case, let’s check the 5%ile for this data set:

```> quantile(df\$value, 0.05)
5%
2```

Cool! Now at least we have an integer value rather than the 0.35 we got earlier. Finally let’s do some bootstrapping over our new distribution and see what 5%ile we come up with:

```resampled = bootstrap_5th_percentile(df\$value, 10000)
byValue = data.frame(value = resampled) %>% count(value)

> byValue
Source: local data frame [3 x 2]

value    n
1   1.0    3
2   1.7    2
3   2.0 9995

ggplot(aes(x = value, y = n), data = byValue) + geom_bar(stat = "identity")``` ‘2’ is by far the most popular 5%ile here although it seems weighted more towards that value than with Julia’s Python version, which I imagine is because we seem to have sampled from a slightly different distribution.

 Reference: R: Bootstrap confidence intervals from our JCG partner Mark Needham at the Mark Needham Blog blog.

# Do you want to know how to develop your skillset to become a Java Rockstar?

## Subscribe to our newsletter to start Rocking right now!

### and many more ....

Subscribe
Notify of This site uses Akismet to reduce spam. Learn how your comment data is processed.

1 Comment
Inline Feedbacks 