We can simulate what can be a course like us

epilepsy <- function(size) {
    return(sample(c("Yes","No"), size=size, 
            replace = TRUE, prob = c(0.01, 0.99)))
}

Here we fix the real probability to 1%

As we know, R has four basic data type: numeric, factor, logic and character. If we mix them, they can change. For example

c(TRUE, FALSE)

[1]  TRUE FALSE

c(TRUE, FALSE, 2)

[1] 1 0 2

Mixing logic and numeric gives numbers

Since TRUE becomes 1 and FALSE becomes 0, it is easy to count how many times the result was TRUE

sum(courses[,723]=="Yes")

[1] 2

Now we just have to do this for each column of courses

(that is, for each simulation of a course)

We see near 30% of courses have nobody with epilepsy

Only one third of times we see one person in 97 with epilepsy

How can this be 1% of world population? Did we do something wrong?

mean(courses=="Yes")

[1] 0.01006186

No problem here. Overall population is ~1%

These two problems are related:

1. From reality to theory: we do experiments to learn about the population.

   In particular about the probabilities of each outcome. That is, the probability distribution

2. From theory to reality: assuming some probabilities, we can use the computer to explore what can be the consequences.

   If consequences do not match reality, we have to change our hypothesis

cases[723]>0

or directly

sum(courses[,723]=="Yes")>0

Even more easy and clear

any(courses[,723]=="Yes")

We store the official result in the vector official

We store our bet into our_bet

If official and bet are sorted we do

all(official==our_bet)

s <- sum_two_dice()
s==2 | s==4 | s==6 | s==8 | s==10 | s==12

It is better if we use reminder

s %% 2 == 0

We normally represent events with a function

event <- function(outcome) {...}

Then we can calculate the empirical frequency of an event

result <- replicate(N, event(simulation(...))
freq <- sum(result)/N 

We have to choose N carefully

We will test different probabilities. Let’s start with 0.5.

p <- 0.5

We can simulate n experiments with success probability p

empirical_freq <- function(n, p) {
  experiments <- sample(c(TRUE,FALSE), prob = c(p, 1-p),
                        size=n, replace = TRUE)
  return(sum(experiments)/n)
}

• If the real probability of an event is \(p\)
• and the sample size is \(n\)
• we calculate \(\sigma=\sqrt{p(1-p)/n}\)
• then the empirical frequency \(m/n\) is in the range
  • \(p\pm 2\sigma\) with probability 95%
  • \(p\pm 3\sigma\) with probability 99%
  • \(p\pm 5\sigma\) with probability 99.999%

The result depends on a value \(z\). Let \[\begin{aligned}\tilde{n} &= n + z^2\\ {\tilde {p}}&=\frac{m+z^2/2}{\tilde{n}}=\frac{m+z^2/2}{n+z^2}\end{aligned}\] Then \({\tilde {\sigma}(z)}={\sqrt{\tilde {p}\left(1-{\tilde {p}}\right)/\tilde {n}}}\), and the real probability \(p\) is in the range

• \(\tilde {p}\pm 2\cdot\tilde {\sigma}(2)\) with probability 95%
• \(\tilde {p}\pm 3\cdot\tilde {\sigma}(3)\) with probability 99%
• \(\tilde {p}\pm 5\cdot\tilde {\sigma}(5)\) with probability 99.999%

Let \(z=2\). Then \(\tilde{n} = n + 4\) and \({\tilde {p}}=(m+2)/(n + 4).\)

The confidence interval for \(p\) is given by \[\frac{m+2}{n + 4}\pm 2{\sqrt{\frac{m+2}{(n + 4)^2}\left(1-\frac{m+2}{n + 4}\right)}}\]

Real probability is in this range with 95% confidence.

This is the formula we have to remember

• Points are your experimental result (event frequency)
• Vertical lines represent the calculated 95% interval
• Horizontal line represent the real probability

• Event: function from outcome to logic
• You can use any(...)
• You can use all(...)
• Event count: sum(replicate(N, event(...)))
• Event frequency: sum(replicate(N, event(...)))/N
• Real probability is somewhere in \(\tilde {p}\pm 2\cdot\tilde {\sigma}(2)\) with probability 95%

Agresti, Alan; Coull, Brent A. (1998). “Approximate is better than ‘exact’ for interval estimation of binomial proportions”. The American Statistician. 52: 119–126. doi:10.2307/2685469. JSTOR 2685469. MR 1628435.